More Robust Integrator cpp

Question

I'm trying to produce a relativistic voigt distribution, the convolution of a relativistic Breit-Wigner distribution and a gaussian function.

MWE:

double relativisticBreitwigner_pdf(double energy, double width, double mass){
  double massSquare = pow(mass,2);
  double widthSquare = pow(width,2);
  double gamma = sqrt(massSquare*(massSquare+widthSquare));
  double k = (2*sqrt(2)*mass*width*gamma)/(M_PI*sqrt(massSquare + gamma) );
  return k/(pow((pow(energy,2)-massSquare),2) + massSquare*widthSquare);
}

double gaussian_pdf(double energy, double sigma, double mass){
  return (1.0/(sigma*sqrt(2*M_PI)))*exp(-(1.0/2.0)*pow((energy-mass)/sigma,2.0));
}

double relativisticVoigt_pdf(double energy, double width, double mass, double sigma, double range=100.0){

  auto f = [&](double dummy) { return ( relativisticBreitwigner_pdf(dummy+mass,width,mass)*gaussian_pdf(energy-dummy,sigma,mass) );};
  boost::math::quadrature::tanh_sinh<double> integrator;
  return integrator.integrate(f,-range,range);
  //  return boost::math::quadrature::trapezoidal(f,-range,range,sqrt(std::numeric_limits<double>::epsilon()),10000);    }

This function relativisticVoigt_pdf(...) correctly produces a relativistic voigt distribution, however there are many 'dips' in the distribution which are incorrect which are due to the value returned from integrator.integrate(f,-range,range); not being correct.

If I reduce the integration range the size/number of these dips is smaller, but then the range of the relativistic voigt distribution is cutoff.

Attached a screenshot showing the relativistic voigt in black with that problem (compared to a nonrelativistic voigt in pink that doesn't have this problem just to check that the values outside of the dips are valid, the black curve should be close to the pink curve but slightly above it to the left of the peak and slightly below it to the right of the peak, as can be seen happens).

I assume the problem is to do with rounding errors in the integration method due to the many small numbers involved, but this is just a guess. Is there a more robust/reliable integrator that works well in this regime?

tanh_sinh integration range=100

trapezoidal integration range=250

Answer 1

For anyone in the future that has a similar problem, by dropping the tolerance of boost::math::quadrature::trapezoidal(...) from the default (sqrt(std::numeric_limits::epsilon()=1.48e-8, which in my code above I put in explicitly but this is the default value if it is not put in) to a larger value 1e-6 I manage to get the trapezoidal integrator to work over the full range.

I do not understand why this works, particularly since the dips do not all go to zero some of them just go a bit below the actual value. If they all went to zero I would understand it as perhaps it just not finding the integral to the required precision and hence returning zero. If anyone understands why having a more precise tolerance results in the integral being wrong I'd like to know.

Update:

While the above helped, it did not remove the problem. To solve the problem you can take the L1 parameter that is returned by trapezoidal as below:

double error; 
double L1; 
boost::math::quadrature::trapezoidal(f,-range,range,1e-6,10000,&error,&L1);

Then check if L1==0 or very small, if it is, vary the range until it isn't.