ROOT   6.14/05 Reference Guide
TMath.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: Rene Brun, Anna Kreshuk, Eddy Offermann, Fons Rademakers 29/07/95
3
4 /*************************************************************************
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11
12 ////////////////////////////////////////////////////////////////////////////////
13 // TMath
14 //
15 // Encapsulate math routines.
16
17 #include "TMath.h"
18 #include "TError.h"
19 #include <math.h>
20 #include <string.h>
21 #include <algorithm>
22 #include "Riostream.h"
23 #include "TString.h"
24
25 #include <Math/SpecFuncMathCore.h>
26 #include <Math/PdfFuncMathCore.h>
27 #include <Math/ProbFuncMathCore.h>
28
29 //const Double_t
30 // TMath::Pi = 3.14159265358979323846,
31 // TMath::E = 2.7182818284590452354;
32
33
34 // Without this macro the THtml doc for TMath can not be generated
35 #if !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
37 #endif
38
39 namespace TMath {
40
44  void VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt);
45
46 }
47
48 ////////////////////////////////////////////////////////////////////////////////
49
51 {
52  return (Long_t) (hypot((Double_t)x, (Double_t)y) + 0.5);
53 }
54
55 ////////////////////////////////////////////////////////////////////////////////
56
58 {
59  return hypot(x, y);
60 }
61
62 ////////////////////////////////////////////////////////////////////////////////
63
65 {
66 #if defined(WIN32)
67  if(x==0.0) return 0.0;
68  Double_t ax = Abs(x);
69  return log(x+ax*sqrt(1.+1./(ax*ax)));
70 #else
71  return asinh(x);
72 #endif
73 }
74
75 ////////////////////////////////////////////////////////////////////////////////
76
78 {
79 #if defined(WIN32)
80  if(x==0.0) return 0.0;
81  Double_t ax = Abs(x);
82  return log(x+ax*sqrt(1.-1./(ax*ax)));
83 #else
84  return acosh(x);
85 #endif
86 }
87
88 ////////////////////////////////////////////////////////////////////////////////
89
91 {
92 #if defined(WIN32)
93  return log((1+x)/(1-x))/2;
94 #else
95  return atanh(x);
96 #endif
97 }
98
99 ////////////////////////////////////////////////////////////////////////////////
100
102 {
103  return log(x)/log(2.0);
104 }
105
106 ////////////////////////////////////////////////////////////////////////////////
107 /// The DiLogarithm function
108 /// Code translated by R.Brun from CERNLIB DILOG function C332
109
111 {
112  const Double_t hf = 0.5;
113  const Double_t pi = TMath::Pi();
114  const Double_t pi2 = pi*pi;
115  const Double_t pi3 = pi2/3;
116  const Double_t pi6 = pi2/6;
117  const Double_t pi12 = pi2/12;
118  const Double_t c[20] = {0.42996693560813697, 0.40975987533077105,
119  -0.01858843665014592, 0.00145751084062268,-0.00014304184442340,
120  0.00001588415541880,-0.00000190784959387, 0.00000024195180854,
121  -0.00000003193341274, 0.00000000434545063,-0.00000000060578480,
122  0.00000000008612098,-0.00000000001244332, 0.00000000000182256,
123  -0.00000000000027007, 0.00000000000004042,-0.00000000000000610,
124  0.00000000000000093,-0.00000000000000014, 0.00000000000000002};
125
126  Double_t t,h,y,s,a,alfa,b1,b2,b0;
127  t=h=y=s=a=alfa=b1=b2=b0=0.;
128
129  if (x == 1) {
130  h = pi6;
131  } else if (x == -1) {
132  h = -pi12;
133  } else {
134  t = -x;
135  if (t <= -2) {
136  y = -1/(1+t);
137  s = 1;
138  b1= TMath::Log(-t);
139  b2= TMath::Log(1+1/t);
140  a = -pi3+hf*(b1*b1-b2*b2);
141  } else if (t < -1) {
142  y = -1-t;
143  s = -1;
144  a = TMath::Log(-t);
145  a = -pi6+a*(a+TMath::Log(1+1/t));
146  } else if (t <= -0.5) {
147  y = -(1+t)/t;
148  s = 1;
149  a = TMath::Log(-t);
150  a = -pi6+a*(-hf*a+TMath::Log(1+t));
151  } else if (t < 0) {
152  y = -t/(1+t);
153  s = -1;
154  b1= TMath::Log(1+t);
155  a = hf*b1*b1;
156  } else if (t <= 1) {
157  y = t;
158  s = 1;
159  a = 0;
160  } else {
161  y = 1/t;
162  s = -1;
163  b1= TMath::Log(t);
164  a = pi6+hf*b1*b1;
165  }
166  h = y+y-1;
167  alfa = h+h;
168  b1 = 0;
169  b2 = 0;
170  for (Int_t i=19;i>=0;i--){
171  b0 = c[i] + alfa*b1-b2;
172  b2 = b1;
173  b1 = b0;
174  }
175  h = -(s*(b0-h*b2)+a);
176  }
177  return h;
178 }
179
180 ////////////////////////////////////////////////////////////////////////////////
181 /// Computation of the error function erf(x).
182 /// Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x
183
185 {
187 }
188
189 ////////////////////////////////////////////////////////////////////////////////
190 /// Compute the complementary error function erfc(x).
191 /// Erfc(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between x and infinity
192 ///
193
195 {
197 }
198
199 ////////////////////////////////////////////////////////////////////////////////
200 /// returns the inverse error function
201 /// x must be <-1<x<1
202
204 {
205  Int_t kMaxit = 50;
206  Double_t kEps = 1e-14;
207  Double_t kConst = 0.8862269254527579; // sqrt(pi)/2.0
208
209  if(TMath::Abs(x) <= kEps) return kConst*x;
210
211  // Newton iterations
212  Double_t erfi, derfi, y0,y1,dy0,dy1;
213  if(TMath::Abs(x) < 1.0) {
214  erfi = kConst*TMath::Abs(x);
215  y0 = TMath::Erf(0.9*erfi);
216  derfi = 0.1*erfi;
217  for (Int_t iter=0; iter<kMaxit; iter++) {
218  y1 = 1. - TMath::Erfc(erfi);
219  dy1 = TMath::Abs(x) - y1;
220  if (TMath::Abs(dy1) < kEps) {if (x < 0) return -erfi; else return erfi;}
221  dy0 = y1 - y0;
222  derfi *= dy1/dy0;
223  y0 = y1;
224  erfi += derfi;
225  if(TMath::Abs(derfi/erfi) < kEps) {if (x < 0) return -erfi; else return erfi;}
226  }
227  }
228  return 0; //did not converge
229 }
230
231 ////////////////////////////////////////////////////////////////////////////////
232 /// returns the inverse of the complementary error function
233 /// x must be 0<x<2
234 /// implement using the quantile of the normal distribution
235 /// instead of ErfInverse for better numerical precision for large x
236
238 {
239
240  // erfc-1(x) = - 1/sqrt(2) * normal_quantile( 0.5 * x)
241  return - 0.70710678118654752440 * TMath::NormQuantile( 0.5 * x);
242 }
243
244 ////////////////////////////////////////////////////////////////////////////////
245 /// Compute factorial(n).
246
248 {
249  if (n <= 0) return 1.;
250  Double_t x = 1;
251  Int_t b = 0;
252  do {
253  b++;
254  x *= b;
255  } while (b != n);
256  return x;
257 }
258
259 ////////////////////////////////////////////////////////////////////////////////
260 /// Computation of the normal frequency function freq(x).
261 /// Freq(x) = (1/sqrt(2pi)) Integral(exp(-t^2/2))dt between -infinity and x.
262 ///
263 /// Translated from CERNLIB C300 by Rene Brun.
264
266 {
267  const Double_t c1 = 0.56418958354775629;
268  const Double_t w2 = 1.41421356237309505;
269
270  const Double_t p10 = 2.4266795523053175e+2, q10 = 2.1505887586986120e+2,
271  p11 = 2.1979261618294152e+1, q11 = 9.1164905404514901e+1,
272  p12 = 6.9963834886191355e+0, q12 = 1.5082797630407787e+1,
273  p13 =-3.5609843701815385e-2, q13 = 1;
274
275  const Double_t p20 = 3.00459261020161601e+2, q20 = 3.00459260956983293e+2,
276  p21 = 4.51918953711872942e+2, q21 = 7.90950925327898027e+2,
277  p22 = 3.39320816734343687e+2, q22 = 9.31354094850609621e+2,
278  p23 = 1.52989285046940404e+2, q23 = 6.38980264465631167e+2,
279  p24 = 4.31622272220567353e+1, q24 = 2.77585444743987643e+2,
280  p25 = 7.21175825088309366e+0, q25 = 7.70001529352294730e+1,
281  p26 = 5.64195517478973971e-1, q26 = 1.27827273196294235e+1,
282  p27 =-1.36864857382716707e-7, q27 = 1;
283
284  const Double_t p30 =-2.99610707703542174e-3, q30 = 1.06209230528467918e-2,
285  p31 =-4.94730910623250734e-2, q31 = 1.91308926107829841e-1,
286  p32 =-2.26956593539686930e-1, q32 = 1.05167510706793207e+0,
287  p33 =-2.78661308609647788e-1, q33 = 1.98733201817135256e+0,
288  p34 =-2.23192459734184686e-2, q34 = 1;
289
290  Double_t v = TMath::Abs(x)/w2;
291  Double_t vv = v*v;
292  Double_t ap, aq, h, hc, y;
293  if (v < 0.5) {
294  y=vv;
295  ap=p13;
296  aq=q13;
297  ap = p12 +y*ap;
298  ap = p11 +y*ap;
299  ap = p10 +y*ap;
300  aq = q12 +y*aq;
301  aq = q11 +y*aq;
302  aq = q10 +y*aq;
303  h = v*ap/aq;
304  hc = 1-h;
305  } else if (v < 4) {
306  ap = p27;
307  aq = q27;
308  ap = p26 +v*ap;
309  ap = p25 +v*ap;
310  ap = p24 +v*ap;
311  ap = p23 +v*ap;
312  ap = p22 +v*ap;
313  ap = p21 +v*ap;
314  ap = p20 +v*ap;
315  aq = q26 +v*aq;
316  aq = q25 +v*aq;
317  aq = q24 +v*aq;
318  aq = q23 +v*aq;
319  aq = q22 +v*aq;
320  aq = q21 +v*aq;
321  aq = q20 +v*aq;
322  hc = TMath::Exp(-vv)*ap/aq;
323  h = 1-hc;
324  } else {
325  y = 1/vv;
326  ap = p34;
327  aq = q34;
328  ap = p33 +y*ap;
329  ap = p32 +y*ap;
330  ap = p31 +y*ap;
331  ap = p30 +y*ap;
332  aq = q33 +y*aq;
333  aq = q32 +y*aq;
334  aq = q31 +y*aq;
335  aq = q30 +y*aq;
336  hc = TMath::Exp(-vv)*(c1+y*ap/aq)/v;
337  h = 1-hc;
338  }
339  if (x > 0) return 0.5 +0.5*h;
340  else return 0.5*hc;
341 }
342
343 ////////////////////////////////////////////////////////////////////////////////
344 /// Computation of gamma(z) for all z.
345 ///
346 /// C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
347
349 {
351 }
352
353 ////////////////////////////////////////////////////////////////////////////////
354 /// Computation of the normalized lower incomplete gamma function P(a,x) as defined in the
355 /// Handbook of Mathematical Functions by Abramowitz and Stegun, formula 6.5.1 on page 260 .
356 /// Its normalization is such that TMath::Gamma(a,+infinity) = 1 .
357 ///
358 /// \f[
359 /// P(a, x) = \frac{1}{\Gamma(a)} \int_{0}^{x} t^{a-1} e^{-t} dt
360 /// \f]
361 ///
362 /// \author NvE 14-nov-1998 UU-SAP Utrecht
363
365 {
367 }
368
369 ////////////////////////////////////////////////////////////////////////////////
370 /// Computation of the incomplete gamma function P(a,x)
371 /// via its continued fraction representation.
372 ///
373 /// \author NvE 14-nov-1998 UU-SAP Utrecht
374
376 {
377  Int_t itmax = 100; // Maximum number of iterations
378  Double_t eps = 3.e-14; // Relative accuracy
379  Double_t fpmin = 1.e-30; // Smallest Double_t value allowed here
380
381  if (a <= 0 || x <= 0) return 0;
382
383  Double_t gln = LnGamma(a);
384  Double_t b = x+1-a;
385  Double_t c = 1/fpmin;
386  Double_t d = 1/b;
387  Double_t h = d;
388  Double_t an,del;
389  for (Int_t i=1; i<=itmax; i++) {
390  an = Double_t(-i)*(Double_t(i)-a);
391  b += 2;
392  d = an*d+b;
393  if (Abs(d) < fpmin) d = fpmin;
394  c = b+an/c;
395  if (Abs(c) < fpmin) c = fpmin;
396  d = 1/d;
397  del = d*c;
398  h = h*del;
399  if (Abs(del-1) < eps) break;
400  //if (i==itmax) std::cout << "*GamCf(a,x)* a too large or itmax too small" << std::endl;
401  }
402  Double_t v = Exp(-x+a*Log(x)-gln)*h;
403  return (1-v);
404 }
405
406 ////////////////////////////////////////////////////////////////////////////////
407 /// Computation of the incomplete gamma function P(a,x)
408 /// via its series representation.
409 ///
410 /// \author NvE 14-nov-1998 UU-SAP Utrecht
411
413 {
414  Int_t itmax = 100; // Maximum number of iterations
415  Double_t eps = 3.e-14; // Relative accuracy
416
417  if (a <= 0 || x <= 0) return 0;
418
419  Double_t gln = LnGamma(a);
420  Double_t ap = a;
421  Double_t sum = 1/a;
422  Double_t del = sum;
423  for (Int_t n=1; n<=itmax; n++) {
424  ap += 1;
425  del = del*x/ap;
426  sum += del;
427  if (TMath::Abs(del) < Abs(sum*eps)) break;
428  //if (n==itmax) std::cout << "*GamSer(a,x)* a too large or itmax too small" << std::endl;
429  }
430  Double_t v = sum*Exp(-x+a*Log(x)-gln);
431  return v;
432 }
433
434 ////////////////////////////////////////////////////////////////////////////////
435 /// Calculate a Breit Wigner function with mean and gamma.
436
438 {
439  Double_t bw = gamma/((x-mean)*(x-mean) + gamma*gamma/4);
440  return bw/(2*Pi());
441 }
442
443 ////////////////////////////////////////////////////////////////////////////////
444 /// Calculate a gaussian function with mean and sigma.
445 /// If norm=kTRUE (default is kFALSE) the result is divided
446 /// by sqrt(2*Pi)*sigma.
447
449 {
450  if (sigma == 0) return 1.e30;
451  Double_t arg = (x-mean)/sigma;
452  // for |arg| > 39 result is zero in double precision
453  if (arg < -39.0 || arg > 39.0) return 0.0;
454  Double_t res = TMath::Exp(-0.5*arg*arg);
455  if (!norm) return res;
456  return res/(2.50662827463100024*sigma); //sqrt(2*Pi)=2.50662827463100024
457 }
458
459 ////////////////////////////////////////////////////////////////////////////////
460 /// The LANDAU function.
461 ///
462 /// mu is a location parameter and correspond approximately to the most probable value
463 /// and sigma is a scale parameter (not the sigma of the full distribution which is not defined)
464 /// Note that for mu=0 and sigma=1 (default values) the exact location of the maximum of the distribution
465 /// (most proper value) is at x = -0.22278
466 /// This function has been adapted from the CERNLIB routine G110 denlan.
467 /// If norm=kTRUE (default is kFALSE) the result is divided by sigma
468
470 {
471  if (sigma <= 0) return 0;
472  Double_t den = ::ROOT::Math::landau_pdf( (x-mu)/sigma );
473  if (!norm) return den;
474  return den/sigma;
475 }
476
477 ////////////////////////////////////////////////////////////////////////////////
478 /// Computation of ln[gamma(z)] for all z.
479 ///
480 /// C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
481 ///
482 /// The accuracy of the result is better than 2e-10.
483 ///
484 /// \author NvE 14-nov-1998 UU-SAP Utrecht
485
487 {
489 }
490
491 ////////////////////////////////////////////////////////////////////////////////
492 /// Normalize a vector v in place.
493 /// Returns the norm of the original vector.
494
496 {
497  Float_t d = Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
498  if (d != 0) {
499  v[0] /= d;
500  v[1] /= d;
501  v[2] /= d;
502  }
503  return d;
504 }
505
506 ////////////////////////////////////////////////////////////////////////////////
507 /// Normalize a vector v in place.
508 /// Returns the norm of the original vector.
509 /// This implementation (thanks Kevin Lynch <krlynch@bu.edu>) is protected
510 /// against possible overflows.
511
513 {
514  // Find the largest element, and divide that one out.
515
516  Double_t av0 = Abs(v[0]), av1 = Abs(v[1]), av2 = Abs(v[2]);
517
518  Double_t amax, foo, bar;
519  // 0 >= {1, 2}
520  if( av0 >= av1 && av0 >= av2 ) {
521  amax = av0;
522  foo = av1;
523  bar = av2;
524  }
525  // 1 >= {0, 2}
526  else if (av1 >= av0 && av1 >= av2) {
527  amax = av1;
528  foo = av0;
529  bar = av2;
530  }
531  // 2 >= {0, 1}
532  else {
533  amax = av2;
534  foo = av0;
535  bar = av1;
536  }
537
538  if (amax == 0.0)
539  return 0.;
540
541  Double_t foofrac = foo/amax, barfrac = bar/amax;
542  Double_t d = amax * Sqrt(1.+foofrac*foofrac+barfrac*barfrac);
543
544  v[0] /= d;
545  v[1] /= d;
546  v[2] /= d;
547  return d;
548 }
549
550 ////////////////////////////////////////////////////////////////////////////////
551 /// Compute the Poisson distribution function for (x,par)
552 /// The Poisson PDF is implemented by means of Euler's Gamma-function
553 /// (for the factorial), so for any x integer argument it is correct.
554 /// BUT for non-integer x values, it IS NOT equal to the Poisson distribution.
555 /// see TMath::PoissonI to get a non-smooth function.
556 /// Note that for large values of par, it is better to call
557 ///
558 /// TMath::Gaus(x,par,sqrt(par),kTRUE)
559 ///
560 /// Begin_Macro
561 /// {
562 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
563 /// TF1 *poisson = new TF1("poisson", "TMath::Poisson(x, 5)", 0, 15);
564 /// poisson->Draw("L");
565 /// }
566 /// End_Macro
567
569 {
570  if (x<0)
571  return 0;
572  else if (x == 0.0)
573  return 1./Exp(par);
574  else {
575  Double_t lnpoisson = x*log(par)-par-LnGamma(x+1.);
576  return Exp(lnpoisson);
577  }
578  // An alternative strategy is to transition to a Gaussian approximation for
579  // large values of par ...
580  // else {
581  // return Gaus(x,par,Sqrt(par),kTRUE);
582  // }
583 }
584
585 ////////////////////////////////////////////////////////////////////////////////
586 /// Compute the Poisson distribution function for (x,par)
587 /// This is a non-smooth function.
588 /// This function is equivalent to ROOT::Math::poisson_pdf
589 ///
590 /// Begin_Macro
591 /// {
592 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
593 /// TF1 *poissoni = new TF1("poissoni", "TMath::PoissonI(x, 5)", 0, 15);
594 /// poissoni->SetNpx(1000);
595 /// poissoni->Draw("L");
596 /// }
597 /// End_Macro
598
600 {
601  Int_t ix = Int_t(x);
602  return Poisson(ix,par);
603 }
604
605 ////////////////////////////////////////////////////////////////////////////////
606 /// Computation of the probability for a certain Chi-squared (chi2)
607 /// and number of degrees of freedom (ndf).
608 ///
609 /// Calculations are based on the incomplete gamma function P(a,x),
610 /// where a=ndf/2 and x=chi2/2.
611 ///
612 /// P(a,x) represents the probability that the observed Chi-squared
613 /// for a correct model should be less than the value chi2.
614 ///
615 /// The returned probability corresponds to 1-P(a,x),
616 /// which denotes the probability that an observed Chi-squared exceeds
617 /// the value chi2 by chance, even for a correct model.
618 ///
619 /// \author NvE 14-nov-1998 UU-SAP Utrecht
620
622 {
623  if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0
624
625  if (chi2 <= 0) {
626  if (chi2 < 0) return 0;
627  else return 1;
628  }
629
631 }
632
633 ////////////////////////////////////////////////////////////////////////////////
634 /// Calculates the Kolmogorov distribution function,
635 ///
636 /// \f[
637 /// P(z) = 2 \sum_{j=1}^{\infty} (-1)^{j-1} e^{-2 j^2 z^2}
638 /// \f]
639 ///
640 /// which gives the probability that Kolmogorov's test statistic will exceed
641 /// the value z assuming the null hypothesis. This gives a very powerful
642 /// test for comparing two one-dimensional distributions.
643 /// see, for example, Eadie et al, "statistical Methods in Experimental
644 /// Physics', pp 269-270).
645 ///
646 /// This function returns the confidence level for the null hypothesis, where:
647 /// - \f$z = dn \sqrt{n} \f$, and
648 /// - \f$dn \f$ is the maximum deviation between a hypothetical distribution
649 /// function and an experimental distribution with
650 /// - \f$n \f$ events
651 ///
652 /// NOTE: To compare two experimental distributions with m and n events,
653 /// use \f$z = \sqrt{m n/(m+n)) dn} \f$
654 ///
655 /// Accuracy: The function is far too accurate for any imaginable application.
656 /// Probabilities less than \f$10^{-15} \f$ are returned as zero.
657 /// However, remember that the formula is only valid for "large" n.
658 ///
659 /// Theta function inversion formula is used for z <= 1
660 ///
661 /// This function was translated by Rene Brun from PROBKL in CERNLIB.
662
664 {
665  Double_t fj[4] = {-2,-8,-18,-32}, r[4];
666  const Double_t w = 2.50662827;
667  // c1 - -pi**2/8, c2 = 9*c1, c3 = 25*c1
668  const Double_t c1 = -1.2337005501361697;
669  const Double_t c2 = -11.103304951225528;
670  const Double_t c3 = -30.842513753404244;
671
672  Double_t u = TMath::Abs(z);
673  Double_t p;
674  if (u < 0.2) {
675  p = 1;
676  } else if (u < 0.755) {
677  Double_t v = 1./(u*u);
678  p = 1 - w*(TMath::Exp(c1*v) + TMath::Exp(c2*v) + TMath::Exp(c3*v))/u;
679  } else if (u < 6.8116) {
680  r[1] = 0;
681  r[2] = 0;
682  r[3] = 0;
683  Double_t v = u*u;
684  Int_t maxj = TMath::Max(1,TMath::Nint(3./u));
685  for (Int_t j=0; j<maxj;j++) {
686  r[j] = TMath::Exp(fj[j]*v);
687  }
688  p = 2*(r[0] - r[1] +r[2] - r[3]);
689  } else {
690  p = 0;
691  }
692  return p;
693  }
694
695 ////////////////////////////////////////////////////////////////////////////////
696 /// Statistical test whether two one-dimensional sets of points are compatible
697 /// with coming from the same parent distribution, using the Kolmogorov test.
698 /// That is, it is used to compare two experimental distributions of unbinned data.
699 ///
700 /// ### Input:
701 /// a,b: One-dimensional arrays of length na, nb, respectively.
702 /// The elements of a and b must be given in ascending order.
703 /// option is a character string to specify options
704 /// "D" Put out a line of "Debug" printout
705 /// "M" Return the Maximum Kolmogorov distance instead of prob
706 ///
707 /// ### Output:
708 /// The returned value prob is a calculated confidence level which gives a
709 /// statistical test for compatibility of a and b.
710 /// Values of prob close to zero are taken as indicating a small probability
711 /// of compatibility. For two point sets drawn randomly from the same parent
712 /// distribution, the value of prob should be uniformly distributed between
713 /// zero and one.
714 /// in case of error the function return -1
715 /// If the 2 sets have a different number of points, the minimum of
716 /// the two sets is used.
717 ///
718 /// ### Method:
719 /// The Kolmogorov test is used. The test statistic is the maximum deviation
720 /// between the two integrated distribution functions, multiplied by the
721 /// normalizing factor (rdmax*sqrt(na*nb/(na+nb)).
722 ///
723 /// Code adapted by Rene Brun from CERNLIB routine TKOLMO (Fred James)
724 /// (W.T. Eadie, D. Drijard, F.E. James, M. Roos and B. Sadoulet,
725 /// Statistical Methods in Experimental Physics, (North-Holland,
726 /// Amsterdam 1971) 269-271)
727 ///
728 /// ### Method Improvement by Jason A Detwiler (JADetwiler@lbl.gov)
729 ///
730 /// The nuts-and-bolts of the TMath::KolmogorovTest() algorithm is a for-loop
731 /// over the two sorted arrays a and b representing empirical distribution
732 /// functions. The for-loop handles 3 cases: when the next points to be
733 /// evaluated satisfy a>b, a<b, or a=b:
734 ///
735 /// ~~~ {cpp}
736 /// for (Int_t i=0;i<na+nb;i++) {
737 /// if (a[ia-1] < b[ib-1]) {
738 /// rdiff -= sa;
739 /// ia++;
740 /// if (ia > na) {ok = kTRUE; break;}
741 /// } else if (a[ia-1] > b[ib-1]) {
742 /// rdiff += sb;
743 /// ib++;
744 /// if (ib > nb) {ok = kTRUE; break;}
745 /// } else {
746 /// rdiff += sb - sa;
747 /// ia++;
748 /// ib++;
749 /// if (ia > na) {ok = kTRUE; break;}
750 /// if (ib > nb) {ok = kTRUE; break;}
751 /// }
752 /// rdmax = TMath::Max(rdmax,TMath::Abs(rdiff));
753 /// }
754 /// ~~~
755 ///
756 /// For the last case, a=b, the algorithm advances each array by one index in an
757 /// attempt to move through the equality. However, this is incorrect when one or
758 /// the other of a or b (or both) have a repeated value, call it x. For the KS
759 /// statistic to be computed properly, rdiff needs to be calculated after all of
760 /// the a and b at x have been tallied (this is due to the definition of the
761 /// empirical distribution function; another way to convince yourself that the
762 /// old CERNLIB method is wrong is that it implies that the function defined as the
763 /// difference between a and b is multi-valued at x -- besides being ugly, this
764 /// would invalidate Kolmogorov's theorem).
765 ///
766 /// The solution is to just add while-loops into the equality-case handling to
767 /// perform the tally:
768 ///
769 /// ~~~ {cpp}
770 /// } else {
771 /// double x = a[ia-1];
772 /// while(a[ia-1] == x && ia <= na) {
773 /// rdiff -= sa;
774 /// ia++;
775 /// }
776 /// while(b[ib-1] == x && ib <= nb) {
777 /// rdiff += sb;
778 /// ib++;
779 /// }
780 /// if (ia > na) {ok = kTRUE; break;}
781 /// if (ib > nb) {ok = kTRUE; break;}
782 /// }
783 /// ~~~
784 ///
785 /// ### Note:
786 /// A good description of the Kolmogorov test can be seen at:
787 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda35g.htm
788
789 Double_t TMath::KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option)
790 {
791 // LM: Nov 2010: clean up and returns now a zero distance when vectors are the same
792
793  TString opt = option;
794  opt.ToUpper();
795
796  Double_t prob = -1;
797 // Require at least two points in each graph
798  if (!a || !b || na <= 2 || nb <= 2) {
799  Error("KolmogorovTest","Sets must have more than 2 points");
800  return prob;
801  }
802 // Constants needed
803  Double_t rna = na;
804  Double_t rnb = nb;
805  Double_t sa = 1./rna;
806  Double_t sb = 1./rnb;
807  Double_t rdiff = 0;
808  Double_t rdmax = 0;
809  Int_t ia = 0;
810  Int_t ib = 0;
811
812 // Main loop over point sets to find max distance
813 // rdiff is the running difference, and rdmax the max.
814  Bool_t ok = kFALSE;
815  for (Int_t i=0;i<na+nb;i++) {
816  if (a[ia] < b[ib]) {
817  rdiff -= sa;
818  ia++;
819  if (ia >= na) {ok = kTRUE; break;}
820  } else if (a[ia] > b[ib]) {
821  rdiff += sb;
822  ib++;
823  if (ib >= nb) {ok = kTRUE; break;}
824  } else {
825  // special cases for the ties
826  double x = a[ia];
827  while(ia < na && a[ia] == x) {
828  rdiff -= sa;
829  ia++;
830  }
831  while(ib < nb && b[ib] == x) {
832  rdiff += sb;
833  ib++;
834  }
835  if (ia >= na) {ok = kTRUE; break;}
836  if (ib >= nb) {ok = kTRUE; break;}
837  }
838  rdmax = TMath::Max(rdmax,TMath::Abs(rdiff));
839  }
840  // Should never terminate this loop with ok = kFALSE!
841  R__ASSERT(ok);
842
843  if (ok) {
844  rdmax = TMath::Max(rdmax,TMath::Abs(rdiff));
845  Double_t z = rdmax * TMath::Sqrt(rna*rnb/(rna+rnb));
846  prob = TMath::KolmogorovProb(z);
847  }
848  // debug printout
849  if (opt.Contains("D")) {
850  printf(" Kolmogorov Probability = %g, Max Dist = %g\n",prob,rdmax);
851  }
852  if(opt.Contains("M")) return rdmax;
853  else return prob;
854 }
855
856
857 ////////////////////////////////////////////////////////////////////////////////
858 /// Computation of Voigt function (normalised).
859 /// Voigt is a convolution of the two functions:
860 /// \f[
861 /// gauss(xx) = \frac{1}{(\sqrt{2\pi} sigma)} e^{\frac{xx^{2}}{(2 sigma{^2})}}
862 /// \f]
863 /// and
864 /// \f[
865 /// lorentz(xx) = \frac{ \frac{1}{\pi} \frac{lg}{2} }{ (xx^{2} + \frac{lg^{2}}{4}) }
866 /// \f]
867 /// .
868 ///
869 /// The Voigt function is known to be the real part of Faddeeva function also
870 /// called complex error function [2].
871 ///
872 /// The algoritm was developed by J. Humlicek [1].
873 /// This code is based on fortran code presented by R. J. Wells [2].
874 /// Translated and adapted by Miha D. Puc
875 ///
876 /// To calculate the Faddeeva function with relative error less than 10^(-r).
877 /// r can be set by the the user subject to the constraints 2 <= r <= 5.
878 ///
879 /// - [1] J. Humlicek, JQSRT, 21, 437 (1982).
880 /// - [2] [R.J. Wells "Rapid Approximation to the Voigt/Faddeeva Function and its Derivatives" JQSRT 62 (1999), pp 29-48.](http://www-atm.physics.ox.ac.uk/user/wells/voigt.html)
881
883 {
884  if ((sigma < 0 || lg < 0) || (sigma==0 && lg==0)) {
885  return 0; // Not meant to be for those who want to be thinner than 0
886  }
887
888  if (sigma == 0) {
889  return lg * 0.159154943 / (xx*xx + lg*lg /4); //pure Lorentz
890  }
891
892  if (lg == 0) { //pure gauss
893  return 0.39894228 / sigma * TMath::Exp(-xx*xx / (2*sigma*sigma));
894  }
895
896  Double_t x, y, k;
897  x = xx / sigma / 1.41421356;
898  y = lg / 2 / sigma / 1.41421356;
899
900  Double_t r0, r1;
901
902  if (r < 2) r = 2;
903  if (r > 5) r = 5;
904
905  r0=1.51 * exp(1.144 * (Double_t)r);
906  r1=1.60 * exp(0.554 * (Double_t)r);
907
908  // Constants
909
910  const Double_t rrtpi = 0.56418958; // 1/SQRT(pi)
911
912  Double_t y0, y0py0, y0q; // for CPF12 algorithm
913  y0 = 1.5;
914  y0py0 = y0 + y0;
915  y0q = y0 * y0;
916
917  Double_t c[6] = { 1.0117281, -0.75197147, 0.012557727, 0.010022008, -0.00024206814, 0.00000050084806};
918  Double_t s[6] = { 1.393237, 0.23115241, -0.15535147, 0.0062183662, 0.000091908299, -0.00000062752596};
919  Double_t t[6] = { 0.31424038, 0.94778839, 1.5976826, 2.2795071, 3.0206370, 3.8897249};
920
921  // Local variables
922
923  int j; // Loop variables
924  int rg1, rg2, rg3; // y polynomial flags
925  Double_t abx, xq, yq, yrrtpi; // --x--, x^2, y^2, y/SQRT(pi)
926  Double_t xlim0, xlim1, xlim2, xlim3, xlim4; // --x-- on region boundaries
927  Double_t a0=0, d0=0, d2=0, e0=0, e2=0, e4=0, h0=0, h2=0, h4=0, h6=0;// W4 temporary variables
928  Double_t p0=0, p2=0, p4=0, p6=0, p8=0, z0=0, z2=0, z4=0, z6=0, z8=0;
929  Double_t xp[6], xm[6], yp[6], ym[6]; // CPF12 temporary values
930  Double_t mq[6], pq[6], mf[6], pf[6];
931  Double_t d, yf, ypy0, ypy0q;
932
933  //***** Start of executable code *****************************************
934
935  rg1 = 1; // Set flags
936  rg2 = 1;
937  rg3 = 1;
938  yq = y * y; // y^2
939  yrrtpi = y * rrtpi; // y/SQRT(pi)
940
941  // Region boundaries when both k and L are required or when R<>4
942
943  xlim0 = r0 - y;
944  xlim1 = r1 - y;
945  xlim3 = 3.097 * y - 0.45;
946  xlim2 = 6.8 - y;
947  xlim4 = 18.1 * y + 1.65;
948  if ( y <= 1e-6 ) { // When y<10^-6 avoid W4 algorithm
949  xlim1 = xlim0;
950  xlim2 = xlim0;
951  }
952
953  abx = fabs(x); // |x|
954  xq = abx * abx; // x^2
955  if ( abx > xlim0 ) { // Region 0 algorithm
956  k = yrrtpi / (xq + yq);
957  } else if ( abx > xlim1 ) { // Humlicek W4 Region 1
958  if ( rg1 != 0 ) { // First point in Region 1
959  rg1 = 0;
960  a0 = yq + 0.5; // Region 1 y-dependents
961  d0 = a0*a0;
962  d2 = yq + yq - 1.0;
963  }
964  d = rrtpi / (d0 + xq*(d2 + xq));
965  k = d * y * (a0 + xq);
966  } else if ( abx > xlim2 ) { // Humlicek W4 Region 2
967  if ( rg2 != 0 ) { // First point in Region 2
968  rg2 = 0;
969  h0 = 0.5625 + yq * (4.5 + yq * (10.5 + yq * (6.0 + yq)));
970  // Region 2 y-dependents
971  h2 = -4.5 + yq * (9.0 + yq * ( 6.0 + yq * 4.0));
972  h4 = 10.5 - yq * (6.0 - yq * 6.0);
973  h6 = -6.0 + yq * 4.0;
974  e0 = 1.875 + yq * (8.25 + yq * (5.5 + yq));
975  e2 = 5.25 + yq * (1.0 + yq * 3.0);
976  e4 = 0.75 * h6;
977  }
978  d = rrtpi / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
979  k = d * y * (e0 + xq * (e2 + xq * (e4 + xq)));
980  } else if ( abx < xlim3 ) { // Humlicek W4 Region 3
981  if ( rg3 != 0 ) { // First point in Region 3
982  rg3 = 0;
983  z0 = 272.1014 + y * (1280.829 + y *
984  (2802.870 + y *
985  (3764.966 + y *
986  (3447.629 + y *
987  (2256.981 + y *
988  (1074.409 + y *
989  (369.1989 + y *
990  (88.26741 + y *
991  (13.39880 + y)
992  )))))))); // Region 3 y-dependents
993  z2 = 211.678 + y * (902.3066 + y *
994  (1758.336 + y *
995  (2037.310 + y *
996  (1549.675 + y *
997  (793.4273 + y *
998  (266.2987 + y *
999  (53.59518 + y * 5.0)
1000  ))))));
1001  z4 = 78.86585 + y * (308.1852 + y *
1002  (497.3014 + y *
1003  (479.2576 + y *
1004  (269.2916 + y *
1005  (80.39278 + y * 10.0)
1006  ))));
1007  z6 = 22.03523 + y * (55.02933 + y *
1008  (92.75679 + y *
1009  (53.59518 + y * 10.0)
1010  ));
1011  z8 = 1.496460 + y * (13.39880 + y * 5.0);
1012  p0 = 153.5168 + y * (549.3954 + y *
1013  (919.4955 + y *
1014  (946.8970 + y *
1015  (662.8097 + y *
1016  (328.2151 + y *
1017  (115.3772 + y *
1018  (27.93941 + y *
1019  (4.264678 + y * 0.3183291)
1020  )))))));
1021  p2 = -34.16955 + y * (-1.322256+ y *
1022  (124.5975 + y *
1023  (189.7730 + y *
1024  (139.4665 + y *
1025  (56.81652 + y *
1026  (12.79458 + y * 1.2733163)
1027  )))));
1028  p4 = 2.584042 + y * (10.46332 + y *
1029  (24.01655 + y *
1030  (29.81482 + y *
1031  (12.79568 + y * 1.9099744)
1032  )));
1033  p6 = -0.07272979 + y * (0.9377051 + y *
1034  (4.266322 + y * 1.273316));
1035  p8 = 0.0005480304 + y * 0.3183291;
1036  }
1037  d = 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 + xq * (z8 + xq)))));
1038  k = d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
1039  } else { // Humlicek CPF12 algorithm
1040  ypy0 = y + y0;
1041  ypy0q = ypy0 * ypy0;
1042  k = 0.0;
1043  for (j = 0; j <= 5; j++) {
1044  d = x - t[j];
1045  mq[j] = d * d;
1046  mf[j] = 1.0 / (mq[j] + ypy0q);
1047  xm[j] = mf[j] * d;
1048  ym[j] = mf[j] * ypy0;
1049  d = x + t[j];
1050  pq[j] = d * d;
1051  pf[j] = 1.0 / (pq[j] + ypy0q);
1052  xp[j] = pf[j] * d;
1053  yp[j] = pf[j] * ypy0;
1054  }
1055  if ( abx <= xlim4 ) { // Humlicek CPF12 Region I
1056  for (j = 0; j <= 5; j++) {
1057  k = k + c[j]*(ym[j]+yp[j]) - s[j]*(xm[j]-xp[j]) ;
1058  }
1059  } else { // Humlicek CPF12 Region II
1060  yf = y + y0py0;
1061  for ( j = 0; j <= 5; j++) {
1062  k = k + (c[j] *
1063  (mq[j] * mf[j] - y0 * ym[j])
1064  + s[j] * yf * xm[j]) / (mq[j]+y0q)
1065  + (c[j] * (pq[j] * pf[j] - y0 * yp[j])
1066  - s[j] * yf * xp[j]) / (pq[j]+y0q);
1067  }
1068  k = y * k + exp( -xq );
1069  }
1070  }
1071  return k / 2.506628 / sigma; // Normalize by dividing by sqrt(2*pi)*sigma.
1072 }
1073
1074 ////////////////////////////////////////////////////////////////////////////////
1075 /// Calculates roots of polynomial of 3rd order a*x^3 + b*x^2 + c*x + d, where
1076 /// - a == coef[3],
1077 /// - b == coef[2],
1078 /// - c == coef[1],
1079 /// - d == coef[0]
1080 ///
1081 ///coef[3] must be different from 0
1082 ///
1083 /// If the boolean returned by the method is false:
1084 /// ==> there are 3 real roots a,b,c
1085
1086 /// If the boolean returned by the method is true:
1087 /// ==> there is one real root a and 2 complex conjugates roots (b+i*c,b-i*c)
1088 ///
1089 /// \author Francois-Xavier Gentit
1090
1092 {
1093  Bool_t complex = kFALSE;
1094  Double_t r,s,t,p,q,d,ps3,ps33,qs2,u,v,tmp,lnu,lnv,su,sv,y1,y2,y3;
1095  a = 0;
1096  b = 0;
1097  c = 0;
1098  if (coef[3] == 0) return complex;
1099  r = coef[2]/coef[3];
1100  s = coef[1]/coef[3];
1101  t = coef[0]/coef[3];
1102  p = s - (r*r)/3;
1103  ps3 = p/3;
1104  q = (2*r*r*r)/27.0 - (r*s)/3 + t;
1105  qs2 = q/2;
1106  ps33 = ps3*ps3*ps3;
1107  d = ps33 + qs2*qs2;
1108  if (d>=0) {
1109  complex = kTRUE;
1110  d = TMath::Sqrt(d);
1111  u = -qs2 + d;
1112  v = -qs2 - d;
1113  tmp = 1./3.;
1114  lnu = TMath::Log(TMath::Abs(u));
1115  lnv = TMath::Log(TMath::Abs(v));
1116  su = TMath::Sign(1.,u);
1117  sv = TMath::Sign(1.,v);
1118  u = su*TMath::Exp(tmp*lnu);
1119  v = sv*TMath::Exp(tmp*lnv);
1120  y1 = u + v;
1121  y2 = -y1/2;
1122  y3 = ((u-v)*TMath::Sqrt(3.))/2;
1123  tmp = r/3;
1124  a = y1 - tmp;
1125  b = y2 - tmp;
1126  c = y3;
1127  } else {
1128  Double_t phi,cphi,phis3,c1,c2,c3,pis3;
1129  ps3 = -ps3;
1130  ps33 = -ps33;
1131  cphi = -qs2/TMath::Sqrt(ps33);
1132  phi = TMath::ACos(cphi);
1133  phis3 = phi/3;
1134  pis3 = TMath::Pi()/3;
1135  c1 = TMath::Cos(phis3);
1136  c2 = TMath::Cos(pis3 + phis3);
1137  c3 = TMath::Cos(pis3 - phis3);
1138  tmp = TMath::Sqrt(ps3);
1139  y1 = 2*tmp*c1;
1140  y2 = -2*tmp*c2;
1141  y3 = -2*tmp*c3;
1142  tmp = r/3;
1143  a = y1 - tmp;
1144  b = y2 - tmp;
1145  c = y3 - tmp;
1146  }
1147  return complex;
1148 }
1149
1150 ////////////////////////////////////////////////////////////////////////////////
1151 ///Computes sample quantiles, corresponding to the given probabilities
1152 ///
1153 /// \param[in] x the data sample
1154 /// \param[in] n its size
1155 /// \param[out] quantiles computed quantiles are returned in there
1156 /// \param[in] prob probabilities where to compute quantiles
1157 /// \param[in] nprob size of prob array
1158 /// \param[in] isSorted is the input array x sorted ?
1159 /// \param[in] type method to compute (from 1 to 9).
1160 ///
1161 /// #### NOTE:
1162 /// When the input is not sorted, an array of integers of size n needs
1163 /// to be allocated. It can be passed by the user in parameter index,
1164 /// or, if not passed, it will be allocated inside the function
1165 ///
1166 /// ### Following types are provided:
1167 /// - Discontinuous:
1168 /// - type=1 - inverse of the empirical distribution function
1169 /// - type=2 - like type 1, but with averaging at discontinuities
1170 /// - type=3 - SAS definition: nearest even order statistic
1171 /// - Piecewise linear continuous:
1172 /// - In this case, sample quantiles can be obtained by linear interpolation
1173 /// between the k-th order statistic and p(k).
1174 /// -type=4 - linear interpolation of empirical cdf, p(k)=k/n;
1175 /// - type=5 - a very popular definition, p(k) = (k-0.5)/n;
1176 /// - type=6 - used by Minitab and SPSS, p(k) = k/(n+1);
1177 /// - type=7 - used by S-Plus and R, p(k) = (k-1)/(n-1);
1178 /// - type=8 - resulting sample quantiles are approximately median unbiased
1179 /// regardless of the distribution of x. p(k) = (k-1/3)/(n+1/3);
1180 /// - type=9 - resulting sample quantiles are approximately unbiased, when
1181 /// the sample comes from Normal distribution. p(k)=(k-3/8)/(n+1/4);
1182 ///
1183 /// default type = 7
1184 ///
1185 /// ### References:
1186 /// 1. Hyndman, R.J and Fan, Y, (1996) "Sample quantiles in statistical packages"
1187 /// American Statistician, 50, 361-365
1188 /// 2. R Project documentation for the function quantile of package {stats}
1189
1190 void TMath::Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted, Int_t *index, Int_t type)
1191 {
1192
1193  if (type<1 || type>9){
1194  printf("illegal value of type\n");
1195  return;
1196  }
1197  Int_t *ind = 0;
1198  Bool_t isAllocated = kFALSE;
1199  if (!isSorted){
1200  if (index) ind = index;
1201  else {
1202  ind = new Int_t[n];
1203  isAllocated = kTRUE;
1204  }
1205  }
1206
1207  // re-implemented by L.M (9/11/2011) to fix bug https://savannah.cern.ch/bugs/?87251
1208  // following now exactly R implementation (available in library/stats/R/quantile.R )
1209  // which follows precisely Hyndman-Fan paper
1210  // (older implementation had a bug for type =3)
1211
1212  for (Int_t i=0; i<nprob; i++){
1213
1214  Double_t nppm = 0;
1215  Double_t gamma = 0;
1216  Int_t j = 0;
1217
1218  //Discontinuous functions
1219  // type = 1,2, or 3
1220  if (type < 4 ){
1221  if (type == 3)
1222  nppm = n*prob[i]-0.5; // use m = -0.5
1223  else
1224  nppm = n*prob[i]; // use m = 0
1225
1226  // be careful with machine precision
1227  double eps = 4 * TMath::Limits<Double_t>::Epsilon();
1228  j = TMath::FloorNint(nppm + eps);
1229
1230  // LM : fix for numerical problems if nppm is actually equal to j, but results different for numerical error
1231  // g in the paper is nppm -j
1232  if (type == 1)
1233  gamma = ( (nppm -j) > j*TMath::Limits<Double_t>::Epsilon() ) ? 1 : 0;
1234  else if (type == 2)
1235  gamma = ( (nppm -j) > j*TMath::Limits<Double_t>::Epsilon() ) ? 1 : 0.5;
1236  else if (type == 3)
1237  // gamma = 0 for g=0 and j even
1238  gamma = ( TMath::Abs(nppm -j) <= j*TMath::Limits<Double_t>::Epsilon() && TMath::Even(j) ) ? 0 : 1;
1239
1240  }
1241  else {
1242  // for continuous quantiles type 4 to 9)
1243  // define alpha and beta
1244  double a = 0;
1245  double b = 0;
1246  if (type == 4) { a = 0; b = 1; }
1247  else if (type == 5) { a = 0.5; b = 0.5; }
1248  else if (type == 6) { a = 0.; b = 0.; }
1249  else if (type == 7) { a = 1.; b = 1.; }
1250  else if (type == 8) { a = 1./3.; b = a; }
1251  else if (type == 9) { a = 3./8.; b = a; }
1252
1253  // be careful with machine precision
1254  double eps = 4 * TMath::Limits<Double_t>::Epsilon();
1255  nppm = a + prob[i] * ( n + 1 -a -b); // n * p + m
1256  j = TMath::FloorNint(nppm + eps);
1257  gamma = nppm - j;
1258  if (gamma < eps) gamma = 0;
1259  }
1260
1261  // since index j starts from 1 first is j-1 and second is j
1262  // add protection to keep indices in range [0,n-1]
1263  int first = (j > 0 && j <=n) ? j-1 : ( j <=0 ) ? 0 : n-1;
1264  int second = (j > 0 && j < n) ? j : ( j <=0 ) ? 0 : n-1;
1265
1266  Double_t xj, xjj;
1267  if (isSorted){
1268  xj = x[first];
1269  xjj = x[second];
1270  } else {
1271  xj = TMath::KOrdStat(n, x, first, ind);
1272  xjj = TMath::KOrdStat(n, x, second, ind);
1273  }
1274  quantiles[i] = (1-gamma)*xj + gamma*xjj;
1275  // printf("x[j] = %12f x[j+1] = %12f \n",xj, xjj);
1276
1277  }
1278
1279
1280
1281  if (isAllocated)
1282  delete [] ind;
1283 }
1284
1285 ////////////////////////////////////////////////////////////////////////////////
1286 /// Bubble sort variant to obtain the order of an array's elements into
1287 /// an index in order to do more useful things than the standard built
1288 /// in functions.
1289 /// \param[in] *arr1 is unchanged;
1290 /// \param[in] *arr2 is the array of indicies corresponding to the descending value
1291 /// of arr1 with arr2[0] corresponding to the largest arr1 value and
1292 /// arr2[Narr] the smallest.
1293 ///
1294 /// \author Adrian Bevan (bevan@slac.stanford.edu)
1295
1296 void TMath::BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
1297 {
1298  if (Narr <= 0) return;
1299  double *localArr1 = new double[Narr];
1300  int *localArr2 = new int[Narr];
1301  int iEl;
1302  int iEl2;
1303
1304  for(iEl = 0; iEl < Narr; iEl++) {
1305  localArr1[iEl] = arr1[iEl];
1306  localArr2[iEl] = iEl;
1307  }
1308
1309  for (iEl = 0; iEl < Narr; iEl++) {
1310  for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1311  if (localArr1[iEl2-1] < localArr1[iEl2]) {
1312  double tmp = localArr1[iEl2-1];
1313  localArr1[iEl2-1] = localArr1[iEl2];
1314  localArr1[iEl2] = tmp;
1315
1316  int tmp2 = localArr2[iEl2-1];
1317  localArr2[iEl2-1] = localArr2[iEl2];
1318  localArr2[iEl2] = tmp2;
1319  }
1320  }
1321  }
1322
1323  for (iEl = 0; iEl < Narr; iEl++) {
1324  arr2[iEl] = localArr2[iEl];
1325  }
1326  delete [] localArr2;
1327  delete [] localArr1;
1328 }
1329
1330 ////////////////////////////////////////////////////////////////////////////////
1331 /// Opposite ordering of the array arr2[] to that of BubbleHigh.
1332 ///
1333 /// \author Adrian Bevan (bevan@slac.stanford.edu)
1334
1335 void TMath::BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
1336 {
1337  if (Narr <= 0) return;
1338  double *localArr1 = new double[Narr];
1339  int *localArr2 = new int[Narr];
1340  int iEl;
1341  int iEl2;
1342
1343  for (iEl = 0; iEl < Narr; iEl++) {
1344  localArr1[iEl] = arr1[iEl];
1345  localArr2[iEl] = iEl;
1346  }
1347
1348  for (iEl = 0; iEl < Narr; iEl++) {
1349  for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1350  if (localArr1[iEl2-1] > localArr1[iEl2]) {
1351  double tmp = localArr1[iEl2-1];
1352  localArr1[iEl2-1] = localArr1[iEl2];
1353  localArr1[iEl2] = tmp;
1354
1355  int tmp2 = localArr2[iEl2-1];
1356  localArr2[iEl2-1] = localArr2[iEl2];
1357  localArr2[iEl2] = tmp2;
1358  }
1359  }
1360  }
1361
1362  for (iEl = 0; iEl < Narr; iEl++) {
1363  arr2[iEl] = localArr2[iEl];
1364  }
1365  delete [] localArr2;
1366  delete [] localArr1;
1367 }
1368
1369
1370 ////////////////////////////////////////////////////////////////////////////////
1371 /// Calculates hash index from any char string.
1372 /// Based on pre-calculated table of 256 specially selected numbers.
1373 /// These numbers are selected in such a way, that for string
1374 /// length == 4 (integer number) the hash is unambiguous, i.e.
1375 /// from hash value we can recalculate input (no degeneration).
1376 ///
1377 /// The quality of hash method is good enough, that
1378 /// "random" numbers made as R = Hash(1), Hash(2), ...Hash(N)
1379 /// tested by <R>, <R*R>, <Ri*Ri+1> gives the same result
1380 /// as for libc rand().
1381 ///
1382 /// - For string: i = TMath::Hash(string,nstring);
1383 /// - For int: i = TMath::Hash(&intword,sizeof(int));
1384 /// - For pointer: i = TMath::Hash(&pointer,sizeof(void*));
1385 ///
1386 /// V.Perev
1387 /// This function is kept for back compatibility. The code previously in this function
1388 /// has been moved to the static function TString::Hash
1389
1390 ULong_t TMath::Hash(const void *txt, Int_t ntxt)
1391 {
1392  return TString::Hash(txt,ntxt);
1393 }
1394
1395
1396 ////////////////////////////////////////////////////////////////////////////////
1397
1398 ULong_t TMath::Hash(const char *txt)
1399 {
1400  return Hash(txt, Int_t(strlen(txt)));
1401 }
1402
1403 ////////////////////////////////////////////////////////////////////////////////
1404 /// Compute the modified Bessel function I_0(x) for any real x.
1405 ///
1406 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1407
1409 {
1410  // Parameters of the polynomial approximation
1411  const Double_t p1=1.0, p2=3.5156229, p3=3.0899424,
1412  p4=1.2067492, p5=0.2659732, p6=3.60768e-2, p7=4.5813e-3;
1413
1414  const Double_t q1= 0.39894228, q2= 1.328592e-2, q3= 2.25319e-3,
1415  q4=-1.57565e-3, q5= 9.16281e-3, q6=-2.057706e-2,
1416  q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3;
1417
1418  const Double_t k1 = 3.75;
1419  Double_t ax = TMath::Abs(x);
1420
1421  Double_t y=0, result=0;
1422
1423  if (ax < k1) {
1424  Double_t xx = x/k1;
1425  y = xx*xx;
1426  result = p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))));
1427  } else {
1428  y = k1/ax;
1429  result = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
1430  }
1431  return result;
1432 }
1433
1434 ////////////////////////////////////////////////////////////////////////////////
1435 /// Compute the modified Bessel function K_0(x) for positive real x.
1436 ///
1437 /// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1438 /// Applied Mathematics Series vol. 55 (1964), Washington.
1439 ///
1440 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1441
1443 {
1444  // Parameters of the polynomial approximation
1445  const Double_t p1=-0.57721566, p2=0.42278420, p3=0.23069756,
1446  p4= 3.488590e-2, p5=2.62698e-3, p6=1.0750e-4, p7=7.4e-6;
1447
1448  const Double_t q1= 1.25331414, q2=-7.832358e-2, q3= 2.189568e-2,
1449  q4=-1.062446e-2, q5= 5.87872e-3, q6=-2.51540e-3, q7=5.3208e-4;
1450
1451  if (x <= 0) {
1452  Error("TMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
1453  return 0;
1454  }
1455
1456  Double_t y=0, result=0;
1457
1458  if (x <= 2) {
1459  y = x*x/4;
1460  result = (-log(x/2.)*TMath::BesselI0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1461  } else {
1462  y = 2/x;
1463  result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
1464  }
1465  return result;
1466 }
1467
1468 ////////////////////////////////////////////////////////////////////////////////
1469 /// Compute the modified Bessel function I_1(x) for any real x.
1470 ///
1471 /// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1472 /// Applied Mathematics Series vol. 55 (1964), Washington.
1473 ///
1474 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1475
1477 {
1478  // Parameters of the polynomial approximation
1479  const Double_t p1=0.5, p2=0.87890594, p3=0.51498869,
1480  p4=0.15084934, p5=2.658733e-2, p6=3.01532e-3, p7=3.2411e-4;
1481
1482  const Double_t q1= 0.39894228, q2=-3.988024e-2, q3=-3.62018e-3,
1483  q4= 1.63801e-3, q5=-1.031555e-2, q6= 2.282967e-2,
1484  q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3;
1485
1486  const Double_t k1 = 3.75;
1487  Double_t ax = TMath::Abs(x);
1488
1489  Double_t y=0, result=0;
1490
1491  if (ax < k1) {
1492  Double_t xx = x/k1;
1493  y = xx*xx;
1494  result = x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1495  } else {
1496  y = k1/ax;
1497  result = (exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
1498  if (x < 0) result = -result;
1499  }
1500  return result;
1501 }
1502
1503 ////////////////////////////////////////////////////////////////////////////////
1504 /// Compute the modified Bessel function K_1(x) for positive real x.
1505 ///
1506 /// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1507 /// Applied Mathematics Series vol. 55 (1964), Washington.
1508 ///
1509 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1510
1512 {
1513  // Parameters of the polynomial approximation
1514  const Double_t p1= 1., p2= 0.15443144, p3=-0.67278579,
1515  p4=-0.18156897, p5=-1.919402e-2, p6=-1.10404e-3, p7=-4.686e-5;
1516
1517  const Double_t q1= 1.25331414, q2= 0.23498619, q3=-3.655620e-2,
1518  q4= 1.504268e-2, q5=-7.80353e-3, q6= 3.25614e-3, q7=-6.8245e-4;
1519
1520  if (x <= 0) {
1521  Error("TMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
1522  return 0;
1523  }
1524
1525  Double_t y=0,result=0;
1526
1527  if (x <= 2) {
1528  y = x*x/4;
1529  result = (log(x/2.)*TMath::BesselI1(x))+(1./x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1530  } else {
1531  y = 2/x;
1532  result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
1533  }
1534  return result;
1535 }
1536
1537 ////////////////////////////////////////////////////////////////////////////////
1538 /// Compute the Integer Order Modified Bessel function K_n(x)
1539 /// for n=0,1,2,... and positive real x.
1540 ///
1541 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1542
1544 {
1545  if (x <= 0 || n < 0) {
1546  Error("TMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
1547  return 0;
1548  }
1549
1550  if (n==0) return TMath::BesselK0(x);
1551  if (n==1) return TMath::BesselK1(x);
1552
1553  // Perform upward recurrence for all x
1554  Double_t tox = 2/x;
1555  Double_t bkm = TMath::BesselK0(x);
1556  Double_t bk = TMath::BesselK1(x);
1557  Double_t bkp = 0;
1558  for (Int_t j=1; j<n; j++) {
1559  bkp = bkm+Double_t(j)*tox*bk;
1560  bkm = bk;
1561  bk = bkp;
1562  }
1563  return bk;
1564 }
1565
1566 ////////////////////////////////////////////////////////////////////////////////
1567 /// Compute the Integer Order Modified Bessel function I_n(x)
1568 /// for n=0,1,2,... and any real x.
1569 ///
1570 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1571
1573 {
1574  Int_t iacc = 40; // Increase to enhance accuracy
1575  const Double_t kBigPositive = 1.e10;
1576  const Double_t kBigNegative = 1.e-10;
1577
1578  if (n < 0) {
1579  Error("TMath::BesselI", "*I* Invalid argument (n,x) = (%d, %g)\n",n,x);
1580  return 0;
1581  }
1582
1583  if (n==0) return TMath::BesselI0(x);
1584  if (n==1) return TMath::BesselI1(x);
1585
1586  if (x == 0) return 0;
1587  if (TMath::Abs(x) > kBigPositive) return 0;
1588
1589  Double_t tox = 2/TMath::Abs(x);
1590  Double_t bip = 0, bim = 0;
1591  Double_t bi = 1;
1592  Double_t result = 0;
1593  Int_t m = 2*((n+Int_t(sqrt(Float_t(iacc*n)))));
1594  for (Int_t j=m; j>=1; j--) {
1595  bim = bip+Double_t(j)*tox*bi;
1596  bip = bi;
1597  bi = bim;
1598  // Renormalise to prevent overflows
1599  if (TMath::Abs(bi) > kBigPositive) {
1600  result *= kBigNegative;
1601  bi *= kBigNegative;
1602  bip *= kBigNegative;
1603  }
1604  if (j==n) result=bip;
1605  }
1606
1607  result *= TMath::BesselI0(x)/bi; // Normalise with BesselI0(x)
1608  if ((x < 0) && (n%2 == 1)) result = -result;
1609
1610  return result;
1611 }
1612
1613 ////////////////////////////////////////////////////////////////////////////////
1614 /// Returns the Bessel function J0(x) for any real x.
1615
1617 {
1618  Double_t ax,z;
1619  Double_t xx,y,result,result1,result2;
1620  const Double_t p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
1621  const Double_t p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
1622  const Double_t p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
1623  const Double_t p10 = 59272.64853, p11 = 267.8532712;
1624
1625  const Double_t q1 = 0.785398164;
1626  const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1627  const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1628  const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1629  const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1630  const Double_t q10 = 0.934935152e-7, q11 = 0.636619772;
1631
1632  if ((ax=fabs(x)) < 8) {
1633  y=x*x;
1634  result1 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6))));
1635  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1636  result = result1/result2;
1637  } else {
1638  z = 8/ax;
1639  y = z*z;
1640  xx = ax-q1;
1641  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1642  result2 = q6 + y*(q7 + y*(q8 + y*(q9 - y*q10)));
1643  result = sqrt(q11/ax)*(cos(xx)*result1-z*sin(xx)*result2);
1644  }
1645  return result;
1646 }
1647
1648 ////////////////////////////////////////////////////////////////////////////////
1649 /// Returns the Bessel function J1(x) for any real x.
1650
1652 {
1653  Double_t ax,z;
1654  Double_t xx,y,result,result1,result2;
1655  const Double_t p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
1656  const Double_t p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
1657  const Double_t p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
1658  const Double_t p10 = 99447.43394, p11 = 376.9991397;
1659
1660  const Double_t q1 = 2.356194491;
1661  const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1662  const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1663  const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1664  const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1665  const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1666
1667  if ((ax=fabs(x)) < 8) {
1668  y=x*x;
1669  result1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6)))));
1670  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1671  result = result1/result2;
1672  } else {
1673  z = 8/ax;
1674  y = z*z;
1675  xx = ax-q1;
1676  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1677  result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1678  result = sqrt(q11/ax)*(cos(xx)*result1-z*sin(xx)*result2);
1679  if (x < 0) result = -result;
1680  }
1681  return result;
1682 }
1683
1684 ////////////////////////////////////////////////////////////////////////////////
1685 /// Returns the Bessel function Y0(x) for positive x.
1686
1688 {
1689  Double_t z,xx,y,result,result1,result2;
1690  const Double_t p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
1691  const Double_t p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
1692  const Double_t p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
1693  const Double_t p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
1694
1695  const Double_t q1 = 0.785398164;
1696  const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1697  const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1698  const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1699  const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1700  const Double_t q10 = -0.934945152e-7, q11 = 0.636619772;
1701
1702  if (x < 8) {
1703  y = x*x;
1704  result1 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6))));
1705  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1706  result = (result1/result2) + p12*TMath::BesselJ0(x)*log(x);
1707  } else {
1708  z = 8/x;
1709  y = z*z;
1710  xx = x-q1;
1711  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1712  result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1713  result = sqrt(q11/x)*(sin(xx)*result1+z*cos(xx)*result2);
1714  }
1715  return result;
1716 }
1717
1718 ////////////////////////////////////////////////////////////////////////////////
1719 /// Returns the Bessel function Y1(x) for positive x.
1720
1722 {
1723  Double_t z,xx,y,result,result1,result2;
1724  const Double_t p1 = -0.4900604943e13, p2 = 0.1275274390e13;
1725  const Double_t p3 = -0.5153438139e11, p4 = 0.7349264551e9;
1726  const Double_t p5 = -0.4237922726e7, p6 = 0.8511937935e4;
1727  const Double_t p7 = 0.2499580570e14, p8 = 0.4244419664e12;
1728  const Double_t p9 = 0.3733650367e10, p10 = 0.2245904002e8;
1729  const Double_t p11 = 0.1020426050e6, p12 = 0.3549632885e3;
1730  const Double_t p13 = 0.636619772;
1731  const Double_t q1 = 2.356194491;
1732  const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1733  const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1734  const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1735  const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1736  const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1737
1738  if (x < 8) {
1739  y=x*x;
1740  result1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6)))));
1741  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y*(p12+y)))));
1742  result = (result1/result2) + p13*(TMath::BesselJ1(x)*log(x)-1/x);
1743  } else {
1744  z = 8/x;
1745  y = z*z;
1746  xx = x-q1;
1747  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1748  result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1749  result = sqrt(q11/x)*(sin(xx)*result1+z*cos(xx)*result2);
1750  }
1751  return result;
1752 }
1753
1754 ////////////////////////////////////////////////////////////////////////////////
1755 /// Struve Functions of Order 0
1756 ///
1757 /// Converted from CERNLIB M342 by Rene Brun.
1758
1760 {
1761  const Int_t n1 = 15;
1762  const Int_t n2 = 25;
1763  const Double_t c1[16] = { 1.00215845609911981, -1.63969292681309147,
1764  1.50236939618292819, -.72485115302121872,
1765  .18955327371093136, -.03067052022988,
1766  .00337561447375194, -2.6965014312602e-4,
1767  1.637461692612e-5, -7.8244408508e-7,
1768  3.021593188e-8, -9.6326645e-10,
1769  2.579337e-11, -5.8854e-13,
1770  1.158e-14, -2e-16 };
1771  const Double_t c2[26] = { .99283727576423943, -.00696891281138625,
1772  1.8205103787037e-4, -1.063258252844e-5,
1773  9.8198294287e-7, -1.2250645445e-7,
1774  1.894083312e-8, -3.44358226e-9,
1775  7.1119102e-10, -1.6288744e-10,
1776  4.065681e-11, -1.091505e-11,
1777  3.12005e-12, -9.4202e-13,
1778  2.9848e-13, -9.872e-14,
1779  3.394e-14, -1.208e-14,
1780  4.44e-15, -1.68e-15,
1781  6.5e-16, -2.6e-16,
1782  1.1e-16, -4e-17,
1783  2e-17, -1e-17 };
1784
1785  const Double_t c0 = 2/TMath::Pi();
1786
1787  Int_t i;
1788  Double_t alfa, h, r, y, b0, b1, b2;
1789  Double_t v = TMath::Abs(x);
1790
1791  v = TMath::Abs(x);
1792  if (v < 8) {
1793  y = v/8;
1794  h = 2*y*y -1;
1795  alfa = h + h;
1796  b0 = 0;
1797  b1 = 0;
1798  b2 = 0;
1799  for (i = n1; i >= 0; --i) {
1800  b0 = c1[i] + alfa*b1 - b2;
1801  b2 = b1;
1802  b1 = b0;
1803  }
1804  h = y*(b0 - h*b2);
1805  } else {
1806  r = 1/v;
1807  h = 128*r*r -1;
1808  alfa = h + h;
1809  b0 = 0;
1810  b1 = 0;
1811  b2 = 0;
1812  for (i = n2; i >= 0; --i) {
1813  b0 = c2[i] + alfa*b1 - b2;
1814  b2 = b1;
1815  b1 = b0;
1816  }
1817  h = TMath::BesselY0(v) + r*c0*(b0 - h*b2);
1818  }
1819  if (x < 0) h = -h;
1820  return h;
1821 }
1822
1823 ////////////////////////////////////////////////////////////////////////////////
1824 /// Struve Functions of Order 1
1825 ///
1826 /// Converted from CERNLIB M342 by Rene Brun.
1827
1829 {
1830  const Int_t n3 = 16;
1831  const Int_t n4 = 22;
1832  const Double_t c3[17] = { .5578891446481605, -.11188325726569816,
1833  -.16337958125200939, .32256932072405902,
1834  -.14581632367244242, .03292677399374035,
1835  -.00460372142093573, 4.434706163314e-4,
1836  -3.142099529341e-5, 1.7123719938e-6,
1837  -7.416987005e-8, 2.61837671e-9,
1838  -7.685839e-11, 1.9067e-12,
1839  -4.052e-14, 7.5e-16,
1840  -1e-17 };
1841  const Double_t c4[23] = { 1.00757647293865641, .00750316051248257,
1842  -7.043933264519e-5, 2.66205393382e-6,
1843  -1.8841157753e-7, 1.949014958e-8,
1844  -2.6126199e-9, 4.236269e-10,
1845  -7.955156e-11, 1.679973e-11,
1846  -3.9072e-12, 9.8543e-13,
1847  -2.6636e-13, 7.645e-14,
1848  -2.313e-14, 7.33e-15,
1849  -2.42e-15, 8.3e-16,
1850  -3e-16, 1.1e-16,
1851  -4e-17, 2e-17,-1e-17 };
1852
1853  const Double_t c0 = 2/TMath::Pi();
1854  const Double_t cc = 2/(3*TMath::Pi());
1855
1856  Int_t i, i1;
1857  Double_t alfa, h, r, y, b0, b1, b2;
1858  Double_t v = TMath::Abs(x);
1859
1860  if (v == 0) {
1861  h = 0;
1862  } else if (v <= 0.3) {
1863  y = v*v;
1864  r = 1;
1865  h = 1;
1866  i1 = (Int_t)(-8. / TMath::Log10(v));
1867  for (i = 1; i <= i1; ++i) {
1868  h = -h*y / ((2*i+ 1)*(2*i + 3));
1869  r += h;
1870  }
1871  h = cc*y*r;
1872  } else if (v < 8) {
1873  h = v*v/32 -1;
1874  alfa = h + h;
1875  b0 = 0;
1876  b1 = 0;
1877  b2 = 0;
1878  for (i = n3; i >= 0; --i) {
1879  b0 = c3[i] + alfa*b1 - b2;
1880  b2 = b1;
1881  b1 = b0;
1882  }
1883  h = b0 - h*b2;
1884  } else {
1885  h = 128/(v*v) -1;
1886  alfa = h + h;
1887  b0 = 0;
1888  b1 = 0;
1889  b2 = 0;
1890  for (i = n4; i >= 0; --i) {
1891  b0 = c4[i] + alfa*b1 - b2;
1892  b2 = b1;
1893  b1 = b0;
1894  }
1895  h = TMath::BesselY1(v) + c0*(b0 - h*b2);
1896  }
1897  return h;
1898 }
1899
1900 ////////////////////////////////////////////////////////////////////////////////
1901 /// Modified Struve Function of Order 0.
1902 /// By Kirill Filimonov.
1903
1905 {
1906  const Double_t pi=TMath::Pi();
1907
1908  Double_t s=1.0;
1909  Double_t r=1.0;
1910
1911  Double_t a0,sl0,a1,bi0;
1912
1913  Int_t km,i;
1914
1915  if (x<=20.) {
1916  a0=2.0*x/pi;
1917  for (i=1; i<=60;i++) {
1918  r*=(x/(2*i+1))*(x/(2*i+1));
1919  s+=r;
1920  if(TMath::Abs(r/s)<1.e-12) break;
1921  }
1922  sl0=a0*s;
1923  } else {
1924  km=int(5*(x+1.0));
1925  if(x>=50.0)km=25;
1926  for (i=1; i<=km; i++) {
1927  r*=(2*i-1)*(2*i-1)/x/x;
1928  s+=r;
1929  if(TMath::Abs(r/s)<1.0e-12) break;
1930  }
1931  a1=TMath::Exp(x)/TMath::Sqrt(2*pi*x);
1932  r=1.0;
1933  bi0=1.0;
1934  for (i=1; i<=16; i++) {
1935  r=0.125*r*(2.0*i-1.0)*(2.0*i-1.0)/(i*x);
1936  bi0+=r;
1937  if(TMath::Abs(r/bi0)<1.0e-12) break;
1938  }
1939
1940  bi0=a1*bi0;
1941  sl0=-2.0/(pi*x)*s+bi0;
1942  }
1943  return sl0;
1944 }
1945
1946 ////////////////////////////////////////////////////////////////////////////////
1947 /// Modified Struve Function of Order 1.
1948 /// By Kirill Filimonov.
1949
1951 {
1952  const Double_t pi=TMath::Pi();
1953  Double_t a1,sl1,bi1,s;
1954  Double_t r=1.0;
1955  Int_t km,i;
1956
1957  if (x<=20.) {
1958  s=0.0;
1959  for (i=1; i<=60;i++) {
1960  r*=x*x/(4.0*i*i-1.0);
1961  s+=r;
1962  if(TMath::Abs(r)<TMath::Abs(s)*1.e-12) break;
1963  }
1964  sl1=2.0/pi*s;
1965  } else {
1966  s=1.0;
1967  km=int(0.5*x);
1968  if(x>50.0)km=25;
1969  for (i=1; i<=km; i++) {
1970  r*=(2*i+3)*(2*i+1)/x/x;
1971  s+=r;
1972  if(TMath::Abs(r/s)<1.0e-12) break;
1973  }
1974  sl1=2.0/pi*(-1.0+1.0/(x*x)+3.0*s/(x*x*x*x));
1975  a1=TMath::Exp(x)/TMath::Sqrt(2*pi*x);
1976  r=1.0;
1977  bi1=1.0;
1978  for (i=1; i<=16; i++) {
1979  r=-0.125*r*(4.0-(2.0*i-1.0)*(2.0*i-1.0))/(i*x);
1980  bi1+=r;
1981  if(TMath::Abs(r/bi1)<1.0e-12) break;
1982  }
1983  sl1+=a1*bi1;
1984  }
1985  return sl1;
1986 }
1987
1988 ////////////////////////////////////////////////////////////////////////////////
1989 /// Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
1990
1992 {
1994 }
1995
1996 ////////////////////////////////////////////////////////////////////////////////
1997 /// Continued fraction evaluation by modified Lentz's method
1998 /// used in calculation of incomplete Beta function.
1999
2001 {
2002  Int_t itmax = 500;
2003  Double_t eps = 3.e-14;
2004  Double_t fpmin = 1.e-30;
2005
2006  Int_t m, m2;
2007  Double_t aa, c, d, del, qab, qam, qap;
2008  Double_t h;
2009  qab = a+b;
2010  qap = a+1.0;
2011  qam = a-1.0;
2012  c = 1.0;
2013  d = 1.0 - qab*x/qap;
2014  if (TMath::Abs(d)<fpmin) d=fpmin;
2015  d=1.0/d;
2016  h=d;
2017  for (m=1; m<=itmax; m++) {
2018  m2=m*2;
2019  aa = m*(b-m)*x/((qam+ m2)*(a+m2));
2020  d = 1.0 +aa*d;
2021  if(TMath::Abs(d)<fpmin) d = fpmin;
2022  c = 1 +aa/c;
2023  if (TMath::Abs(c)<fpmin) c = fpmin;
2024  d=1.0/d;
2025  h*=d*c;
2026  aa = -(a+m)*(qab +m)*x/((a+m2)*(qap+m2));
2027  d=1.0+aa*d;
2028  if(TMath::Abs(d)<fpmin) d = fpmin;
2029  c = 1.0 +aa/c;
2030  if (TMath::Abs(c)<fpmin) c = fpmin;
2031  d=1.0/d;
2032  del = d*c;
2033  h*=del;
2034  if (TMath::Abs(del-1)<=eps) break;
2035  }
2036  if (m>itmax) {
2037  Info("TMath::BetaCf", "a or b too big, or itmax too small, a=%g, b=%g, x=%g, h=%g, itmax=%d",
2038  a,b,x,h,itmax);
2039  }
2040  return h;
2041 }
2042
2043 ////////////////////////////////////////////////////////////////////////////////
2044 /// Computes the probability density function of the Beta distribution
2045 /// (the distribution function is computed in BetaDistI).
2046 /// The first argument is the point, where the function will be
2047 /// computed, second and third are the function parameters.
2048 /// Since the Beta distribution is bounded on both sides, it's often
2049 /// used to represent processes with natural lower and upper limits.
2050
2052 {
2053  if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2054  Error("TMath::BetaDist", "parameter value outside allowed range");
2055  return 0;
2056  }
2057  Double_t beta = TMath::Beta(p, q);
2058  Double_t r = TMath::Power(x, p-1)*TMath::Power(1-x, q-1)/beta;
2059  return r;
2060 }
2061
2062 ////////////////////////////////////////////////////////////////////////////////
2063 /// Computes the distribution function of the Beta distribution.
2064 /// The first argument is the point, where the function will be
2065 /// computed, second and third are the function parameters.
2066 /// Since the Beta distribution is bounded on both sides, it's often
2067 /// used to represent processes with natural lower and upper limits.
2068
2070 {
2071  if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2072  Error("TMath::BetaDistI", "parameter value outside allowed range");
2073  return 0;
2074  }
2075  Double_t betai = TMath::BetaIncomplete(x, p, q);
2076  return betai;
2077 }
2078
2079 ////////////////////////////////////////////////////////////////////////////////
2080 /// Calculates the incomplete Beta-function.
2081
2083 {
2085 }
2086
2087 ////////////////////////////////////////////////////////////////////////////////
2088 /// Calculate the binomial coefficient n over k.
2089
2091 {
2092  if (n<0 || k<0 || n<k) return TMath::SignalingNaN();
2093  if (k==0 || n==k) return 1;
2094
2095  Int_t k1=TMath::Min(k,n-k);
2096  Int_t k2=n-k1;
2097  Double_t fact=k2+1;
2098  for (Double_t i=k1;i>1.;--i)
2099  fact *= (k2+i)/i;
2100  return fact;
2101 }
2102
2103 ////////////////////////////////////////////////////////////////////////////////
2104 /// Suppose an event occurs with probability _p_ per trial
2105 /// Then the probability P of its occurring _k_ or more times
2106 /// in _n_ trials is termed a cumulative binomial probability
2107 /// the formula is P = sum_from_j=k_to_n(TMath::Binomial(n, j)*
2108 /// *TMath::Power(p, j)*TMath::Power(1-p, n-j)
2109 /// For _n_ larger than 12 BetaIncomplete is a much better way
2110 /// to evaluate the sum than would be the straightforward sum calculation
2111 /// for _n_ smaller than 12 either method is acceptable
2112 /// ("Numerical Recipes")
2113 /// --implementation by Anna Kreshuk
2114
2116 {
2117  if(k <= 0) return 1.0;
2118  if(k > n) return 0.0;
2119  if(k == n) return TMath::Power(p, n);
2120
2121  return BetaIncomplete(p, Double_t(k), Double_t(n-k+1));
2122 }
2123
2124 ////////////////////////////////////////////////////////////////////////////////
2125 /// Computes the density of Cauchy distribution at point x
2126 /// by default, standard Cauchy distribution is used (t=0, s=1)
2127 /// - t is the location parameter
2128 /// - s is the scale parameter
2129 ///
2130 /// The Cauchy distribution, also called Lorentzian distribution,
2131 /// is a continuous distribution describing resonance behavior
2132 /// The mean and standard deviation of the Cauchy distribution are undefined.
2133 /// The practical meaning of this is that collecting 1,000 data points gives
2134 /// no more accurate an estimate of the mean and standard deviation than
2135 /// does a single point.
2136 /// The formula was taken from "Engineering Statistics Handbook" on site
2137 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3663.htm
2138 /// Implementation by Anna Kreshuk.
2139 ///
2140 /// Example:
2141 ///
2142 /// ~~~ {cpp}
2143 /// TF1* fc = new TF1("fc", "TMath::CauchyDist(x, [0], [1])", -5, 5);
2144 /// fc->SetParameters(0, 1);
2145 /// fc->Draw();
2146 /// ~~~
2147
2149 {
2150  Double_t temp= (x-t)*(x-t)/(s*s);
2151  Double_t result = 1/(s*TMath::Pi()*(1+temp));
2152  return result;
2153 }
2154
2155 ////////////////////////////////////////////////////////////////////////////////
2156 /// Evaluate the quantiles of the chi-squared probability distribution function.
2157 /// Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
2158 /// implemented by Anna Kreshuk.
2159 /// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
2160 ///
2161 /// \param[in] p the probability value, at which the quantile is computed
2162 /// \param[in] ndf number of degrees of freedom
2163
2165 {
2166  Double_t c[]={0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
2167  4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
2168  84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
2169  210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
2170  462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
2171  932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
2172  2520.0, 5040.0};
2173  Double_t e = 5e-7;
2174  Double_t aa = 0.6931471806;
2175  Int_t maxit = 20;
2176  Double_t ch, p1, p2, q, t, a, b, x;
2177  Double_t s1, s2, s3, s4, s5, s6;
2178
2179  if (ndf <= 0) return 0;
2180
2181  Double_t g = TMath::LnGamma(0.5*ndf);
2182
2183  Double_t xx = 0.5 * ndf;
2184  Double_t cp = xx - 1;
2185  if (ndf >= TMath::Log(p)*(-c[5])){
2186  //starting approximation for ndf less than or equal to 0.32
2187  if (ndf > c[3]) {
2188  x = TMath::NormQuantile(p);
2189  //starting approximation using Wilson and Hilferty estimate
2190  p1=c[2]/ndf;
2191  ch = ndf*TMath::Power((x*TMath::Sqrt(p1) + 1 - p1), 3);
2192  if (ch > c[6]*ndf + 6)
2193  ch = -2 * (TMath::Log(1-p) - cp * TMath::Log(0.5 * ch) + g);
2194  } else {
2195  ch = c[4];
2196  a = TMath::Log(1-p);
2197  do{
2198  q = ch;
2199  p1 = 1 + ch * (c[7]+ch);
2200  p2 = ch * (c[9] + ch * (c[8] + ch));
2201  t = -0.5 + (c[7] + 2 * ch) / p1 - (c[9] + ch * (c[10] + 3 * ch)) / p2;
2202  ch = ch - (1 - TMath::Exp(a + g + 0.5 * ch + cp * aa) *p2 / p1) / t;
2203  }while (TMath::Abs(q/ch - 1) > c[1]);
2204  }
2205  } else {
2206  ch = TMath::Power((p * xx * TMath::Exp(g + xx * aa)),(1./xx));
2207  if (ch < e) return ch;
2208  }
2209 //call to algorithm AS 239 and calculation of seven term Taylor series
2210  for (Int_t i=0; i<maxit; i++){
2211  q = ch;
2212  p1 = 0.5 * ch;
2213  p2 = p - TMath::Gamma(xx, p1);
2214
2215  t = p2 * TMath::Exp(xx * aa + g + p1 - cp * TMath::Log(ch));
2216  b = t / ch;
2217  a = 0.5 * t - b * cp;
2218  s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] +c[11] * a))))) / c[24];
2219  s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
2220  s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
2221  s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
2222  s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
2223  s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
2224  ch = ch + t * (1 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
2225  if (TMath::Abs(q / ch - 1) > e) break;
2226  }
2227  return ch;
2228 }
2229
2230 ////////////////////////////////////////////////////////////////////////////////
2231 /// Computes the density function of F-distribution
2232 /// (probability function, integral of density, is computed in FDistI).
2233 ///
2234 /// Parameters N and M stand for degrees of freedom of chi-squares
2235 /// mentioned above parameter F is the actual variable x of the
2236 /// density function p(x) and the point at which the density function
2237 /// is calculated.
2238 ///
2239 /// ### About F distribution:
2240 /// F-distribution arises in testing whether two random samples
2241 /// have the same variance. It is the ratio of two chi-square
2242 /// distributions, with N and M degrees of freedom respectively,
2243 /// where each chi-square is first divided by it's number of degrees
2244 /// of freedom.
2245 /// Implementation by Anna Kreshuk.
2246
2248 {
2250 }
2251
2252 ////////////////////////////////////////////////////////////////////////////////
2253 /// Calculates the cumulative distribution function of F-distribution,
2254 /// this function occurs in the statistical test of whether two observed
2255 /// samples have the same variance. For this test a certain statistic F,
2256 /// the ratio of observed dispersion of the first sample to that of the
2257 /// second sample, is calculated. N and M stand for numbers of degrees
2258 /// of freedom in the samples 1-FDistI() is the significance level at
2259 /// which the hypothesis "1 has smaller variance than 2" can be rejected.
2260 /// A small numerical value of 1 - FDistI() implies a very significant
2261 /// rejection, in turn implying high confidence in the hypothesis
2262 /// "1 has variance greater than 2".
2263 ///
2264 /// Implementation by Anna Kreshuk.
2265
2267 {
2269  return fi;
2270 }
2271
2272 ////////////////////////////////////////////////////////////////////////////////
2273 /// Computes the density function of Gamma distribution at point x.
2274 ///
2275 /// \param[in] gamma shape parameter
2276 /// \param[in] mu location parameter
2277 /// \param[in] beta scale parameter
2278 ///
2279 /// The definition can be found in "Engineering Statistics Handbook" on site
2280 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
2281 /// use now implementation in ROOT::Math::gamma_pdf
2282 ///
2283 /// Begin_Macro
2284 /// {
2285 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2286 ///
2287 /// c1->SetLogy();
2288 /// c1->SetGridx();
2289 /// c1->SetGridy();
2290 ///
2291 /// TF1 *gdist = new TF1("gdist", "TMath::GammaDist(x, [0], [1], [2])", 0, 10);
2292 ///
2293 /// gdist->SetParameters(0.5, 0., 1.);
2294 /// gdist->SetLineColor(2);
2295 /// TF1 *gdist1 = gdist->DrawCopy("L");
2296 /// gdist->SetParameters(1.0, 0., 1.);
2297 /// gdist->SetLineColor(3);
2298 /// TF1 *gdist2 = gdist->DrawCopy("LSAME");
2299 /// gdist->SetParameters(2.0, 0., 1.);
2300 /// gdist->SetLineColor(4);
2301 /// TF1 *gdist3 = gdist->DrawCopy("LSAME");
2302 /// gdist->SetParameters(5.0, 0., 1.);
2303 /// gdist->SetLineColor(6);
2304 /// TF1 *gdist4 = gdist->DrawCopy("LSAME");
2305 ///
2306 /// legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2307 /// legend->AddEntry(gdist1, "gamma = 0.5, mu = 0, beta = 1", "L");
2308 /// legend->AddEntry(gdist2, "gamma = 1.0, mu = 0, beta = 1", "L");
2309 /// legend->AddEntry(gdist3, "gamma = 2.0, mu = 0, beta = 1", "L");
2310 /// legend->AddEntry(gdist4, "gamma = 5.0, mu = 0, beta = 1", "L");
2311 /// legend->Draw();
2312 /// }
2313 /// End_Macro
2314
2316 {
2317  if ((x<mu) || (gamma<=0) || (beta <=0)) {
2318  Error("TMath::GammaDist", "illegal parameter values x = %f , gamma = %f beta = %f",x,gamma,beta);
2319  return 0;
2320  }
2321  return ::ROOT::Math::gamma_pdf(x, gamma, beta, mu);
2322 }
2323
2324 ////////////////////////////////////////////////////////////////////////////////
2325 /// Computes the probability density function of Laplace distribution
2326 /// at point x, with location parameter alpha and shape parameter beta.
2327 /// By default, alpha=0, beta=1
2328 /// This distribution is known under different names, most common is
2329 /// double exponential distribution, but it also appears as
2330 /// the two-tailed exponential or the bilateral exponential distribution
2331
2333 {
2334  Double_t temp;
2335  temp = TMath::Exp(-TMath::Abs((x-alpha)/beta));
2336  temp /= (2.*beta);
2337  return temp;
2338 }
2339
2340 ////////////////////////////////////////////////////////////////////////////////
2341 /// Computes the distribution function of Laplace distribution
2342 /// at point x, with location parameter alpha and shape parameter beta.
2343 /// By default, alpha=0, beta=1
2344 /// This distribution is known under different names, most common is
2345 /// double exponential distribution, but it also appears as
2346 /// the two-tailed exponential or the bilateral exponential distribution
2347
2349 {
2350  Double_t temp;
2351  if (x<=alpha){
2352  temp = 0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2353  } else {
2354  temp = 1-0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2355  }
2356  return temp;
2357 }
2358
2359 ////////////////////////////////////////////////////////////////////////////////
2360 /// Computes the density of LogNormal distribution at point x.
2361 /// Variable X has lognormal distribution if Y=Ln(X) has normal distribution
2362 ///
2363 /// \param[in] sigma is the shape parameter
2364 /// \param[in] theta is the location parameter
2365 /// \param[in] m is the scale parameter
2366 ///
2367 /// The formula was taken from "Engineering Statistics Handbook" on site
2368 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm
2369 /// Implementation using ROOT::Math::lognormal_pdf
2370 ///
2371 /// Begin_Macro
2372 /// {
2373 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2374 ///
2375 /// c1->SetLogy();
2376 /// c1->SetGridx();
2377 /// c1->SetGridy();
2378 ///
2379 /// TF1 *logn = new TF1("logn", "TMath::LogNormal(x, [0], [1], [2])", 0, 5);
2380 /// logn->SetMinimum(1e-3);
2381 ///
2382 /// logn->SetParameters(0.5, 0., 1.);
2383 /// logn->SetLineColor(2);
2384 /// TF1 *logn1 = logn->DrawCopy("L");
2385 /// logn->SetParameters(1.0, 0., 1.);
2386 /// logn->SetLineColor(3);
2387 /// TF1 *logn2 = logn->DrawCopy("LSAME");
2388 /// logn->SetParameters(2.0, 0., 1.);
2389 /// logn->SetLineColor(4);
2390 /// TF1 *logn3 = logn->DrawCopy("LSAME");
2391 /// logn->SetParameters(5.0, 0., 1.);
2392 /// logn->SetLineColor(6);
2393 /// TF1 *logn4 = logn->DrawCopy("LSAME");
2394 ///
2395 /// legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2396 /// legend->AddEntry(logn1, "sigma = 0.5, theta = 0, m = 1", "L");
2397 /// legend->AddEntry(logn2, "sigma = 1.0, theta = 0, m = 1", "L");
2398 /// legend->AddEntry(logn3, "sigma = 2.0, theta = 0, m = 1", "L");
2399 /// legend->AddEntry(logn4, "sigma = 5.0, theta = 0, m = 1", "L");
2400 /// legend->Draw();
2401 /// }
2402 /// End_Macro
2403
2405 {
2406  if ((x<theta) || (sigma<=0) || (m<=0)) {
2407  Error("TMath::Lognormal", "illegal parameter values");
2408  return 0;
2409  }
2410  // lognormal_pdf uses same definition of http://en.wikipedia.org/wiki/Log-normal_distribution
2411  // where mu = log(m)
2412
2413  return ::ROOT::Math::lognormal_pdf(x, TMath::Log(m), sigma, theta);
2414
2415 }
2416
2417 ////////////////////////////////////////////////////////////////////////////////
2418 /// Computes quantiles for standard normal distribution N(0, 1)
2419 /// at probability p
2420 ///
2421 /// ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.
2422
2424 {
2425  if ((p<=0)||(p>=1)) {
2426  Error("TMath::NormQuantile", "probability outside (0, 1)");
2427  return 0;
2428  }
2429
2430  Double_t a0 = 3.3871328727963666080e0;
2431  Double_t a1 = 1.3314166789178437745e+2;
2432  Double_t a2 = 1.9715909503065514427e+3;
2433  Double_t a3 = 1.3731693765509461125e+4;
2434  Double_t a4 = 4.5921953931549871457e+4;
2435  Double_t a5 = 6.7265770927008700853e+4;
2436  Double_t a6 = 3.3430575583588128105e+4;
2437  Double_t a7 = 2.5090809287301226727e+3;
2438  Double_t b1 = 4.2313330701600911252e+1;
2439  Double_t b2 = 6.8718700749205790830e+2;
2440  Double_t b3 = 5.3941960214247511077e+3;
2441  Double_t b4 = 2.1213794301586595867e+4;
2442  Double_t b5 = 3.9307895800092710610e+4;
2443  Double_t b6 = 2.8729085735721942674e+4;
2444  Double_t b7 = 5.2264952788528545610e+3;
2445  Double_t c0 = 1.42343711074968357734e0;
2446  Double_t c1 = 4.63033784615654529590e0;
2447  Double_t c2 = 5.76949722146069140550e0;
2448  Double_t c3 = 3.64784832476320460504e0;
2449  Double_t c4 = 1.27045825245236838258e0;
2450  Double_t c5 = 2.41780725177450611770e-1;
2451  Double_t c6 = 2.27238449892691845833e-2;
2452  Double_t c7 = 7.74545014278341407640e-4;
2453  Double_t d1 = 2.05319162663775882187e0;
2454  Double_t d2 = 1.67638483018380384940e0;
2455  Double_t d3 = 6.89767334985100004550e-1;
2456  Double_t d4 = 1.48103976427480074590e-1;
2457  Double_t d5 = 1.51986665636164571966e-2;
2458  Double_t d6 = 5.47593808499534494600e-4;
2459  Double_t d7 = 1.05075007164441684324e-9;
2460  Double_t e0 = 6.65790464350110377720e0;
2461  Double_t e1 = 5.46378491116411436990e0;
2462  Double_t e2 = 1.78482653991729133580e0;
2463  Double_t e3 = 2.96560571828504891230e-1;
2464  Double_t e4 = 2.65321895265761230930e-2;
2465  Double_t e5 = 1.24266094738807843860e-3;
2466  Double_t e6 = 2.71155556874348757815e-5;
2467  Double_t e7 = 2.01033439929228813265e-7;
2468  Double_t f1 = 5.99832206555887937690e-1;
2469  Double_t f2 = 1.36929880922735805310e-1;
2470  Double_t f3 = 1.48753612908506148525e-2;
2471  Double_t f4 = 7.86869131145613259100e-4;
2472  Double_t f5 = 1.84631831751005468180e-5;
2473  Double_t f6 = 1.42151175831644588870e-7;
2474  Double_t f7 = 2.04426310338993978564e-15;
2475
2476  Double_t split1 = 0.425;
2477  Double_t split2=5.;
2478  Double_t konst1=0.180625;
2479  Double_t konst2=1.6;
2480
2481  Double_t q, r, quantile;
2482  q=p-0.5;
2483  if (TMath::Abs(q)<split1) {
2484  r=konst1-q*q;
2485  quantile = q* (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
2486  * r + a2) * r + a1) * r + a0) /
2487  (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
2488  * r + b2) * r + b1) * r + 1.);
2489  } else {
2490  if(q<0) r=p;
2491  else r=1-p;
2492  //error case
2493  if (r<=0)
2494  quantile=0;
2495  else {
2496  r=TMath::Sqrt(-TMath::Log(r));
2497  if (r<=split2) {
2498  r=r-konst2;
2499  quantile=(((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
2500  * r + c2) * r + c1) * r + c0) /
2501  (((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
2502  * r + d2) * r + d1) * r + 1);
2503  } else{
2504  r=r-split2;
2505  quantile=(((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
2506  * r + e2) * r + e1) * r + e0) /
2507  (((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
2508  * r + f2) * r + f1) * r + 1);
2509  }
2510  if (q<0) quantile=-quantile;
2511  }
2512  }
2513  return quantile;
2514 }
2515
2516 ////////////////////////////////////////////////////////////////////////////////
2517 /// Simple recursive algorithm to find the permutations of
2518 /// n natural numbers, not necessarily all distinct
2519 /// adapted from CERNLIB routine PERMU.
2520 /// The input array has to be initialised with a non descending
2521 /// sequence. The method returns kFALSE when
2522 /// all combinations are exhausted.
2523
2525 {
2526  Int_t i,itmp;
2527  Int_t i1=-1;
2528
2529  // find rightmost upward transition
2530  for(i=n-2; i>-1; i--) {
2531  if(a[i]<a[i+1]) {
2532  i1=i;
2533  break;
2534  }
2535  }
2536  // no more upward transitions, end of the story
2537  if(i1==-1) return kFALSE;
2538  else {
2539  // find lower right element higher than the lower
2540  // element of the upward transition
2541  for(i=n-1;i>i1;i--) {
2542  if(a[i] > a[i1]) {
2543  // swap the two
2544  itmp=a[i1];
2545  a[i1]=a[i];
2546  a[i]=itmp;
2547  break;
2548  }
2549  }
2550  // order the rest, in fact just invert, as there
2551  // are only downward transitions from here on
2552  for(i=0;i<(n-i1-1)/2;i++) {
2553  itmp=a[i1+i+1];
2554  a[i1+i+1]=a[n-i-1];
2555  a[n-i-1]=itmp;
2556  }
2557  }
2558  return kTRUE;
2559 }
2560
2561 ////////////////////////////////////////////////////////////////////////////////
2562 /// Computes density function for Student's t- distribution
2563 /// (the probability function (integral of density) is computed in StudentI).
2564 ///
2565 /// First parameter stands for x - the actual variable of the
2566 /// density function p(x) and the point at which the density is calculated.
2567 /// Second parameter stands for number of degrees of freedom.
2568 ///
2570 /// Student's t-distribution is used for many significance tests, for example,
2571 /// for the Student's t-tests for the statistical significance of difference
2572 /// between two sample means and for confidence intervals for the difference
2573 /// between two population means.
2574 ///
2575 /// Example: suppose we have a random sample of size n drawn from normal
2576 /// distribution with mean Mu and st.deviation Sigma. Then the variable
2577 ///
2578 /// t = (sample_mean - Mu)/(sample_deviation / sqrt(n))
2579 ///
2580 /// has Student's t-distribution with n-1 degrees of freedom.
2581 ///
2582 /// NOTE that this function's second argument is number of degrees of freedom,
2583 /// not the sample size.
2584 ///
2585 /// As the number of degrees of freedom grows, t-distribution approaches
2586 /// Normal(0,1) distribution.
2587 ///
2588 /// Implementation by Anna Kreshuk.
2589
2591 {
2592  if (ndf < 1) {
2593  return 0;
2594  }
2595
2596  Double_t r = ndf;
2597  Double_t rh = 0.5*r;
2598  Double_t rh1 = rh + 0.5;
2599  Double_t denom = TMath::Sqrt(r*TMath::Pi())*TMath::Gamma(rh)*TMath::Power(1+T*T/r, rh1);
2600  return TMath::Gamma(rh1)/denom;
2601 }
2602
2603 ////////////////////////////////////////////////////////////////////////////////
2604 /// Calculates the cumulative distribution function of Student's
2605 /// t-distribution second parameter stands for number of degrees of freedom,
2606 /// not for the number of samples
2607 /// if x has Student's t-distribution, the function returns the probability of
2608 /// x being less than T.
2609 ///
2610 /// Implementation by Anna Kreshuk.
2611
2613 {
2614  Double_t r = ndf;
2615
2616  Double_t si = (T>0) ?
2617  (1 - 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5)) :
2618  0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5);
2619  return si;
2620 }
2621
2622 ////////////////////////////////////////////////////////////////////////////////
2623 /// Computes quantiles of the Student's t-distribution
2624 /// 1st argument is the probability, at which the quantile is computed
2625 /// 2nd argument - the number of degrees of freedom of the
2626 /// Student distribution
2627 /// When the 3rd argument lower_tail is kTRUE (default)-
2628 /// the algorithm returns such x0, that
2629 ///
2630 /// P(x < x0)=p
2631 ///
2632 /// upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that
2633 ///
2634 /// P(x > x0)=p
2635 ///
2636 /// the algorithm was taken from:
2637 /// G.W.Hill, "Algorithm 396, Student's t-quantiles"
2638 /// "Communications of the ACM", 13(10), October 1970
2639
2641 {
2642  Double_t quantile;
2643  Double_t temp;
2644  Bool_t neg;
2645  Double_t q;
2646  if (ndf<1 || p>=1 || p<=0) {
2647  Error("TMath::StudentQuantile", "illegal parameter values");
2648  return 0;
2649  }
2650  if ((lower_tail && p>0.5)||(!lower_tail && p<0.5)){
2651  neg=kFALSE;
2652  q=2*(lower_tail ? (1-p) : p);
2653  } else {
2654  neg=kTRUE;
2655  q=2*(lower_tail? p : (1-p));
2656  }
2657
2658  if ((ndf-1)<1e-8) {
2659  temp=TMath::PiOver2()*q;
2660  quantile = TMath::Cos(temp)/TMath::Sin(temp);
2661  } else {
2662  if ((ndf-2)<1e-8) {
2663  quantile = TMath::Sqrt(2./(q*(2-q))-2);
2664  } else {
2665  Double_t a=1./(ndf-0.5);
2666  Double_t b=48./(a*a);
2667  Double_t c=((20700*a/b -98)*a-16)*a+96.36;
2668  Double_t d=((94.5/(b+c)-3.)/b+1)*TMath::Sqrt(a*TMath::PiOver2())*ndf;
2669  Double_t x=q*d;
2670  Double_t y=TMath::Power(x, (2./ndf));
2671  if (y>0.05+a){
2672  //asymptotic inverse expansion about normal
2673  x=TMath::NormQuantile(q*0.5);
2674  y=x*x;
2675  if (ndf<5) c+=0.3*(ndf-4.5)*(x+0.6);
2676  c+=(((0.05*d*x-5.)*x-7.)*x-2.)*x +b;
2677  y=(((((0.4*y+6.3)*y+36.)*y+94.5)/c - y-3.)/b+1)*x;
2678  y=a*y*y;
2679  if(y>0.002) y = TMath::Exp(y)-1;
2680  else y += 0.5*y*y;
2681  } else {
2682  y=((1./(((ndf+6.)/(ndf*y)-0.089*d-0.822)*(ndf+2.)*3)+0.5/(ndf+4.))*y-1.)*
2683  (ndf+1.)/(ndf+2.)+1/y;
2684  }
2685  quantile = TMath::Sqrt(ndf*y);
2686  }
2687  }
2688  if(neg) quantile=-quantile;
2689  return quantile;
2690 }
2691
2692 ////////////////////////////////////////////////////////////////////////////////
2693 /// Returns the value of the Vavilov density function
2694 ///
2695 /// \param[in] x the point were the density function is evaluated
2696 /// \param[in] kappa value of kappa (distribution parameter)
2697 /// \param[in] beta2 value of beta2 (distribution parameter)
2698 ///
2699 /// The algorithm was taken from the CernLib function vavden(G115)
2700 /// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2701 /// Nucl.Instr. and Meth. B47(1990), 215-224
2702 ///
2703 /// Accuracy: quote from the reference above:
2704 ///
2705 /// "The results of our code have been compared with the values of the Vavilov
2706 /// density function computed numerically in an accurate way: our approximation
2707 /// shows a difference of less than 3% around the peak of the density function, slowly
2708 /// increasing going towards the extreme tails to the right and to the left"
2709 ///
2710 /// Begin_Macro
2711 /// {
2712 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2713 ///
2714 /// c1->SetGridx();
2715 /// c1->SetGridy();
2716 ///
2717 /// TF1 *vavilov = new TF1("vavilov", "TMath::Vavilov(x, [0], [1])", -3, 11);
2718 ///
2719 /// vavilov->SetParameters(0.5, 0.);
2720 /// vavilov->SetLineColor(2);
2721 /// TF1 *vavilov1 = vavilov->DrawCopy("L");
2722 /// vavilov->SetParameters(0.3, 0.);
2723 /// vavilov->SetLineColor(3);
2724 /// TF1 *vavilov2 = vavilov->DrawCopy("LSAME");
2725 /// vavilov->SetParameters(0.2, 0.);
2726 /// vavilov->SetLineColor(4);
2727 /// TF1 *vavilov3 = vavilov->DrawCopy("LSAME");
2728 /// vavilov->SetParameters(0.1, 0.);
2729 /// vavilov->SetLineColor(6);
2730 /// TF1 *vavilov4 = vavilov->DrawCopy("LSAME");
2731 ///
2732 /// legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2733 /// legend->AddEntry(vavilov1, "kappa = 0.5, beta2 = 0", "L");
2734 /// legend->AddEntry(vavilov2, "kappa = 0.3, beta2 = 0", "L");
2735 /// legend->AddEntry(vavilov3, "kappa = 0.2, beta2 = 0", "L");
2736 /// legend->AddEntry(vavilov4, "kappa = 0.1, beta2 = 0", "L");
2737 /// legend->Draw();
2738 /// }
2739 /// End_Macro
2740
2742 {
2743  Double_t *ac = new Double_t[14];
2744  Double_t *hc = new Double_t[9];
2745
2746  Int_t itype;
2747  Int_t npt;
2748  TMath::VavilovSet(kappa, beta2, 0, 0, ac, hc, itype, npt);
2749  Double_t v = TMath::VavilovDenEval(x, ac, hc, itype);
2750  delete [] ac;
2751  delete [] hc;
2752  return v;
2753 }
2754
2755 ////////////////////////////////////////////////////////////////////////////////
2756 /// Returns the value of the Vavilov distribution function
2757 ///
2758 /// \param[in] x the point were the density function is evaluated
2759 /// \param[in] kappa value of kappa (distribution parameter)
2760 /// \param[in] beta2 value of beta2 (distribution parameter)
2761 ///
2762 /// The algorithm was taken from the CernLib function vavden(G115)
2763 ///
2764 /// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2765 /// Nucl.Instr. and Meth. B47(1990), 215-224
2766 ///
2767 /// Accuracy: quote from the reference above:
2768 ///
2769 /// "The results of our code have been compared with the values of the Vavilov
2770 /// density function computed numerically in an accurate way: our approximation
2771 /// shows a difference of less than 3% around the peak of the density function, slowly
2772 /// increasing going towards the extreme tails to the right and to the left"
2773
2775 {
2776  Double_t *ac = new Double_t[14];
2777  Double_t *hc = new Double_t[9];
2778  Double_t *wcm = new Double_t[201];
2779  Int_t itype;
2780  Int_t npt;
2781  Int_t k;
2782  Double_t xx, v;
2783  TMath::VavilovSet(kappa, beta2, 1, wcm, ac, hc, itype, npt);
2784  if (x < ac[0]) v = 0;
2785  else if (x >=ac[8]) v = 1;
2786  else {
2787  xx = x - ac[0];
2788  k = Int_t(xx*ac[10]);
2789  v = TMath::Min(wcm[k] + (xx - k*ac[9])*(wcm[k+1]-wcm[k])*ac[10], 1.);
2790  }
2791  delete [] ac;
2792  delete [] hc;
2793  delete [] wcm;
2794  return v;
2795 }
2796
2797 ////////////////////////////////////////////////////////////////////////////////
2798 /// Returns the value of the Landau distribution function at point x.
2799 /// The algorithm was taken from the Cernlib function dislan(G110)
2800 /// Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
2801 /// distribution", Computer Phys.Comm., 31(1984), 97-111
2802
2804 {
2806 }
2807
2808
2809 ////////////////////////////////////////////////////////////////////////////////
2810 /// Internal function, called by Vavilov and VavilovI
2811
2812 void TMath::VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt)
2813 {
2814
2815  Double_t BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
2816  BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
2817  BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
2818
2819  Double_t FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
2820  FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
2821  FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
2822
2823  Double_t FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
2824
2825  Double_t EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
2826  0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
2827
2828  Double_t U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
2829  0. , 0.24880692e-1, 0.47404356e-2,
2830  -0.74445130e-3, 0.73225731e-2, 0. ,
2831  0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
2832
2833  Double_t U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
2834  0. , 0.42127077e-1, 0.73167928e-2,
2835  -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
2836  0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
2837
2838  Double_t U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
2839  0. , 0.23834176e-1, 0.21624675e-2,
2840  -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
2841  0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
2842
2843  Double_t U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
2844  -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
2845  0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
2846  -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
2847
2848  Double_t U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
2849  0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
2850  0.48736023e-3, 0.34850854e-2, 0. ,
2851  -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
2852
2853  Double_t U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
2854  -0.30188807e-2, -0.84479939e-3, 0. ,
2855  0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
2856  -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
2857
2858  Double_t U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
2859  -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
2860  0. , 0.50505298e+0};
2861  Double_t U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
2862  -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
2863  -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
2864  0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
2865
2866  Double_t V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
2867  0. , 0.45091424e-1, 0.80559636e-2,
2868  -0.38974523e-2, 0. , -0.30634124e-2,
2869  0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
2870
2871  Double_t V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
2872  0. , 0.12693873e+0, 0.22999801e-1,
2873  -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
2874  0. , 0.19716857e-1, 0.32596742e-2};
2875
2876  Double_t V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
2877  -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
2878  -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
2879  0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
2880
2881  Double_t V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
2882  0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
2883  0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
2884  -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
2885
2886  Double_t V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
2887  -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
2888  0.56856517e-2, -0.13438433e-1, 0. ,
2889  0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
2890
2891  Double_t V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
2892  0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
2893  0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
2894  -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
2895
2896  Double_t V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
2897  0. , -0.12218009e+1, -0.32324120e+0,
2898  -0.27373591e-1, 0.12173464e+0, 0. ,
2899  0. , 0.40917471e-1};
2900
2901  Double_t V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
2902  0. , -0.18160479e+1, -0.50919193e+0,
2903  -0.51384654e-1, 0.21413992e+0, 0. ,
2904  0. , 0.66596366e-1};
2905
2906  Double_t W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
2907  -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
2908  -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
2909  0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
2910
2911  Double_t W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
2912  -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
2913  0. , 0.23020158e-1, 0.50574313e-2,
2914  0.94550140e-2, 0.19300232e-1};
2915
2916  Double_t W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
2917  -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
2918  -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
2919  0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
2920
2921  Double_t W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
2922  0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
2923  0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
2924  -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
2925
2926  Double_t W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
2927  0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
2928  0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
2929  -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
2930
2931  Double_t W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
2932  0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
2933  0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
2934  -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
2935
2936  Double_t W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
2937  0. , -0.14540925e+1, -0.39529833e+0,
2938  -0.44293243e-1, 0.88741049e-1};
2939
2940  itype = 0;
2941  if (rkappa <0.01 || rkappa >12) {
2942  Error("Vavilov distribution", "illegal value of kappa");
2943  return;
2944  }
2945
2946  Double_t DRK[6];
2947  Double_t DSIGM[6];
2948  Double_t ALFA[8];
2949  Int_t j;
2950  Double_t x, y, xx, yy, x2, x3, y2, y3, xy, p2, p3, q2, q3, pq;
2951  if (rkappa >= 0.29) {
2952  itype = 1;
2953  npt = 100;
2954  Double_t wk = 1./TMath::Sqrt(rkappa);
2955
2956  AC[0] = (-0.032227*beta2-0.074275)*rkappa + (0.24533*beta2+0.070152)*wk + (-0.55610*beta2-3.1579);
2957  AC[8] = (-0.013483*beta2-0.048801)*rkappa + (-1.6921*beta2+8.3656)*wk + (-0.73275*beta2-3.5226);
2958  DRK[1] = wk*wk;
2959  DSIGM[1] = TMath::Sqrt(rkappa/(1-0.5*beta2));
2960  for (j=1; j<=4; j++) {
2961  DRK[j+1] = DRK[1]*DRK[j];
2962  DSIGM[j+1] = DSIGM[1]*DSIGM[j];
2963  ALFA[j+1] = (FNINV[j]-beta2*FNINV[j+1])*DRK[j];
2964  }
2965  HC[0]=TMath::Log(rkappa)+beta2+0.42278434;
2966  HC[1]=DSIGM[1];
2967  HC[2]=ALFA[3]*DSIGM[3];
2968  HC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
2969  HC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*HC[2];
2970  HC[5]=HC[2]*HC[2];
2971  HC[6]=HC[2]*HC[3];
2972  HC[7]=HC[2]*HC[5];
2973  for (j=2; j<=7; j++)
2974  HC[j]*=EDGEC[j];
2975  HC[8]=0.39894228*HC[1];
2976  }
2977  else if (rkappa >=0.22) {
2978  itype = 2;
2979  npt = 150;
2980  x = 1+(rkappa-BKMXX3)*FBKX3;
2981  y = 1+(TMath::Sqrt(beta2)-BKMXY3)*FBKY3;
2982  xx = 2*x;
2983  yy = 2*y;
2984  x2 = xx*x-1;
2985  x3 = xx*x2-x;
2986  y2 = yy*y-1;
2987  y3 = yy*y2-y;
2988  xy = x*y;
2989  p2 = x2*y;
2990  p3 = x3*y;
2991  q2 = y2*x;
2992  q3 = y3*x;
2993  pq = x2*y2;
2994  AC[1] = W1[1] + W1[2]*x + W1[4]*x3 + W1[5]*y + W1[6]*y2 + W1[7]*y3 +
2995  W1[8]*xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
2996  AC[2] = W2[1] + W2[2]*x + W2[3]*x2 + W2[4]*x3 + W2[5]*y + W2[6]*y2 +
2997  W2[8]*xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
2998  AC[3] = W3[1] + W3[3]*x2 + W3[4]*x3 + W3[5]*y + W3[6]*y2 + W3[7]*y3 +
2999  W3[8]*xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
3000  AC[4] = W4[1] + W4[2]*x + W4[3]*x2 + W4[4]*x3 + W4[5]*y + W4[6]*y2 + W4[7]*y3 +
3001  W4[8]*xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
3002  AC[5] = W5[1] + W5[2]*x + W5[4]*x3 + W5[5]*y + W5[6]*y2 + W5[7]*y3 +
3003  W5[8]*xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
3004  AC[6] = W6[1] + W6[2]*x + W6[3]*x2 + W6[4]*x3 + W6[5]*y + W6[6]*y2 + W6[7]*y3 +
3005  W6[8]*xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
3006  AC[8] = W8[1] + W8[2]*x + W8[3]*x2 + W8[5]*y + W8[6]*y2 + W8[7]*y3 + W8[8]*xy;
3007  AC[0] = -3.05;
3008  } else if (rkappa >= 0.12) {
3009  itype = 3;
3010  npt = 200;
3011  x = 1 + (rkappa-BKMXX2)*FBKX2;
3012  y = 1 + (TMath::Sqrt(beta2)-BKMXY2)*FBKY2;
3013  xx = 2*x;
3014  yy = 2*y;
3015  x2 = xx*x-1;
3016  x3 = xx*x2-x;
3017  y2 = yy*y-1;
3018  y3 = yy*y2-y;
3019  xy = x*y;
3020  p2 = x2*y;
3021  p3 = x3*y;
3022  q2 = y2*x;
3023  q3 = y3*x;
3024  pq = x2*y2;
3025  AC[1] = V1[1] + V1[2]*x + V1[3]*x2 + V1[5]*y + V1[6]*y2 + V1[7]*y3 +
3026  V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
3027  AC[2] = V2[1] + V2[2]*x + V2[3]*x2 + V2[5]*y + V2[6]*y2 + V2[7]*y3 +
3028  V2[8]*xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
3029  AC[3] = V3[1] + V3[2]*x + V3[3]*x2 + V3[4]*x3 + V3[5]*y + V3[6]*y2 + V3[7]*y3 +
3030  V3[8]*xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
3031  AC[4] = V4[1] + V4[2]*x + V4[3]*x2 + V4[4]*x3 + V4[5]*y + V4[6]*y2 + V4[7]*y3 +
3032  V4[8]*xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
3033  AC[5] = V5[1] + V5[2]*x + V5[3]*x2 + V5[4]*x3 + V5[5]*y + V5[6]*y2 + V5[7]*y3 +
3034  V5[8]*xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
3035  AC[6] = V6[1] + V6[2]*x + V6[3]*x2 + V6[4]*x3 + V6[5]*y + V6[6]*y2 + V6[7]*y3 +
3036  V6[8]*xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
3037  AC[7] = V7[1] + V7[2]*x + V7[3]*x2 + V7[5]*y + V7[6]*y2 + V7[7]*y3 +
3038  V7[8]*xy + V7[11]*q2;
3039  AC[8] = V8[1] + V8[2]*x + V8[3]*x2 + V8[5]*y + V8[6]*y2 + V8[7]*y3 +
3040  V8[8]*xy + V8[11]*q2;
3041  AC[0] = -3.04;
3042  } else {
3043  itype = 4;
3044  if (rkappa >=0.02) itype = 3;
3045  npt = 200;
3046  x = 1+(rkappa-BKMXX1)*FBKX1;
3047  y = 1+(TMath::Sqrt(beta2)-BKMXY1)*FBKY1;
3048  xx = 2*x;
3049  yy = 2*y;
3050  x2 = xx*x-1;
3051  x3 = xx*x2-x;
3052  y2 = yy*y-1;
3053  y3 = yy*y2-y;
3054  xy = x*y;
3055  p2 = x2*y;
3056  p3 = x3*y;
3057  q2 = y2*x;
3058  q3 = y3*x;
3059  pq = x2*y2;
3060  if (itype==3){
3061  AC[1] = U1[1] + U1[2]*x + U1[3]*x2 + U1[5]*y + U1[6]*y2 + U1[7]*y3 +
3062  U1[8]*xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
3063  AC[2] = U2[1] + U2[2]*x + U2[3]*x2 + U2[5]*y + U2[6]*y2 + U2[7]*y3 +
3064  U2[8]*xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
3065  AC[3] = U3[1] + U3[2]*x + U3[3]*x2 + U3[5]*y + U3[6]*y2 + U3[7]*y3 +
3066  U3[8]*xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
3067  AC[4] = U4[1] + U4[2]*x + U4[3]*x2 + U4[4]*x3 + U4[5]*y + U4[6]*y2 + U4[7]*y3 +
3068  U4[8]*xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
3069  AC[5] = U5[1] + U5[2]*x + U5[3]*x2 + U5[4]*x3 + U5[5]*y + U5[6]*y2 + U5[7]*y3 +
3070  U5[8]*xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
3071  AC[6] = U6[1] + U6[2]*x + U6[3]*x2 + U6[4]*x3 + U6[5]*y + U6[7]*y3 +
3072  U6[8]*xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
3073  AC[7] = U7[1] + U7[2]*x + U7[3]*x2 + U7[4]*x3 + U7[5]*y + U7[6]*y2 + U7[8]*xy;
3074  }
3075  AC[8] = U8[1] + U8[2]*x + U8[3]*x2 + U8[4]*x3 + U8[5]*y + U8[6]*y2 + U8[7]*y3 +
3076  U8[8]*xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
3077  AC[0] = -3.03;
3078  }
3079
3080  AC[9] = (AC[8] - AC[0])/npt;
3081  AC[10] = 1./AC[9];
3082  if (itype == 3) {
3083  x = (AC[7]-AC[8])/(AC[7]*AC[8]);
3084  y = 1./TMath::Log(AC[8]/AC[7]);
3085  p2 = AC[7]*AC[7];
3086  AC[11] = p2*(AC[1]*TMath::Exp(-AC[2]*(AC[7]+AC[5]*p2)-
3087  AC[3]*TMath::Exp(-AC[4]*(AC[7]+AC[6]*p2)))-0.045*y/AC[7])/(1+x*y*AC[7]);
3088  AC[12] = (0.045+x*AC[11])*y;
3089  }
3090  if (itype == 4) AC[13] = 0.995/LandauI(AC[8]);
3091
3092  if (mode==0) return;
3093  //
3094  x = AC[0];
3095  WCM[0] = 0;
3096  Double_t fl, fu;
3097  Int_t k;
3098  fl = TMath::VavilovDenEval(x, AC, HC, itype);
3099  for (k=1; k<=npt; k++) {
3100  x += AC[9];
3101  fu = TMath::VavilovDenEval(x, AC, HC, itype);
3102  WCM[k] = WCM[k-1] + fl + fu;
3103  fl = fu;
3104  }
3105  x = 0.5*AC[9];
3106  for (k=1; k<=npt; k++)
3107  WCM[k]*=x;
3108 }
3109
3110 ////////////////////////////////////////////////////////////////////////////////
3111 /// Internal function, called by Vavilov and VavilovSet
3112
3114 {
3115  Double_t v = 0;
3116  if (rlam < AC[0] || rlam > AC[8])
3117  return 0;
3118  Int_t k;
3119  Double_t x, fn, s;
3120  Double_t h[10];
3121  if (itype ==1 ) {
3122  fn = 1;
3123  x = (rlam + HC[0])*HC[1];
3124  h[1] = x;
3125  h[2] = x*x -1;
3126  for (k=2; k<=8; k++) {
3127  fn++;
3128  h[k+1] = x*h[k]-fn*h[k-1];
3129  }
3130  s = 1 + HC[7]*h[9];
3131  for (k=2; k<=6; k++)
3132  s+=HC[k]*h[k+1];
3133  v = HC[8]*TMath::Exp(-0.5*x*x)*TMath::Max(s, 0.);
3134  }
3135  else if (itype == 2) {
3136  x = rlam*rlam;
3137  v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x) - AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3138  }
3139  else if (itype == 3) {
3140  if (rlam < AC[7]) {
3141  x = rlam*rlam;
3142  v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x)-AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3143  } else {
3144  x = 1./rlam;
3145  v = (AC[11]*x + AC[12])*x;
3146  }
3147  }
3148  else if (itype == 4) {
3149  v = AC[13]*TMath::Landau(rlam);
3150  }
3151  return v;
3152 }
3153
3154
3155 //explicitly instantiate template functions from VecCore
3156 #ifdef R__HAS_VECCORE
3157 #include <Math/Types.h>
3163 // missing in veccore
3164 // template ROOT::Double_v vecCore::math::ACos(const ROOT::Double_v & x);
3165 // template ROOT::Double_v vecCore::math::Sinh(const ROOT::Double_v & x);
3166 // template ROOT::Double_v vecCore::math::Cosh(const ROOT::Double_v & x);
3167 // template ROOT::Double_v vecCore::math::Tanh(const ROOT::Double_v & x);
3168 // template ROOT::Double_v vecCore::math::Cbrt(const ROOT::Double_v & x);
3169 #endif
double lognormal_pdf(double x, double m, double s, double x0=0)
Probability density function of the lognormal distribution.
Double_t ACosH(Double_t)
Definition: TMath.cxx:77
Double_t BesselI(Int_t n, Double_t x)
constexpr Double_t HC()
in
Definition: TMath.h:232
Double_t FDist(Double_t F, Double_t N, Double_t M)
double erf(double x)
Error function encountered in integrating the normal distribution.
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
Definition: TMath.cxx:469
Double_t ErfInverse(Double_t x)
returns the inverse error function x must be <-1<x<1
Definition: TMath.cxx:203
static long int sum(long int i)
Definition: Factory.cxx:2258
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculate a Breit Wigner function with mean and gamma.
Definition: TMath.cxx:437
Double_t FDistI(Double_t F, Double_t N, Double_t M)
Double_t PoissonI(Double_t x, Double_t par)
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition: TMath.cxx:495
static constexpr double pi
auto * m
Definition: textangle.C:8
static double p3(double t, double a, double b, double c, double d)
T1 Sign(T1 a, T2 b)
Definition: TMathBase.h:153
Double_t Log(Double_t x)
Definition: TMath.h:759
float Float_t
Definition: RtypesCore.h:53
const char Option_t
Definition: RtypesCore.h:62
Double_t KolmogorovProb(Double_t z)
return c1
Definition: legend1.C:41
double T(double x)
Definition: ChebyshevPol.h:34
#define g(i)
Definition: RSha256.hxx:105
Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1)
Double_t StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail=kTRUE)
Double_t NormQuantile(Double_t p)
static constexpr double km
Double_t BetaDist(Double_t x, Double_t p, Double_t q)
Double_t LaplaceDist(Double_t x, Double_t alpha=0, Double_t beta=1)
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1113
double inc_beta(double x, double a, double b)
Calculates the normalized (regularized) incomplete beta function.
Double_t Log2(Double_t x)
Definition: TMath.cxx:101
#define R__ASSERT(e)
Definition: TError.h:96
#define N
static constexpr double bar
Basic string class.
Definition: TString.h:131
Double_t StruveH1(Double_t x)
Struve functions of order 0.
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
Int_t FloorNint(Double_t x)
Definition: TMath.h:706
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:348
double tgamma(double x)
The gamma function is defined to be the extension of the factorial to real numbers.
Double_t BesselJ1(Double_t x)
Bessel function J0(x) for any real x.
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
Double_t GamCf(Double_t a, Double_t x)
Computation of the incomplete gamma function P(a,x) via its continued fraction representation.
Definition: TMath.cxx:375
double cos(double)
Double_t Prob(Double_t chi2, Int_t ndf)
Double_t StudentI(Double_t T, Double_t ndf)
UInt_t Hash(ECaseCompare cmp=kExact) const
Return hash value.
Definition: TString.cxx:626
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
double beta(double x, double y)
Calculates the beta function.
double erfc(double x)
Complementary error function.
double sqrt(double)
Double_t BesselY1(Double_t x)
Bessel function Y0(x) for positive x.
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution: with where .
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
Double_t ASinH(Double_t)
Definition: TMath.cxx:64
Double_t Log10(Double_t x)
Definition: TMath.h:763
Double_t BesselI1(Double_t x)
modified Bessel function K_0(x)
static double p2(double t, double a, double b, double c)
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition: TMath.cxx:265
void Info(const char *location, const char *msgfmt,...)
Double_t VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t itype)
Double_t GamSer(Double_t a, Double_t x)
Computation of the incomplete gamma function P(a,x) via its series representation.
Definition: TMath.cxx:412
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 gamma_pdf(double x, double alpha, double theta, double x0=0)
Probability density function of the gamma distribution.
Double_t ATan2(Double_t, Double_t)
Definition: TMath.h:678
constexpr Double_t Pi()
Definition: TMath.h:38
double sin(double)
void BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
void Error(const char *location, const char *msgfmt,...)
Double_t Erfc(Double_t x)
Compute the complementary error function erfc(x).
Definition: TMath.cxx:194
static constexpr double second
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
Bool_t Even(Long_t a)
Definition: TMathBase.h:100
Double_t DiLog(Double_t x)
Modified Struve functions of order 1.
Definition: TMath.cxx:110
Double_t Double_v
Definition: Types.h:50
#define F(x, y, z)
double gamma(double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
ROOT::R::TRInterface & r
Definition: Object.C:4
Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2)
SVector< double, 2 > v
Definition: Dict.h:5
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
Double_t BesselY0(Double_t x)
Bessel function J1(x) for any real x.
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:184
auto * a
Definition: textangle.C:12
Double_t BetaDistI(Double_t x, Double_t p, Double_t q)
Double_t StruveL1(Double_t x)
Modified Struve functions of order 0.
Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r=4)
void VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt)
Bool_t Permute(Int_t n, Int_t *a)
Double_t StruveL0(Double_t x)
Struve functions of order 1.
double fdistribution_pdf(double x, double n, double m, double x0=0)
Probability density function of the F-distribution.
Double_t KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option)
#define NamespaceImp(name)
Definition: Rtypes.h:374
Double_t Binomial(Int_t n, Int_t k)
Double_t LandauI(Double_t x)
#define s1(x)
Definition: RSha256.hxx:91
Double_t ACos(Double_t)
Definition: TMath.h:667
static constexpr double m2
Double_t BesselJ0(Double_t x)
modified Bessel function K_1(x)
static double p1(double t, double a, double b)
Double_t ErfcInverse(Double_t x)
returns the inverse of the complementary error function x must be 0<x<2 implement using the quantile ...
Definition: TMath.cxx:237
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN) ...
Definition: TMath.h:907
ULong_t Hash(const void *txt, Int_t ntxt)
Double_t Poisson(Double_t x, Double_t par)
#define h(i)
Definition: RSha256.hxx:106
Double_t Student(Double_t T, Double_t ndf)
Double_t Beta(Double_t p, Double_t q)
Double_t Cos(Double_t)
Definition: TMath.h:640
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculate a gaussian function with mean and sigma.
Definition: TMath.cxx:448
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t BesselI0(Double_t x)
integer order modified Bessel function K_n(x)
Double_t BesselK(Int_t n, Double_t x)
integer order modified Bessel function I_n(x)
long Long_t
Definition: RtypesCore.h:50
#define d(i)
Definition: RSha256.hxx:102
static constexpr double pi2
Double_t Exp(Double_t x)
Definition: TMath.h:726
return c2
Definition: legend2.C:14
double Double_t
Definition: RtypesCore.h:55
Double_t StruveH0(Double_t x)
Bessel function Y1(x) for positive x.
int type
Definition: TGX11.cxx:120
unsigned long ULong_t
Definition: RtypesCore.h:51
Double_t y[n]
Definition: legend1.C:17
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
static constexpr double s
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:57
Double_t BesselK0(Double_t x)
modified Bessel function I_0(x)
Double_t BinomialI(Double_t p, Int_t n, Int_t k)
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Double_t BetaCf(Double_t x, Double_t a, Double_t b)
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Returns k_th order statistic of the array a of size n (k_th smallest element out of n elements)...
Definition: TMath.h:1453
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_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1)
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
Double_t Factorial(Int_t i)
Compute factorial(n).
Definition: TMath.cxx:247
Double_t Sin(Double_t)
Definition: TMath.h:636
TF1 * f1
Definition: legend1.C:11
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_t ASin(Double_t)
Definition: TMath.h:660
Double_t BesselK1(Double_t x)
modified Bessel function I_1(x)
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:486
#define c(i)
Definition: RSha256.hxx:101
Double_t LogNormal(Double_t x, Double_t sigma, Double_t theta=0, Double_t m=1)
Definition: first.py:1
Bool_t RootsCubic(const Double_t coef[4], Double_t &a, Double_t &b, Double_t &c)
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
double lgamma(double x)
Calculates the logarithm of the gamma function.
double exp(double)
Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b)
float * q
Definition: THbookFile.cxx:87
TMath.
Definition: TMathBase.h:34
Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1)
double inc_gamma(double a, double x)
Calculates the normalized (regularized) lower incomplete gamma function (lower integral) ...
static T Epsilon()
Returns minimum double representation.
Definition: TMath.h:934
const Bool_t kTRUE
Definition: RtypesCore.h:87
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition: TMath.h:712
return c3
Definition: legend3.C:15
const Int_t n
Definition: legend1.C:16
Double_t ATanH(Double_t)
Definition: TMath.cxx:90
constexpr Double_t PiOver2()
Definition: TMath.h:52
double log(double)
Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2)
static const double x3[11]
Double_t ATan(Double_t)
Definition: TMath.h:674