ROOT   6.10/09 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...
static double p3(double t, double a, double b, double c, double d)
Namespace for new ROOT classes and functions.
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)...
TArc * a
Definition: textangle.C:12
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:47
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).
static double p2(double t, double a, double b, double c)
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
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:50
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[]
TRandom2 r(17)
#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)...
SVector< double, 2 > v
Definition: Dict.h:5
double cauchy_cdf(double x, double b, double x0=0)
Cumulative distribution function (lower tail) of the Cauchy distribution which is also Lorentzian dis...
TMarker * m
Definition: textangle.C:8
double expm1(double x)
exp(x) -1 with error cancellation when x is small
Definition: Math.h:112
static double p1(double t, double a, double b)
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).