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