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

1. int to double conversion could be an issue in first function. 2. There's others constants other than `M_PI`, you probably should use instead of calculating it yourself. 3. Add more formatting to your code please. It's ok to write oneliners but they're hard to debug. 4. What are the dip regions input? It's probably good idea to manually calculate/debug what should happen there — Sugar, Feb 11 '21 at 16:52
Did you try other integration methods? Even a basic one? I would be surprised that these dips are due to rounding errors. — Damien, Feb 11 '21 at 16:54
and, yes,as @Damien said, did you tried check precision errors accumulation and overflow/nan values after certain steps? bigger values could produce bigger errors you see on dips — Sugar, Feb 11 '21 at 17:01
Hi, I do not understand what problem with int to double conversion you are referring to, or what constants other than M_PI there are? Either way there is no issue with the 1st or 2nd function, they are correct everywhere. — John Cunningham, Feb 11 '21 at 17:21
I've tried both the integration methods that I have in the MWE (with one commented out) and they both have similar issues, I've attached a picture of the trapezoidal case for range=250 but the dips are in different places depending on the integrator and the integration range chosen The dip at 544 for example with tanh_sinh and range=100 is relativisticVoigt_pdf(544,2.085,501,1) = 3.48874e-06 — John Cunningham, Feb 11 '21 at 17:23
The pink line is from TMath::Voigt((energy-mass),sigma,width) — John Cunningham, Feb 11 '21 at 17:23

John Cunningham · Answer 1 · 2021-02-14T06:03:48.030

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.

Producing more distributions I see this solution just makes the problem more rare, it still occurs for certain values. Is there any solution than properly solves this? — John Cunningham, Feb 14 '21 at 02:51

More Robust Integrator cpp

1 Answers1