Logo ROOT   6.19/01
Reference Guide
ProbFuncMathCore.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: L. Moneta, A. Zsenei 06/2005
3 
4 
5 
6 #include "Math/Math.h"
7 #include "Math/Error.h"
10 #include <stdio.h>
11 #include <limits>
12 using namespace std;
13 namespace ROOT {
14 namespace Math {
15 
16 
17 
18  static const double kSqrt2 = 1.41421356237309515; // sqrt(2.)
19 
20  double beta_cdf_c(double x, double a, double b)
21  {
22  // use the fact that I(x,a,b) = 1. - I(1-x,b,a)
23  return ROOT::Math::inc_beta(1-x, b, a);
24  }
25 
26 
27  double beta_cdf(double x, double a, double b )
28  {
29  return ROOT::Math::inc_beta(x, a, b);
30  }
31 
32 
33  double breitwigner_cdf_c(double x, double gamma, double x0)
34  {
35  return 0.5 - std::atan(2.0 * (x-x0) / gamma) / M_PI;
36  }
37 
38 
39  double breitwigner_cdf(double x, double gamma, double x0)
40  {
41  return 0.5 + std::atan(2.0 * (x-x0) / gamma) / M_PI;
42  }
43 
44 
45  double cauchy_cdf_c(double x, double b, double x0)
46  {
47  return 0.5 - std::atan( (x-x0) / b) / M_PI;
48  }
49 
50 
51  double cauchy_cdf(double x, double b, double x0)
52  {
53  return 0.5 + std::atan( (x-x0) / b) / M_PI;
54  }
55 
56 
57  double chisquared_cdf_c(double x, double r, double x0)
58  {
59  return ROOT::Math::inc_gamma_c ( 0.5 * r , 0.5* (x-x0) );
60  }
61 
62 
63  double chisquared_cdf(double x, double r, double x0)
64  {
65  return ROOT::Math::inc_gamma ( 0.5 * r , 0.5* (x-x0) );
66  }
67 
68 
69  double crystalball_cdf(double x, double alpha, double n, double sigma, double mean )
70  {
71  if (n <= 1.) {
72  MATH_ERROR_MSG("crystalball_cdf","CrystalBall cdf not defined for n <=1");
73  return std::numeric_limits<double>::quiet_NaN();
74  }
75 
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);
80 
81  double integral = crystalball_integral(x,alpha,n,sigma,mean);
82  return (alpha > 0) ? 1. - integral/totIntegral : integral/totIntegral;
83  }
84  double crystalball_cdf_c(double x, double alpha, double n, double sigma, double mean )
85  {
86  if (n <= 1.) {
87  MATH_ERROR_MSG("crystalball_cdf_c","CrystalBall cdf not defined for n <=1");
88  return std::numeric_limits<double>::quiet_NaN();
89  }
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);
94 
95  double integral = crystalball_integral(x,alpha,n,sigma,mean);
96  return (alpha > 0) ? integral/totIntegral : 1. - (integral/totIntegral);
97  }
98  double crystalball_integral(double x, double alpha, double n, double sigma, double mean)
99  {
100  // compute the integral of the crystal ball function (ROOT::Math::crystalball_function)
101  // If alpha > 0 the integral is the right tail integral.
102  // If alpha < 0 is the left tail integrals which are always finite for finite x.
103  // parameters:
104  // alpha : is non equal to zero, define the # of sigma from which it becomes a power-law function (from mean-alpha*sigma)
105  // n > 1 : is integrer, is the power of the low tail
106  // add a value xmin for cases when n <=1 the integral diverges
107  if (sigma == 0) return 0;
108  if (alpha==0)
109  {
110  MATH_ERROR_MSG("crystalball_integral","CrystalBall function not defined at alpha=0");
111  return 0.;
112  }
113  bool useLog = (n == 1.0);
114  if (n<=0) MATH_WARN_MSG("crystalball_integral","No physical meaning when n<=0");
115 
116  double z = (x-mean)/sigma;
117  if (alpha < 0 ) z = -z;
118 
119  double abs_alpha = std::abs(alpha);
120 
121  //double D = *(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
122  //double N = 1./(sigma*(C+D));
123  double intgaus = 0.;
124  double intpow = 0.;
125 
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.);
129  if (z <= -abs_alpha)
130  {
131  double A = std::pow(n/abs_alpha,n) * std::exp(-0.5 * alpha*alpha);
132  double B = n/abs_alpha - abs_alpha;
133 
134  if (!useLog) {
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) ;
137  }
138  else {
139  // for n=1 the primitive of 1/x is log(x)
140  intpow = -A * std::log( n / abs_alpha ) + A * std::log( B -z );
141  }
142  intgaus = sqrtpiover2*(1.+ROOT::Math::erf(abs_alpha*oneoversqrt2));
143  }
144  else
145  {
146  intgaus = ROOT::Math::gaussian_cdf_c(z, 1);
147  intgaus *= sqrt2pi;
148  intpow = 0;
149  }
150  return sigma * (intgaus + intpow);
151  }
152 
153 
154  double exponential_cdf_c(double x, double lambda, double x0)
155  {
156  if ((x-x0) < 0) return 1.0;
157  else return std::exp(- lambda * (x-x0));
158  }
159 
160 
161  double exponential_cdf(double x, double lambda, double x0)
162  {
163  if ((x-x0) < 0) return 0.0;
164  else // use expm1 function to avoid errors at small x
165  return - ROOT::Math::expm1( - lambda * (x-x0) ) ;
166  }
167 
168 
169  double fdistribution_cdf_c(double x, double n, double m, double x0)
170  {
171  // f distribution is defined only for both n and m > 0
172  if (n < 0 || m < 0) return std::numeric_limits<double>::quiet_NaN();
173 
174  double z = m/(m + n*(x-x0));
175  // fox z->1 and large a and b IB looses precision use complement function
176  if (z > 0.9 && n > 1 && m > 1) return 1.- fdistribution_cdf(x,n,m,x0);
177 
178  // for the complement use the fact that IB(x,a,b) = 1. - IB(1-x,b,a)
179  return ROOT::Math::inc_beta(m/(m + n*(x-x0)), .5*m, .5*n);
180  }
181 
182 
183  double fdistribution_cdf(double x, double n, double m, double x0)
184  {
185  // f distribution is defined only for both n and m > 0
186  if (n < 0 || m < 0)
187  return std::numeric_limits<double>::quiet_NaN();
188 
189  double z = n*(x-x0)/(m + n*(x-x0));
190  // fox z->1 and large a and b IB looses precision use complement function
191  if (z > 0.9 && n > 1 && m > 1)
192  return 1. - fdistribution_cdf_c(x,n,m,x0);
193 
194  return ROOT::Math::inc_beta(z, .5*n, .5*m);
195  }
196 
197 
198  double gamma_cdf_c(double x, double alpha, double theta, double x0)
199  {
200  return ROOT::Math::inc_gamma_c(alpha, (x-x0)/theta);
201  }
202 
203 
204  double gamma_cdf(double x, double alpha, double theta, double x0)
205  {
206  return ROOT::Math::inc_gamma(alpha, (x-x0)/theta);
207  }
208 
209 
210  double lognormal_cdf_c(double x, double m, double s, double x0)
211  {
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));
215  }
216 
217 
218  double lognormal_cdf(double x, double m, double s, double x0)
219  {
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));
223  }
224 
225 
226  double normal_cdf_c(double x, double sigma, double x0)
227  {
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));
231  }
232 
233 
234  double normal_cdf(double x, double sigma, double x0)
235  {
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));
239  }
240 
241 
242  double tdistribution_cdf_c(double x, double r, double x0)
243  {
244  double p = x-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;
247  }
248 
249 
250  double tdistribution_cdf(double x, double r, double x0)
251  {
252  double p = x-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;
255  }
256 
257 
258  double uniform_cdf_c(double x, double a, double b, double x0)
259  {
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);
263  }
264 
265 
266  double uniform_cdf(double x, double a, double b, double x0)
267  {
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);
271  }
272 
273  /// discrete distributions
274 
275  double poisson_cdf_c(unsigned int n, double mu)
276  {
277  // mu must be >= 0 . Use poisson - gamma relation
278  // Pr ( x <= n) = Pr( y >= a) where x is poisson and y is gamma distributed ( a = n+1)
279  double a = (double) n + 1.0;
280  return ROOT::Math::gamma_cdf(mu, a, 1.0);
281  }
282 
283 
284  double poisson_cdf(unsigned int n, double mu)
285  {
286  // mu must be >= 0 . Use poisson - gamma relation
287  // Pr ( x <= n) = Pr( y >= a) where x is poisson and y is gamma distributed ( a = n+1)
288  double a = (double) n + 1.0;
289  return ROOT::Math::gamma_cdf_c(mu, a, 1.0);
290  }
291 
292 
293  double binomial_cdf_c(unsigned int k, double p, unsigned int n)
294  {
295  // use relation with in beta distribution
296  // p must be 0 <=p <= 1
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);
301  }
302 
303 
304  double binomial_cdf(unsigned int k, double p, unsigned int n)
305  {
306  // use relation with in beta distribution
307  // p must be 0 <=p <= 1
308  if ( k >= n) return 1.0;
309 
310  double a = (double) k + 1.0;
311  double b = (double) n - k;
312  return ROOT::Math::beta_cdf_c(p, a, b);
313  }
314 
315 
316  double negative_binomial_cdf(unsigned int k, double p, double n)
317  {
318  // use relation with in beta distribution
319  // p must be 0 <=p <= 1
320  if (n < 0) return 0;
321  if (p < 0 || p > 1) return 0;
322  return ROOT::Math::beta_cdf(p, n, k+1.0);
323  }
324 
325 
326  double negative_binomial_cdf_c(unsigned int k, double p, double n)
327  {
328  // use relation with in beta distribution
329  // p must be 0 <=p <= 1
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);
333  }
334 
335 
336  double landau_cdf(double x, double xi, double x0)
337  {
338  // implementation of landau distribution (from DISLAN)
339  //The algorithm was taken from the Cernlib function dislan(G110)
340  //Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
341  //distribution", Computer Phys.Comm., 31(1984), 97-111
342 
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};
345 
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};
348 
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};
351 
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};
354 
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};
357 
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};
360 
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};
363 
364  double v = (x - x0)/xi;
365  double u;
366  double lan;
367 
368  if (v < -5.5)
369  {
370  u = std::exp(v+1);
371  lan = 0.3989422803*std::exp(-1./u)*std::sqrt(u)*(1+(a1[1]+(a1[2]+a1[3]*u)*u)*u);
372  }
373  else if (v < -1 )
374  {
375  u = std::exp(-v-1);
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);
378  }
379  else if (v < 1)
380  lan = (p2[0]+(p2[1]+(p2[2]+p2[3]*v)*v)*v)/(q2[0]+(q2[1]+(q2[2]+q2[3]*v)*v)*v);
381 
382  else if (v < 4)
383  lan = (p3[0]+(p3[1]+(p3[2]+p3[3]*v)*v)*v)/(q3[0]+(q3[1]+(q3[2]+q3[3]*v)*v)*v);
384 
385  else if (v < 12)
386  {
387  u = 1./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);
389  }
390  else if (v < 50)
391  {
392  u = 1./v;
393  lan = (p5[0]+(p5[1]+(p5[2]+p5[3]*u)*u)*u)/(q5[0]+(q5[1]+(q5[2]+q5[3]*u)*u)*u);
394  }
395  else if (v < 300)
396  {
397  u = 1./v;
398  lan = (p6[0]+(p6[1]+(p6[2]+p6[3]*u)*u)*u)/(q6[0]+(q6[1]+(q6[2]+q6[3]*u)*u)*u);
399  }
400  else
401  {
402  u = 1./(v-v*std::log(v)/(v+1));
403  lan = 1-(a2[1]+(a2[2]+a2[3]*u)*u)*u;
404  }
405  return lan;
406  }
407 
408 
409  double landau_xm1(double x, double xi, double x0)
410  {
411  // implementation of first momentum of Landau distribution
412  // translated from Cernlib (XM1LAN function) by Benno List
413 
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,
427  0.9026661865E+3};
428  static double q4[4] = { 1.0 , 0.7752562854E+2,-0.5637811998E+3,
429  -0.5513156752E+3};
430  static double p5[4] = { 0.9489335583E+0, 0.5561246706E+3, 0.3208274617E+5,
431  -0.4889926524E+5};
432  static double q5[4] = { 1.0 , 0.6028275940E+3, 0.3716962017E+5,
433  0.3686272898E+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,
437  -0.1641741416E+1};
438  static double a2[5] = { 0, -0.1958333333E+1, 0.5563368056E+1,
439  -0.2111352961E+2, 0.1006946266E+3};
440 
441  double v = (x-x0)/xi;
442  double xm1lan;
443  if (v < -4.5)
444  {
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);
448  }
449  else if (v < -2)
450  {
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);
453  }
454  else if (v < 2)
455  {
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);
458  }
459  else if (v < 10)
460  {
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);
463  }
464  else if (v < 40)
465  {
466  double u = 1/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);
469  }
470  else if (v < 200)
471  {
472  double u = 1/v;
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);
475  }
476  else
477  {
478  double u = v-v*std::log(v)/(v+1);
479  v = 1/(u-u*(u+ std::log(u)-v)/(u+1));
480  u = -std::log(v);
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);
483  }
484  return xm1lan*xi + x0;
485  }
486 
487 
488 
489  double landau_xm2(double x, double xi, double x0)
490  {
491  // implementation of second momentum of Landau distribution
492  // translated from Cernlib (XM2LAN function) by Benno List
493 
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,
503  0.3641361437E+2};
504  static double q3[4] = { 1.0 , 0.8614160194E+1, 0.3118929630E+2,
505  0.1514351300E+0};
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,
511  -0.7075963938E+5};
512  static double q5[4] = { 1.0 , 0.3605950254E+3, 0.1392784158E+5,
513  -0.1881680027E+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,
516  -0.1076714945E+2};
517  static double a1[4] = { 0. ,-0.4583333333E+0, 0.6675347222E+0,
518  -0.1641741416E+1};
519  static double a2[4] = {-0.1958333333E+1, 0.5563368056E+1,-0.2111352961E+2,
520  0.1006946266E+3};
521  static double a3[4] = {-1.0 , 0.4458333333E+1,-0.2116753472E+2,
522  0.1163674359E+3};
523 
524  double v = (x-x0)/xi;
525  double xm2lan;
526  if (v < -4.5)
527  {
528  double u = std::exp(v+1);
529  xm2lan = v*v-2*u*u*
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);
533  }
534  else if (v < -2)
535  {
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);
538  }
539  else if (v < 2)
540  {
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);
543  }
544  else if (v < 5)
545  {
546  double u = 1/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);
549  }
550  else if (v < 50)
551  {
552  double u = 1/v;
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);
555  }
556  else if (v < 200)
557  {
558  double u = 1/v;
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);
561  }
562  else
563  {
564  double u = v-v*std::log(v)/(v+1);
565  v = 1/(u-u*(u+log(u)-v)/(u+1));
566  u = -std::log(v);
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);
569  }
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;
573  }
574 
575 } // namespace Math
576 } // namespace ROOT
577 
578 
579 
static double B[]
double erf(double x)
Error function encountered in integrating the normal distribution.
double tdistribution_cdf(double x, double r, double x0=0)
Cumulative distribution function of Student&#39;s t-distribution (lower tail).
double binomial_cdf(unsigned int k, double p, unsigned int n)
Cumulative distribution function of the Binomial distribution Lower tail of the integral of the binom...
auto * m
Definition: textangle.C:8
VSD Structures.
Definition: StringConv.hxx:21
double crystalball_cdf(double x, double alpha, double n, double sigma, double x0=0)
Cumulative distribution for the Crystal Ball distribution function.
double inc_beta(double x, double a, double b)
Calculates the normalized (regularized) incomplete beta function.
double negative_binomial_cdf_c(unsigned int k, double p, double n)
Complement of the cumulative distribution function of the Negative Binomial distribution.
double poisson_cdf(unsigned int n, double mu)
Cumulative distribution function of the Poisson distribution Lower tail of the integral of the poisso...
double poisson_cdf_c(unsigned int n, double mu)
Complement of the cumulative distribution function of the Poisson distribution.
double binomial_cdf_c(unsigned int k, double p, unsigned int n)
Complement of the cumulative distribution function of the Binomial distribution.
static const double kSqrt2
double gamma_cdf_c(double x, double alpha, double theta, double x0=0)
Complement of the cumulative distribution function of the gamma distribution (upper tail)...
double uniform_cdf_c(double x, double a, double b, double x0=0)
Complement of the cumulative distribution function of the uniform (flat) distribution (upper tail)...
STL namespace.
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:79
static double A[]
double landau_xm1(double x, double xi=1, double x0=0)
First moment (mean) of the truncated Landau distribution.
double lognormal_cdf(double x, double m, double s, double x0=0)
Cumulative distribution function of the lognormal distribution (lower tail).
double erfc(double x)
Complementary error function.
double inc_gamma_c(double a, double x)
Calculates the normalized (regularized) upper incomplete gamma function (upper integral) ...
double negative_binomial_cdf(unsigned int k, double p, double n)
Cumulative distribution function of the Negative Binomial distribution Lower tail of the integral of ...
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
double fdistribution_cdf_c(double x, double n, double m, double x0=0)
Complement of the cumulative distribution function of the F-distribution (upper tail).
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
double pow(double, double)
double cauchy_cdf_c(double x, double b, double x0=0)
Complement of the cumulative distribution function (upper tail) of the Cauchy distribution which is a...
double fdistribution_cdf(double x, double n, double m, double x0=0)
Cumulative distribution function of the F-distribution (lower tail).
const Double_t sigma
static constexpr double s
double uniform_cdf(double x, double a, double b, double x0=0)
Cumulative distribution function of the uniform (flat) distribution (lower tail). ...
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:82
double beta_cdf(double x, double a, double b)
Cumulative distribution function of the beta distribution Upper tail of the integral of the beta_pdf...
double breitwigner_cdf(double x, double gamma, double x0=0)
Cumulative distribution function (lower tail) of the Breit_Wigner distribution and it is similar (jus...
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
double landau_xm2(double x, double xi=1, double x0=0)
Second moment of the truncated Landau distribution.
double gamma(double x)
static double C[]
ROOT::R::TRInterface & r
Definition: Object.C:4
#define M_PI
Definition: Rotated.cxx:105
double exponential_cdf_c(double x, double lambda, double x0=0)
Complement of the cumulative distribution function of the exponential distribution (upper tail)...
auto * a
Definition: textangle.C:12
double cauchy_cdf(double x, double b, double x0=0)
Cumulative distribution function (lower tail) of the Cauchy distribution which is also Lorentzian dis...
double expm1(double x)
exp(x) -1 with error cancellation when x is small
Definition: Math.h:110
double normal_cdf_c(double x, double sigma=1, double x0=0)
Complement of the cumulative distribution function of the normal (Gaussian) distribution (upper tail)...
double tdistribution_cdf_c(double x, double r, double x0=0)
Complement of the cumulative distribution function of Student&#39;s t-distribution (upper tail)...
double gaussian_cdf_c(double x, double sigma=1, double x0=0)
Alternative name for same function.
double atan(double)
Namespace for new Math classes and functions.
double beta_cdf_c(double x, double a, double b)
Complement of the cumulative distribution function of the beta distribution.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
double chisquared_cdf_c(double x, double r, double x0=0)
Complement of the cumulative distribution function of the distribution with degrees of freedom (upp...
double crystalball_integral(double x, double alpha, double n, double sigma, double x0=0)
Integral of the not-normalized Crystal Ball function.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double chisquared_cdf(double x, double r, double x0=0)
Cumulative distribution function of the distribution with degrees of freedom (lower tail)...
double breitwigner_cdf_c(double x, double gamma, double x0=0)
Complement of the cumulative distribution function (upper tail) of the Breit_Wigner distribution and ...
double exp(double)
double inc_gamma(double a, double x)
Calculates the normalized (regularized) lower incomplete gamma function (lower integral) ...
double crystalball_cdf_c(double x, double alpha, double n, double sigma, double x0=0)
Complement of the Cumulative distribution for the Crystal Ball distribution.
const Int_t n
Definition: legend1.C:16
double lognormal_cdf_c(double x, double m, double s, double x0=0)
Complement of the cumulative distribution function of the lognormal distribution (upper tail)...
double log(double)
double exponential_cdf(double x, double lambda, double x0=0)
Cumulative distribution function of the exponential distribution (lower tail).