NIntegrate blowing up/behaving weirdly at the turning point of the integrand

I’m performing (what should be) a straightforward numerical integration (Fourier transform) of a function with no poles / singularities (at least in a particular parameter regime):

<< FourierSeries`

ClearAll["Global`*"]

l = 1;
ϵ = 10^-4;
r = 0.1; (*radial coordinate*)

p = 0;
P = 80;

Data1 = ConstantArray[Null, {P, 2}];

While[p < P, p++;
  
  w = -10 + 20 p/P;
  
  σ = (
   2 (-2 l^2 r^2 Sin[s/(2 Sqrt[l^-2 - r^2])]^2 + (1 - l^2 r^2) 2 Sinh[
        s/(2 Sqrt[l^-2 - r^2]) - I ϵ]^2))/l^2;
  
  RF = NInverseFourierTransform[-1/(4 π^2)/σ, s,
      w, Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
     MinRecursion -> 6] // Chop;
  
  Data1[[p, 1]] = w;
  Data1[[p, 2]] = RF;
  
  ];

ListLinePlot[Data1, PlotRange -> All, 
   LabelStyle -> {FontSize -> 16, FontFamily -> "CMU Serif", Black}, 
   PlotLabel -> StringForm["r = ``.", {r // N}]] // Print;

When the parameter r is small, the numerics seem to work fine/as expected, producing plots like this (r = 0.1):

(x-axis is the energy, y-axis is the transition rate of a two-level system – this shows a roughly thermal Planck spectrum).

When increasing r beyond some critical point, the numerics seem to blow up, giving a weird answer. For example, r = 0.8 yields:

(see especially the magnitude of the y-axis).

I plotted the integrand for these two values below. The numerics seem to blow up when the function bifurcates from having one turning point to two turning points:

r=0.1 (integrand plotted as a function of s)

and r =0.8

Why does NIntegrate not like this seemingly inconspicuous change in behaviour of the integrand? Any help would be greatly appreciated! (If there’s an analytic solution, even better :P)

AbsoluteTiming[
 Table[
  Block[{},
   l = 1;
   [Epsilon] = 10^-4;
   p = 0;
   P = 80;
   
   Data1 = ConstantArray[Null, {P, 2}];
   
   While[p < P,
    p++;
    w = -10 + 20 p/P;
    [Sigma] = (2 (-2 l^2 r^2 Sin[s/(2 Sqrt[l^-2 - r^2])]^2 + (1 - 
             l^2 r^2) 2 Sinh[
             s/(2 Sqrt[l^-2 - r^2]) - I [Epsilon]]^2))/l^2;
    RF = NInverseFourierTransform[-1/(4 [Pi]^2) E^(-I w s)/[Sigma], 
      s, w, Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0, 
        "SingularityHandler" -> None, "MaxErrorIncreases" -> 100000}, 
      MinRecursion -> 6, MaxRecursion -> 100, WorkingPrecision -> 30, 
      PrecisionGoal -> 4, AccuracyGoal -> 4];
    Data1[[p, 1]] = w;
    Data1[[p, 2]] = RF;
    ];
   ListLinePlot[Data1, PlotRange -> All, 
    LabelStyle -> {FontSize -> 16, FontFamily -> "CMU Serif", Black}, 
    PlotLabel -> StringForm["r = ``.", {r // N}]]
   ],
  {r, Range[7, 8, 1/10]/10}]
 ]