18 static const double kSqrt2 = 1.41421356237309515;
20 double beta_cdf_c(
double x,
double a,
double b)
23 return ROOT::Math::inc_beta(1-x, b, a);
27 double beta_cdf(
double x,
double a,
double b )
29 return ROOT::Math::inc_beta(x, a, b);
33 double breitwigner_cdf_c(
double x,
double gamma,
double x0)
35 return 0.5 - std::atan(2.0 * (x-x0) / gamma) / M_PI;
39 double breitwigner_cdf(
double x,
double gamma,
double x0)
41 return 0.5 + std::atan(2.0 * (x-x0) / gamma) / M_PI;
45 double cauchy_cdf_c(
double x,
double b,
double x0)
47 return 0.5 - std::atan( (x-x0) / b) / M_PI;
51 double cauchy_cdf(
double x,
double b,
double x0)
53 return 0.5 + std::atan( (x-x0) / b) / M_PI;
57 double chisquared_cdf_c(
double x,
double r,
double x0)
59 return ROOT::Math::inc_gamma_c ( 0.5 * r , 0.5* (x-x0) );
63 double chisquared_cdf(
double x,
double r,
double x0)
65 return ROOT::Math::inc_gamma ( 0.5 * r , 0.5* (x-x0) );
69 double crystalball_cdf(
double x,
double alpha,
double n,
double sigma,
double mean )
72 MATH_ERROR_MSG(
"crystalball_cdf",
"CrystalBall cdf not defined for n <=1");
73 return std::numeric_limits<double>::quiet_NaN();
76 double abs_alpha = std::abs(alpha);
77 double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
78 double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
79 double totIntegral = sigma*(C+D);
81 double integral = crystalball_integral(x,alpha,n,sigma,mean);
82 return (alpha > 0) ? 1. - integral/totIntegral : integral/totIntegral;
84 double crystalball_cdf_c(
double x,
double alpha,
double n,
double sigma,
double mean )
87 MATH_ERROR_MSG(
"crystalball_cdf_c",
"CrystalBall cdf not defined for n <=1");
88 return std::numeric_limits<double>::quiet_NaN();
90 double abs_alpha = std::abs(alpha);
91 double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
92 double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
93 double totIntegral = sigma*(C+D);
95 double integral = crystalball_integral(x,alpha,n,sigma,mean);
96 return (alpha > 0) ? integral/totIntegral : 1. - (integral/totIntegral);
98 double crystalball_integral(
double x,
double alpha,
double n,
double sigma,
double mean)
107 if (sigma == 0)
return 0;
110 MATH_ERROR_MSG(
"crystalball_integral",
"CrystalBall function not defined at alpha=0");
113 bool useLog = (n == 1.0);
114 if (n<=0) MATH_WARN_MSG(
"crystalball_integral",
"No physical meaning when n<=0");
116 double z = (x-mean)/sigma;
117 if (alpha < 0 ) z = -z;
119 double abs_alpha = std::abs(alpha);
126 const double sqrtpiover2 = std::sqrt(M_PI/2.);
127 const double sqrt2pi = std::sqrt( 2.*M_PI);
128 const double oneoversqrt2 = 1./sqrt(2.);
131 double A = std::pow(n/abs_alpha,n) * std::exp(-0.5 * alpha*alpha);
132 double B = n/abs_alpha - abs_alpha;
135 double C = (n/abs_alpha) * (1./(n-1)) * std::exp(-alpha*alpha/2.);
136 intpow = C - A /(n-1.) * std::pow(B-z,-n+1) ;
140 intpow = -A * std::log( n / abs_alpha ) + A * std::log( B -z );
142 intgaus = sqrtpiover2*(1.+ROOT::Math::erf(abs_alpha*oneoversqrt2));
146 intgaus = ROOT::Math::gaussian_cdf_c(z, 1);
150 return sigma * (intgaus + intpow);
154 double exponential_cdf_c(
double x,
double lambda,
double x0)
156 if ((x-x0) < 0)
return 1.0;
157 else return std::exp(- lambda * (x-x0));
161 double exponential_cdf(
double x,
double lambda,
double x0)
163 if ((x-x0) < 0)
return 0.0;
165 return - ROOT::Math::expm1( - lambda * (x-x0) ) ;
169 double fdistribution_cdf_c(
double x,
double n,
double m,
double x0)
172 if (n < 0 || m < 0)
return std::numeric_limits<double>::quiet_NaN();
174 double z = m/(m + n*(x-x0));
176 if (z > 0.9 && n > 1 && m > 1)
return 1.- fdistribution_cdf(x,n,m,x0);
179 return ROOT::Math::inc_beta(m/(m + n*(x-x0)), .5*m, .5*n);
183 double fdistribution_cdf(
double x,
double n,
double m,
double x0)
187 return std::numeric_limits<double>::quiet_NaN();
189 double z = n*(x-x0)/(m + n*(x-x0));
191 if (z > 0.9 && n > 1 && m > 1)
192 return 1. - fdistribution_cdf_c(x,n,m,x0);
194 return ROOT::Math::inc_beta(z, .5*n, .5*m);
198 double gamma_cdf_c(
double x,
double alpha,
double theta,
double x0)
200 return ROOT::Math::inc_gamma_c(alpha, (x-x0)/theta);
204 double gamma_cdf(
double x,
double alpha,
double theta,
double x0)
206 return ROOT::Math::inc_gamma(alpha, (x-x0)/theta);
210 double lognormal_cdf_c(
double x,
double m,
double s,
double x0)
212 double z = (std::log((x-x0))-m)/(s*kSqrt2);
213 if (z > 1.)
return 0.5*ROOT::Math::erfc(z);
214 else return 0.5*(1.0 - ROOT::Math::erf(z));
218 double lognormal_cdf(
double x,
double m,
double s,
double x0)
220 double z = (std::log((x-x0))-m)/(s*kSqrt2);
221 if (z < -1.)
return 0.5*ROOT::Math::erfc(-z);
222 else return 0.5*(1.0 + ROOT::Math::erf(z));
226 double normal_cdf_c(
double x,
double sigma,
double x0)
228 double z = (x-x0)/(sigma*kSqrt2);
229 if (z > 1.)
return 0.5*ROOT::Math::erfc(z);
230 else return 0.5*(1.-ROOT::Math::erf(z));
234 double normal_cdf(
double x,
double sigma,
double x0)
236 double z = (x-x0)/(sigma*kSqrt2);
237 if (z < -1.)
return 0.5*ROOT::Math::erfc(-z);
238 else return 0.5*(1.0 + ROOT::Math::erf(z));
242 double tdistribution_cdf_c(
double x,
double r,
double x0)
245 double sign = (p>0) ? 1. : -1;
246 return .5 - .5*ROOT::Math::inc_beta(p*p/(r + p*p), .5, .5*r)*sign;
250 double tdistribution_cdf(
double x,
double r,
double x0)
253 double sign = (p>0) ? 1. : -1;
254 return .5 + .5*ROOT::Math::inc_beta(p*p/(r + p*p), .5, .5*r)*sign;
258 double uniform_cdf_c(
double x,
double a,
double b,
double x0)
260 if ((x-x0) < a)
return 1.0;
261 else if ((x-x0) >= b)
return 0.0;
262 else return (b-(x-x0))/(b-a);
266 double uniform_cdf(
double x,
double a,
double b,
double x0)
268 if ((x-x0) < a)
return 0.0;
269 else if ((x-x0) >= b)
return 1.0;
270 else return ((x-x0)-a)/(b-a);
275 double poisson_cdf_c(
unsigned int n,
double mu)
279 double a = (double) n + 1.0;
280 return ROOT::Math::gamma_cdf(mu, a, 1.0);
284 double poisson_cdf(
unsigned int n,
double mu)
288 double a = (double) n + 1.0;
289 return ROOT::Math::gamma_cdf_c(mu, a, 1.0);
293 double binomial_cdf_c(
unsigned int k,
double p,
unsigned int n)
297 if ( k >= n)
return 0;
298 double a = (double) k + 1.0;
299 double b = (double) n - k;
300 return ROOT::Math::beta_cdf(p, a, b);
304 double binomial_cdf(
unsigned int k,
double p,
unsigned int n)
308 if ( k >= n)
return 1.0;
310 double a = (double) k + 1.0;
311 double b = (double) n - k;
312 return ROOT::Math::beta_cdf_c(p, a, b);
316 double negative_binomial_cdf(
unsigned int k,
double p,
double n)
321 if (p < 0 || p > 1)
return 0;
322 return ROOT::Math::beta_cdf(p, n, k+1.0);
326 double negative_binomial_cdf_c(
unsigned int k,
double p,
double n)
330 if ( n < 0)
return 0;
331 if ( p < 0 || p > 1)
return 0;
332 return ROOT::Math::beta_cdf_c(p, n, k+1.0);
336 double landau_cdf(
double x,
double xi,
double x0)
343 static double p1[5] = {0.2514091491e+0,-0.6250580444e-1, 0.1458381230e-1,-0.2108817737e-2, 0.7411247290e-3};
344 static double q1[5] = {1.0 ,-0.5571175625e-2, 0.6225310236e-1,-0.3137378427e-2, 0.1931496439e-2};
346 static double p2[4] = {0.2868328584e+0, 0.3564363231e+0, 0.1523518695e+0, 0.2251304883e-1};
347 static double q2[4] = {1.0 , 0.6191136137e+0, 0.1720721448e+0, 0.2278594771e-1};
349 static double p3[4] = {0.2868329066e+0, 0.3003828436e+0, 0.9950951941e-1, 0.8733827185e-2};
350 static double q3[4] = {1.0 , 0.4237190502e+0, 0.1095631512e+0, 0.8693851567e-2};
352 static double p4[4] = {0.1000351630e+1, 0.4503592498e+1, 0.1085883880e+2, 0.7536052269e+1};
353 static double q4[4] = {1.0 , 0.5539969678e+1, 0.1933581111e+2, 0.2721321508e+2};
355 static double p5[4] = {0.1000006517e+1, 0.4909414111e+2, 0.8505544753e+2, 0.1532153455e+3};
356 static double q5[4] = {1.0 , 0.5009928881e+2, 0.1399819104e+3, 0.4200002909e+3};
358 static double p6[4] = {0.1000000983e+1, 0.1329868456e+3, 0.9162149244e+3,-0.9605054274e+3};
359 static double q6[4] = {1.0 , 0.1339887843e+3, 0.1055990413e+4, 0.5532224619e+3};
361 static double a1[4] = {0 ,-0.4583333333e+0, 0.6675347222e+0,-0.1641741416e+1};
362 static double a2[4] = {0 , 1.0 ,-0.4227843351e+0,-0.2043403138e+1};
364 double v = (x - x0)/xi;
371 lan = 0.3989422803*std::exp(-1./u)*std::sqrt(u)*(1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
376 lan = (std::exp(-u)/std::sqrt(u))*(p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
377 (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
380 lan = (p2[0]+(p2[1]+(p2[2]+p2[3]*v)*v)*v)/(q2[0]+(q2[1]+(q2[2]+q2[3]*v)*v)*v);
383 lan = (p3[0]+(p3[1]+(p3[2]+p3[3]*v)*v)*v)/(q3[0]+(q3[1]+(q3[2]+q3[3]*v)*v)*v);
388 lan = (p4[0]+(p4[1]+(p4[2]+p4[3]*u)*u)*u)/(q4[0]+(q4[1]+(q4[2]+q4[3]*u)*u)*u);
393 lan = (p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/(q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
398 lan = (p6[0]+(p6[1]+(p6[2]+p6[3]*u)*u)*u)/(q6[0]+(q6[1]+(q6[2]+q6[3]*u)*u)*u);
402 u = 1./(v-v*std::log(v)/(v+1));
403 lan = 1-(a2[1]+(a2[2]+a2[3]*u)*u)*u;
409 double landau_xm1(
double x,
double xi,
double x0)
414 static double p1[5] = {-0.8949374280E+0, 0.4631783434E+0,-0.4053332915E-1,
415 0.1580075560E-1,-0.3423874194E-2};
416 static double q1[5] = { 1.0 , 0.1002930749E+0, 0.3575271633E-1,
417 -0.1915882099E-2, 0.4811072364E-4};
418 static double p2[5] = {-0.8933384046E+0, 0.1161296496E+0, 0.1200082940E+0,
419 0.2185699725E-1, 0.2128892058E-2};
420 static double q2[5] = { 1.0 , 0.4935531886E+0, 0.1066347067E+0,
421 0.1250161833E-1, 0.5494243254E-3};
422 static double p3[5] = {-0.8933322067E+0, 0.2339544896E+0, 0.8257653222E-1,
423 0.1411226998E-1, 0.2892240953E-3};
424 static double q3[5] = { 1.0 , 0.3616538408E+0, 0.6628026743E-1,
425 0.4839298984E-2, 0.5248310361E-4};
426 static double p4[4] = { 0.9358419425E+0, 0.6716831438E+2,-0.6765069077E+3,
428 static double q4[4] = { 1.0 , 0.7752562854E+2,-0.5637811998E+3,
430 static double p5[4] = { 0.9489335583E+0, 0.5561246706E+3, 0.3208274617E+5,
432 static double q5[4] = { 1.0 , 0.6028275940E+3, 0.3716962017E+5,
434 static double a0[6] = {-0.4227843351E+0,-0.1544313298E+0, 0.4227843351E+0,
435 0.3276496874E+1, 0.2043403138E+1,-0.8681296500E+1};
436 static double a1[4] = { 0, -0.4583333333E+0, 0.6675347222E+0,
438 static double a2[5] = { 0, -0.1958333333E+1, 0.5563368056E+1,
439 -0.2111352961E+2, 0.1006946266E+3};
441 double v = (x-x0)/xi;
445 double u = std::exp(v+1);
446 xm1lan = v-u*(1+(a2[1]+(a2[2]+(a2[3]+a2[4]*u)*u)*u)*u)/
447 (1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
451 xm1lan = (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
452 (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
456 xm1lan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
457 (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
461 xm1lan = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
462 (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
467 xm1lan = std::log(v)*(p4[0]+(p4[1]+(p4[2]+p4[3]*u)*u)*u)/
468 (q4[0]+(q4[1]+(q4[2]+q4[3]*u)*u)*u);
473 xm1lan = std::log(v)*(p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/
474 (q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
478 double u = v-v*std::log(v)/(v+1);
479 v = 1/(u-u*(u+ std::log(u)-v)/(u+1));
481 xm1lan = (u+a0[0]+(-u+a0[1]+(a0[2]*u+a0[3]+(a0[4]*u+a0[5])*v)*v)*v)/
482 (1-(1-(a0[2]+a0[4]*v)*v)*v);
484 return xm1lan*xi + x0;
489 double landau_xm2(
double x,
double xi,
double x0)
494 static double p1[5] = { 0.1169837582E+1,-0.4834874539E+0, 0.4383774644E+0,
495 0.3287175228E-2, 0.1879129206E-1};
496 static double q1[5] = { 1.0 , 0.1795154326E+0, 0.4612795899E-1,
497 0.2183459337E-2, 0.7226623623E-4};
498 static double p2[5] = { 0.1157939823E+1,-0.3842809495E+0, 0.3317532899E+0,
499 0.3547606781E-1, 0.6725645279E-2};
500 static double q2[5] = { 1.0 , 0.2916824021E+0, 0.5259853480E-1,
501 0.3840011061E-2, 0.9950324173E-4};
502 static double p3[4] = { 0.1178191282E+1, 0.1011623342E+2,-0.1285585291E+2,
504 static double q3[4] = { 1.0 , 0.8614160194E+1, 0.3118929630E+2,
506 static double p4[5] = { 0.1030763698E+1, 0.1216758660E+3, 0.1637431386E+4,
507 -0.2171466507E+4, 0.7010168358E+4};
508 static double q4[5] = { 1.0 , 0.1022487911E+3, 0.1377646350E+4,
509 0.3699184961E+4, 0.4251315610E+4};
510 static double p5[4] = { 0.1010084827E+1, 0.3944224824E+3, 0.1773025353E+5,
512 static double q5[4] = { 1.0 , 0.3605950254E+3, 0.1392784158E+5,
514 static double a0[7] = {-0.2043403138E+1,-0.8455686702E+0,-0.3088626596E+0,
515 0.5821346754E+1, 0.4227843351E+0, 0.6552993748E+1,
517 static double a1[4] = { 0. ,-0.4583333333E+0, 0.6675347222E+0,
519 static double a2[4] = {-0.1958333333E+1, 0.5563368056E+1,-0.2111352961E+2,
521 static double a3[4] = {-1.0 , 0.4458333333E+1,-0.2116753472E+2,
524 double v = (x-x0)/xi;
528 double u = std::exp(v+1);
530 (v/u+a2[0]*v+a3[0]+(a2[1]*v+a3[1]+(a2[2]*v+a3[2]+
531 (a2[3]*v+a3[3])*u)*u)*u)/
532 (1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
536 xm2lan = (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
537 (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
541 xm2lan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
542 (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
547 xm2lan = v*(p3[0]+(p3[1]+(p3[2]+p3[3]*u)*u)*u)/
548 (q3[0]+(q3[1]+(q3[2]+q3[3]*u)*u)*u);
553 xm2lan = v*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
554 (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
559 xm2lan = v*(p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/
560 (q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
564 double u = v-v*std::log(v)/(v+1);
565 v = 1/(u-u*(u+log(u)-v)/(u+1));
567 xm2lan = (1/v+u*u+a0[0]+a0[1]*u+(-u*u+a0[2]*u+a0[3]+
568 (a0[4]*u*u+a0[5]*u+a0[6])*v)*v)/(1-(1-a0[4]*v)*v);
570 if (x0 == 0)
return xm2lan*xi*xi;
571 double xm1lan = ROOT::Math::landau_xm1(x, xi, x0);
572 return xm2lan*xi*xi + (2*xm1lan-x0)*x0;