Logo ROOT   6.13/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...
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
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: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 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).
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).
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 Int_t m
Definition: TRolke.cxx:637
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[]
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...
double expm1(double x)
exp(x) -1 with error cancellation when x is small
Definition: Math.h:112
#define M_PI
Definition: Math.h:38
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)...
* x
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:630
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.
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).