Logo ROOT   6.12/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 /// \author 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 /// \author 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 /// \author 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 /// \author 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 /// \author 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 ///
1080 /// \author Francois-Xavier Gentit
1081 
1083 {
1084  Bool_t complex = kFALSE;
1085  Double_t r,s,t,p,q,d,ps3,ps33,qs2,u,v,tmp,lnu,lnv,su,sv,y1,y2,y3;
1086  a = 0;
1087  b = 0;
1088  c = 0;
1089  if (coef[3] == 0) return complex;
1090  r = coef[2]/coef[3];
1091  s = coef[1]/coef[3];
1092  t = coef[0]/coef[3];
1093  p = s - (r*r)/3;
1094  ps3 = p/3;
1095  q = (2*r*r*r)/27.0 - (r*s)/3 + t;
1096  qs2 = q/2;
1097  ps33 = ps3*ps3*ps3;
1098  d = ps33 + qs2*qs2;
1099  if (d>=0) {
1100  complex = kTRUE;
1101  d = TMath::Sqrt(d);
1102  u = -qs2 + d;
1103  v = -qs2 - d;
1104  tmp = 1./3.;
1105  lnu = TMath::Log(TMath::Abs(u));
1106  lnv = TMath::Log(TMath::Abs(v));
1107  su = TMath::Sign(1.,u);
1108  sv = TMath::Sign(1.,v);
1109  u = su*TMath::Exp(tmp*lnu);
1110  v = sv*TMath::Exp(tmp*lnv);
1111  y1 = u + v;
1112  y2 = -y1/2;
1113  y3 = ((u-v)*TMath::Sqrt(3.))/2;
1114  tmp = r/3;
1115  a = y1 - tmp;
1116  b = y2 - tmp;
1117  c = y3;
1118  } else {
1119  Double_t phi,cphi,phis3,c1,c2,c3,pis3;
1120  ps3 = -ps3;
1121  ps33 = -ps33;
1122  cphi = -qs2/TMath::Sqrt(ps33);
1123  phi = TMath::ACos(cphi);
1124  phis3 = phi/3;
1125  pis3 = TMath::Pi()/3;
1126  c1 = TMath::Cos(phis3);
1127  c2 = TMath::Cos(pis3 + phis3);
1128  c3 = TMath::Cos(pis3 - phis3);
1129  tmp = TMath::Sqrt(ps3);
1130  y1 = 2*tmp*c1;
1131  y2 = -2*tmp*c2;
1132  y3 = -2*tmp*c3;
1133  tmp = r/3;
1134  a = y1 - tmp;
1135  b = y2 - tmp;
1136  c = y3 - tmp;
1137  }
1138  return complex;
1139 }
1140 
1141 ////////////////////////////////////////////////////////////////////////////////
1142 ///Computes sample quantiles, corresponding to the given probabilities
1143 ///Parameters:
1144 /// x -the data sample
1145 /// n - its size
1146 /// quantiles - computed quantiles are returned in there
1147 /// prob - probabilities where to compute quantiles
1148 /// nprob - size of prob array
1149 /// isSorted - is the input array x sorted?
1150 /// NOTE, that when the input is not sorted, an array of integers of size n needs
1151 /// to be allocated. It can be passed by the user in parameter index,
1152 /// or, if not passed, it will be allocated inside the function
1153 ///
1154 /// type - method to compute (from 1 to 9). Following types are provided:
1155 /// Discontinuous:
1156 /// type=1 - inverse of the empirical distribution function
1157 /// type=2 - like type 1, but with averaging at discontinuities
1158 /// type=3 - SAS definition: nearest even order statistic
1159 /// Piecwise linear continuous:
1160 /// In this case, sample quantiles can be obtained by linear interpolation
1161 /// between the k-th order statistic and p(k).
1162 /// type=4 - linear interpolation of empirical cdf, p(k)=k/n;
1163 /// type=5 - a very popular definition, p(k) = (k-0.5)/n;
1164 /// type=6 - used by Minitab and SPSS, p(k) = k/(n+1);
1165 /// type=7 - used by S-Plus and R, p(k) = (k-1)/(n-1);
1166 /// type=8 - resulting sample quantiles are approximately median unbiased
1167 /// regardless of the distribution of x. p(k) = (k-1/3)/(n+1/3);
1168 /// type=9 - resulting sample quantiles are approximately unbiased, when
1169 /// the sample comes from Normal distribution. p(k)=(k-3/8)/(n+1/4);
1170 ///
1171 /// default type = 7
1172 ///
1173 /// References:
1174 /// 1) Hyndman, R.J and Fan, Y, (1996) "Sample quantiles in statistical packages"
1175 /// American Statistician, 50, 361-365
1176 /// 2) R Project documentation for the function quantile of package {stats}
1177 
1178 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)
1179 {
1180 
1181  if (type<1 || type>9){
1182  printf("illegal value of type\n");
1183  return;
1184  }
1185  Int_t *ind = 0;
1186  Bool_t isAllocated = kFALSE;
1187  if (!isSorted){
1188  if (index) ind = index;
1189  else {
1190  ind = new Int_t[n];
1191  isAllocated = kTRUE;
1192  }
1193  }
1194 
1195  // re-implemented by L.M (9/11/2011) to fix bug https://savannah.cern.ch/bugs/?87251
1196  // following now exactly R implementation (available in library/stats/R/quantile.R )
1197  // which follows precisely Hyndman-Fan paper
1198  // (older implementation had a bug for type =3)
1199 
1200  for (Int_t i=0; i<nprob; i++){
1201 
1202  Double_t nppm = 0;
1203  Double_t gamma = 0;
1204  Int_t j = 0;
1205 
1206  //Discontinuous functions
1207  // type = 1,2, or 3
1208  if (type < 4 ){
1209  if (type == 3)
1210  nppm = n*prob[i]-0.5; // use m = -0.5
1211  else
1212  nppm = n*prob[i]; // use m = 0
1213 
1214  // be careful with machine precision
1215  double eps = 4 * TMath::Limits<Double_t>::Epsilon();
1216  j = TMath::FloorNint(nppm + eps);
1217 
1218  // LM : fix for numerical problems if nppm is actually equal to j, but results different for numerical error
1219  // g in the paper is nppm -j
1220  if (type == 1)
1221  gamma = ( (nppm -j) > j*TMath::Limits<Double_t>::Epsilon() ) ? 1 : 0;
1222  else if (type == 2)
1223  gamma = ( (nppm -j) > j*TMath::Limits<Double_t>::Epsilon() ) ? 1 : 0.5;
1224  else if (type == 3)
1225  // gamma = 0 for g=0 and j even
1226  gamma = ( TMath::Abs(nppm -j) <= j*TMath::Limits<Double_t>::Epsilon() && TMath::Even(j) ) ? 0 : 1;
1227 
1228  }
1229  else {
1230  // for continuous quantiles type 4 to 9)
1231  // define alpha and beta
1232  double a = 0;
1233  double b = 0;
1234  if (type == 4) { a = 0; b = 1; }
1235  else if (type == 5) { a = 0.5; b = 0.5; }
1236  else if (type == 6) { a = 0.; b = 0.; }
1237  else if (type == 7) { a = 1.; b = 1.; }
1238  else if (type == 8) { a = 1./3.; b = a; }
1239  else if (type == 9) { a = 3./8.; b = a; }
1240 
1241  // be careful with machine precision
1242  double eps = 4 * TMath::Limits<Double_t>::Epsilon();
1243  nppm = a + prob[i] * ( n + 1 -a -b); // n * p + m
1244  j = TMath::FloorNint(nppm + eps);
1245  gamma = nppm - j;
1246  if (gamma < eps) gamma = 0;
1247  }
1248 
1249  // since index j starts from 1 first is j-1 and second is j
1250  // add protection to keep indices in range [0,n-1]
1251  int first = (j > 0 && j <=n) ? j-1 : ( j <=0 ) ? 0 : n-1;
1252  int second = (j > 0 && j < n) ? j : ( j <=0 ) ? 0 : n-1;
1253 
1254  Double_t xj, xjj;
1255  if (isSorted){
1256  xj = x[first];
1257  xjj = x[second];
1258  } else {
1259  xj = TMath::KOrdStat(n, x, first, ind);
1260  xjj = TMath::KOrdStat(n, x, second, ind);
1261  }
1262  quantiles[i] = (1-gamma)*xj + gamma*xjj;
1263  // printf("x[j] = %12f x[j+1] = %12f \n",xj, xjj);
1264 
1265  }
1266 
1267 
1268 
1269  if (isAllocated)
1270  delete [] ind;
1271 }
1272 
1273 ////////////////////////////////////////////////////////////////////////////////
1274 /// Bubble sort variant to obtain the order of an array's elements into
1275 /// an index in order to do more useful things than the standard built
1276 /// in functions.
1277 /// *arr1 is unchanged;
1278 /// *arr2 is the array of indicies corresponding to the descending value
1279 /// of arr1 with arr2[0] corresponding to the largest arr1 value and
1280 /// arr2[Narr] the smallest.
1281 ///
1282 /// \author Adrian Bevan (bevan@slac.stanford.edu)
1283 
1284 void TMath::BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
1285 {
1286  if (Narr <= 0) return;
1287  double *localArr1 = new double[Narr];
1288  int *localArr2 = new int[Narr];
1289  int iEl;
1290  int iEl2;
1291 
1292  for(iEl = 0; iEl < Narr; iEl++) {
1293  localArr1[iEl] = arr1[iEl];
1294  localArr2[iEl] = iEl;
1295  }
1296 
1297  for (iEl = 0; iEl < Narr; iEl++) {
1298  for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1299  if (localArr1[iEl2-1] < localArr1[iEl2]) {
1300  double tmp = localArr1[iEl2-1];
1301  localArr1[iEl2-1] = localArr1[iEl2];
1302  localArr1[iEl2] = tmp;
1303 
1304  int tmp2 = localArr2[iEl2-1];
1305  localArr2[iEl2-1] = localArr2[iEl2];
1306  localArr2[iEl2] = tmp2;
1307  }
1308  }
1309  }
1310 
1311  for (iEl = 0; iEl < Narr; iEl++) {
1312  arr2[iEl] = localArr2[iEl];
1313  }
1314  delete [] localArr2;
1315  delete [] localArr1;
1316 }
1317 
1318 ////////////////////////////////////////////////////////////////////////////////
1319 /// Opposite ordering of the array arr2[] to that of BubbleHigh.
1320 ///
1321 /// \author Adrian Bevan (bevan@slac.stanford.edu)
1322 
1323 void TMath::BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
1324 {
1325  if (Narr <= 0) return;
1326  double *localArr1 = new double[Narr];
1327  int *localArr2 = new int[Narr];
1328  int iEl;
1329  int iEl2;
1330 
1331  for (iEl = 0; iEl < Narr; iEl++) {
1332  localArr1[iEl] = arr1[iEl];
1333  localArr2[iEl] = iEl;
1334  }
1335 
1336  for (iEl = 0; iEl < Narr; iEl++) {
1337  for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1338  if (localArr1[iEl2-1] > localArr1[iEl2]) {
1339  double tmp = localArr1[iEl2-1];
1340  localArr1[iEl2-1] = localArr1[iEl2];
1341  localArr1[iEl2] = tmp;
1342 
1343  int tmp2 = localArr2[iEl2-1];
1344  localArr2[iEl2-1] = localArr2[iEl2];
1345  localArr2[iEl2] = tmp2;
1346  }
1347  }
1348  }
1349 
1350  for (iEl = 0; iEl < Narr; iEl++) {
1351  arr2[iEl] = localArr2[iEl];
1352  }
1353  delete [] localArr2;
1354  delete [] localArr1;
1355 }
1356 
1357 
1358 ////////////////////////////////////////////////////////////////////////////////
1359 /// Calculates hash index from any char string.
1360 /// Based on pre-calculated table of 256 specially selected numbers.
1361 /// These numbers are selected in such a way, that for string
1362 /// length == 4 (integer number) the hash is unambiguous, i.e.
1363 /// from hash value we can recalculate input (no degeneration).
1364 ///
1365 /// The quality of hash method is good enough, that
1366 /// "random" numbers made as R = Hash(1), Hash(2), ...Hash(N)
1367 /// tested by <R>, <R*R>, <Ri*Ri+1> gives the same result
1368 /// as for libc rand().
1369 ///
1370 /// For string: i = TMath::Hash(string,nstring);
1371 /// For int: i = TMath::Hash(&intword,sizeof(int));
1372 /// For pointer: i = TMath::Hash(&pointer,sizeof(void*));
1373 ///
1374 /// V.Perev
1375 /// This function is kept for back compatibility. The code previously in this function
1376 /// has been moved to the static function TString::Hash
1377 
1378 ULong_t TMath::Hash(const void *txt, Int_t ntxt)
1379 {
1380  return TString::Hash(txt,ntxt);
1381 }
1382 
1383 
1384 ////////////////////////////////////////////////////////////////////////////////
1385 
1386 ULong_t TMath::Hash(const char *txt)
1387 {
1388  return Hash(txt, Int_t(strlen(txt)));
1389 }
1390 
1391 ////////////////////////////////////////////////////////////////////////////////
1392 /// Compute the modified Bessel function I_0(x) for any real x.
1393 ///
1394 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1395 
1397 {
1398  // Parameters of the polynomial approximation
1399  const Double_t p1=1.0, p2=3.5156229, p3=3.0899424,
1400  p4=1.2067492, p5=0.2659732, p6=3.60768e-2, p7=4.5813e-3;
1401 
1402  const Double_t q1= 0.39894228, q2= 1.328592e-2, q3= 2.25319e-3,
1403  q4=-1.57565e-3, q5= 9.16281e-3, q6=-2.057706e-2,
1404  q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3;
1405 
1406  const Double_t k1 = 3.75;
1407  Double_t ax = TMath::Abs(x);
1408 
1409  Double_t y=0, result=0;
1410 
1411  if (ax < k1) {
1412  Double_t xx = x/k1;
1413  y = xx*xx;
1414  result = p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))));
1415  } else {
1416  y = k1/ax;
1417  result = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
1418  }
1419  return result;
1420 }
1421 
1422 ////////////////////////////////////////////////////////////////////////////////
1423 /// Compute the modified Bessel function K_0(x) for positive real x.
1424 ///
1425 /// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1426 /// Applied Mathematics Series vol. 55 (1964), Washington.
1427 ///
1428 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1429 
1431 {
1432  // Parameters of the polynomial approximation
1433  const Double_t p1=-0.57721566, p2=0.42278420, p3=0.23069756,
1434  p4= 3.488590e-2, p5=2.62698e-3, p6=1.0750e-4, p7=7.4e-6;
1435 
1436  const Double_t q1= 1.25331414, q2=-7.832358e-2, q3= 2.189568e-2,
1437  q4=-1.062446e-2, q5= 5.87872e-3, q6=-2.51540e-3, q7=5.3208e-4;
1438 
1439  if (x <= 0) {
1440  Error("TMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
1441  return 0;
1442  }
1443 
1444  Double_t y=0, result=0;
1445 
1446  if (x <= 2) {
1447  y = x*x/4;
1448  result = (-log(x/2.)*TMath::BesselI0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1449  } else {
1450  y = 2/x;
1451  result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
1452  }
1453  return result;
1454 }
1455 
1456 ////////////////////////////////////////////////////////////////////////////////
1457 /// Compute the modified Bessel function I_1(x) for any real x.
1458 ///
1459 /// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1460 /// Applied Mathematics Series vol. 55 (1964), Washington.
1461 ///
1462 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1463 
1465 {
1466  // Parameters of the polynomial approximation
1467  const Double_t p1=0.5, p2=0.87890594, p3=0.51498869,
1468  p4=0.15084934, p5=2.658733e-2, p6=3.01532e-3, p7=3.2411e-4;
1469 
1470  const Double_t q1= 0.39894228, q2=-3.988024e-2, q3=-3.62018e-3,
1471  q4= 1.63801e-3, q5=-1.031555e-2, q6= 2.282967e-2,
1472  q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3;
1473 
1474  const Double_t k1 = 3.75;
1475  Double_t ax = TMath::Abs(x);
1476 
1477  Double_t y=0, result=0;
1478 
1479  if (ax < k1) {
1480  Double_t xx = x/k1;
1481  y = xx*xx;
1482  result = x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1483  } else {
1484  y = k1/ax;
1485  result = (exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
1486  if (x < 0) result = -result;
1487  }
1488  return result;
1489 }
1490 
1491 ////////////////////////////////////////////////////////////////////////////////
1492 /// Compute the modified Bessel function K_1(x) for positive real x.
1493 ///
1494 /// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1495 /// Applied Mathematics Series vol. 55 (1964), Washington.
1496 ///
1497 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1498 
1500 {
1501  // Parameters of the polynomial approximation
1502  const Double_t p1= 1., p2= 0.15443144, p3=-0.67278579,
1503  p4=-0.18156897, p5=-1.919402e-2, p6=-1.10404e-3, p7=-4.686e-5;
1504 
1505  const Double_t q1= 1.25331414, q2= 0.23498619, q3=-3.655620e-2,
1506  q4= 1.504268e-2, q5=-7.80353e-3, q6= 3.25614e-3, q7=-6.8245e-4;
1507 
1508  if (x <= 0) {
1509  Error("TMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
1510  return 0;
1511  }
1512 
1513  Double_t y=0,result=0;
1514 
1515  if (x <= 2) {
1516  y = x*x/4;
1517  result = (log(x/2.)*TMath::BesselI1(x))+(1./x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1518  } else {
1519  y = 2/x;
1520  result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
1521  }
1522  return result;
1523 }
1524 
1525 ////////////////////////////////////////////////////////////////////////////////
1526 /// Compute the Integer Order Modified Bessel function K_n(x)
1527 /// for n=0,1,2,... and positive real x.
1528 ///
1529 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1530 
1532 {
1533  if (x <= 0 || n < 0) {
1534  Error("TMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
1535  return 0;
1536  }
1537 
1538  if (n==0) return TMath::BesselK0(x);
1539  if (n==1) return TMath::BesselK1(x);
1540 
1541  // Perform upward recurrence for all x
1542  Double_t tox = 2/x;
1543  Double_t bkm = TMath::BesselK0(x);
1544  Double_t bk = TMath::BesselK1(x);
1545  Double_t bkp = 0;
1546  for (Int_t j=1; j<n; j++) {
1547  bkp = bkm+Double_t(j)*tox*bk;
1548  bkm = bk;
1549  bk = bkp;
1550  }
1551  return bk;
1552 }
1553 
1554 ////////////////////////////////////////////////////////////////////////////////
1555 /// Compute the Integer Order Modified Bessel function I_n(x)
1556 /// for n=0,1,2,... and any real x.
1557 ///
1558 /// \author NvE 12-mar-2000 UU-SAP Utrecht
1559 
1561 {
1562  Int_t iacc = 40; // Increase to enhance accuracy
1563  const Double_t kBigPositive = 1.e10;
1564  const Double_t kBigNegative = 1.e-10;
1565 
1566  if (n < 0) {
1567  Error("TMath::BesselI", "*I* Invalid argument (n,x) = (%d, %g)\n",n,x);
1568  return 0;
1569  }
1570 
1571  if (n==0) return TMath::BesselI0(x);
1572  if (n==1) return TMath::BesselI1(x);
1573 
1574  if (x == 0) return 0;
1575  if (TMath::Abs(x) > kBigPositive) return 0;
1576 
1577  Double_t tox = 2/TMath::Abs(x);
1578  Double_t bip = 0, bim = 0;
1579  Double_t bi = 1;
1580  Double_t result = 0;
1581  Int_t m = 2*((n+Int_t(sqrt(Float_t(iacc*n)))));
1582  for (Int_t j=m; j>=1; j--) {
1583  bim = bip+Double_t(j)*tox*bi;
1584  bip = bi;
1585  bi = bim;
1586  // Renormalise to prevent overflows
1587  if (TMath::Abs(bi) > kBigPositive) {
1588  result *= kBigNegative;
1589  bi *= kBigNegative;
1590  bip *= kBigNegative;
1591  }
1592  if (j==n) result=bip;
1593  }
1594 
1595  result *= TMath::BesselI0(x)/bi; // Normalise with BesselI0(x)
1596  if ((x < 0) && (n%2 == 1)) result = -result;
1597 
1598  return result;
1599 }
1600 
1601 ////////////////////////////////////////////////////////////////////////////////
1602 /// Returns the Bessel function J0(x) for any real x.
1603 
1605 {
1606  Double_t ax,z;
1607  Double_t xx,y,result,result1,result2;
1608  const Double_t p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
1609  const Double_t p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
1610  const Double_t p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
1611  const Double_t p10 = 59272.64853, p11 = 267.8532712;
1612 
1613  const Double_t q1 = 0.785398164;
1614  const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1615  const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1616  const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1617  const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1618  const Double_t q10 = 0.934935152e-7, q11 = 0.636619772;
1619 
1620  if ((ax=fabs(x)) < 8) {
1621  y=x*x;
1622  result1 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6))));
1623  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1624  result = result1/result2;
1625  } else {
1626  z = 8/ax;
1627  y = z*z;
1628  xx = ax-q1;
1629  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1630  result2 = q6 + y*(q7 + y*(q8 + y*(q9 - y*q10)));
1631  result = sqrt(q11/ax)*(cos(xx)*result1-z*sin(xx)*result2);
1632  }
1633  return result;
1634 }
1635 
1636 ////////////////////////////////////////////////////////////////////////////////
1637 /// Returns the Bessel function J1(x) for any real x.
1638 
1640 {
1641  Double_t ax,z;
1642  Double_t xx,y,result,result1,result2;
1643  const Double_t p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
1644  const Double_t p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
1645  const Double_t p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
1646  const Double_t p10 = 99447.43394, p11 = 376.9991397;
1647 
1648  const Double_t q1 = 2.356194491;
1649  const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1650  const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1651  const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1652  const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1653  const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1654 
1655  if ((ax=fabs(x)) < 8) {
1656  y=x*x;
1657  result1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6)))));
1658  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1659  result = result1/result2;
1660  } else {
1661  z = 8/ax;
1662  y = z*z;
1663  xx = ax-q1;
1664  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1665  result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1666  result = sqrt(q11/ax)*(cos(xx)*result1-z*sin(xx)*result2);
1667  if (x < 0) result = -result;
1668  }
1669  return result;
1670 }
1671 
1672 ////////////////////////////////////////////////////////////////////////////////
1673 /// Returns the Bessel function Y0(x) for positive x.
1674 
1676 {
1677  Double_t z,xx,y,result,result1,result2;
1678  const Double_t p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
1679  const Double_t p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
1680  const Double_t p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
1681  const Double_t p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
1682 
1683  const Double_t q1 = 0.785398164;
1684  const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1685  const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1686  const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1687  const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1688  const Double_t q10 = -0.934945152e-7, q11 = 0.636619772;
1689 
1690  if (x < 8) {
1691  y = x*x;
1692  result1 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6))));
1693  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1694  result = (result1/result2) + p12*TMath::BesselJ0(x)*log(x);
1695  } else {
1696  z = 8/x;
1697  y = z*z;
1698  xx = x-q1;
1699  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1700  result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1701  result = sqrt(q11/x)*(sin(xx)*result1+z*cos(xx)*result2);
1702  }
1703  return result;
1704 }
1705 
1706 ////////////////////////////////////////////////////////////////////////////////
1707 /// Returns the Bessel function Y1(x) for positive x.
1708 
1710 {
1711  Double_t z,xx,y,result,result1,result2;
1712  const Double_t p1 = -0.4900604943e13, p2 = 0.1275274390e13;
1713  const Double_t p3 = -0.5153438139e11, p4 = 0.7349264551e9;
1714  const Double_t p5 = -0.4237922726e7, p6 = 0.8511937935e4;
1715  const Double_t p7 = 0.2499580570e14, p8 = 0.4244419664e12;
1716  const Double_t p9 = 0.3733650367e10, p10 = 0.2245904002e8;
1717  const Double_t p11 = 0.1020426050e6, p12 = 0.3549632885e3;
1718  const Double_t p13 = 0.636619772;
1719  const Double_t q1 = 2.356194491;
1720  const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1721  const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1722  const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1723  const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1724  const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1725 
1726  if (x < 8) {
1727  y=x*x;
1728  result1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6)))));
1729  result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y*(p12+y)))));
1730  result = (result1/result2) + p13*(TMath::BesselJ1(x)*log(x)-1/x);
1731  } else {
1732  z = 8/x;
1733  y = z*z;
1734  xx = x-q1;
1735  result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1736  result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1737  result = sqrt(q11/x)*(sin(xx)*result1+z*cos(xx)*result2);
1738  }
1739  return result;
1740 }
1741 
1742 ////////////////////////////////////////////////////////////////////////////////
1743 /// Struve Functions of Order 0
1744 ///
1745 /// Converted from CERNLIB M342 by Rene Brun.
1746 
1748 {
1749  const Int_t n1 = 15;
1750  const Int_t n2 = 25;
1751  const Double_t c1[16] = { 1.00215845609911981, -1.63969292681309147,
1752  1.50236939618292819, -.72485115302121872,
1753  .18955327371093136, -.03067052022988,
1754  .00337561447375194, -2.6965014312602e-4,
1755  1.637461692612e-5, -7.8244408508e-7,
1756  3.021593188e-8, -9.6326645e-10,
1757  2.579337e-11, -5.8854e-13,
1758  1.158e-14, -2e-16 };
1759  const Double_t c2[26] = { .99283727576423943, -.00696891281138625,
1760  1.8205103787037e-4, -1.063258252844e-5,
1761  9.8198294287e-7, -1.2250645445e-7,
1762  1.894083312e-8, -3.44358226e-9,
1763  7.1119102e-10, -1.6288744e-10,
1764  4.065681e-11, -1.091505e-11,
1765  3.12005e-12, -9.4202e-13,
1766  2.9848e-13, -9.872e-14,
1767  3.394e-14, -1.208e-14,
1768  4.44e-15, -1.68e-15,
1769  6.5e-16, -2.6e-16,
1770  1.1e-16, -4e-17,
1771  2e-17, -1e-17 };
1772 
1773  const Double_t c0 = 2/TMath::Pi();
1774 
1775  Int_t i;
1776  Double_t alfa, h, r, y, b0, b1, b2;
1777  Double_t v = TMath::Abs(x);
1778 
1779  v = TMath::Abs(x);
1780  if (v < 8) {
1781  y = v/8;
1782  h = 2*y*y -1;
1783  alfa = h + h;
1784  b0 = 0;
1785  b1 = 0;
1786  b2 = 0;
1787  for (i = n1; i >= 0; --i) {
1788  b0 = c1[i] + alfa*b1 - b2;
1789  b2 = b1;
1790  b1 = b0;
1791  }
1792  h = y*(b0 - h*b2);
1793  } else {
1794  r = 1/v;
1795  h = 128*r*r -1;
1796  alfa = h + h;
1797  b0 = 0;
1798  b1 = 0;
1799  b2 = 0;
1800  for (i = n2; i >= 0; --i) {
1801  b0 = c2[i] + alfa*b1 - b2;
1802  b2 = b1;
1803  b1 = b0;
1804  }
1805  h = TMath::BesselY0(v) + r*c0*(b0 - h*b2);
1806  }
1807  if (x < 0) h = -h;
1808  return h;
1809 }
1810 
1811 ////////////////////////////////////////////////////////////////////////////////
1812 /// Struve Functions of Order 1
1813 ///
1814 /// Converted from CERNLIB M342 by Rene Brun.
1815 
1817 {
1818  const Int_t n3 = 16;
1819  const Int_t n4 = 22;
1820  const Double_t c3[17] = { .5578891446481605, -.11188325726569816,
1821  -.16337958125200939, .32256932072405902,
1822  -.14581632367244242, .03292677399374035,
1823  -.00460372142093573, 4.434706163314e-4,
1824  -3.142099529341e-5, 1.7123719938e-6,
1825  -7.416987005e-8, 2.61837671e-9,
1826  -7.685839e-11, 1.9067e-12,
1827  -4.052e-14, 7.5e-16,
1828  -1e-17 };
1829  const Double_t c4[23] = { 1.00757647293865641, .00750316051248257,
1830  -7.043933264519e-5, 2.66205393382e-6,
1831  -1.8841157753e-7, 1.949014958e-8,
1832  -2.6126199e-9, 4.236269e-10,
1833  -7.955156e-11, 1.679973e-11,
1834  -3.9072e-12, 9.8543e-13,
1835  -2.6636e-13, 7.645e-14,
1836  -2.313e-14, 7.33e-15,
1837  -2.42e-15, 8.3e-16,
1838  -3e-16, 1.1e-16,
1839  -4e-17, 2e-17,-1e-17 };
1840 
1841  const Double_t c0 = 2/TMath::Pi();
1842  const Double_t cc = 2/(3*TMath::Pi());
1843 
1844  Int_t i, i1;
1845  Double_t alfa, h, r, y, b0, b1, b2;
1846  Double_t v = TMath::Abs(x);
1847 
1848  if (v == 0) {
1849  h = 0;
1850  } else if (v <= 0.3) {
1851  y = v*v;
1852  r = 1;
1853  h = 1;
1854  i1 = (Int_t)(-8. / TMath::Log10(v));
1855  for (i = 1; i <= i1; ++i) {
1856  h = -h*y / ((2*i+ 1)*(2*i + 3));
1857  r += h;
1858  }
1859  h = cc*y*r;
1860  } else if (v < 8) {
1861  h = v*v/32 -1;
1862  alfa = h + h;
1863  b0 = 0;
1864  b1 = 0;
1865  b2 = 0;
1866  for (i = n3; i >= 0; --i) {
1867  b0 = c3[i] + alfa*b1 - b2;
1868  b2 = b1;
1869  b1 = b0;
1870  }
1871  h = b0 - h*b2;
1872  } else {
1873  h = 128/(v*v) -1;
1874  alfa = h + h;
1875  b0 = 0;
1876  b1 = 0;
1877  b2 = 0;
1878  for (i = n4; i >= 0; --i) {
1879  b0 = c4[i] + alfa*b1 - b2;
1880  b2 = b1;
1881  b1 = b0;
1882  }
1883  h = TMath::BesselY1(v) + c0*(b0 - h*b2);
1884  }
1885  return h;
1886 }
1887 
1888 
1889 ////////////////////////////////////////////////////////////////////////////////
1890 /// Modified Struve Function of Order 0.
1891 /// By Kirill Filimonov.
1892 
1894 {
1895  const Double_t pi=TMath::Pi();
1896 
1897  Double_t s=1.0;
1898  Double_t r=1.0;
1899 
1900  Double_t a0,sl0,a1,bi0;
1901 
1902  Int_t km,i;
1903 
1904  if (x<=20.) {
1905  a0=2.0*x/pi;
1906  for (i=1; i<=60;i++) {
1907  r*=(x/(2*i+1))*(x/(2*i+1));
1908  s+=r;
1909  if(TMath::Abs(r/s)<1.e-12) break;
1910  }
1911  sl0=a0*s;
1912  } else {
1913  km=int(5*(x+1.0));
1914  if(x>=50.0)km=25;
1915  for (i=1; i<=km; i++) {
1916  r*=(2*i-1)*(2*i-1)/x/x;
1917  s+=r;
1918  if(TMath::Abs(r/s)<1.0e-12) break;
1919  }
1920  a1=TMath::Exp(x)/TMath::Sqrt(2*pi*x);
1921  r=1.0;
1922  bi0=1.0;
1923  for (i=1; i<=16; i++) {
1924  r=0.125*r*(2.0*i-1.0)*(2.0*i-1.0)/(i*x);
1925  bi0+=r;
1926  if(TMath::Abs(r/bi0)<1.0e-12) break;
1927  }
1928 
1929  bi0=a1*bi0;
1930  sl0=-2.0/(pi*x)*s+bi0;
1931  }
1932  return sl0;
1933 }
1934 
1935 ////////////////////////////////////////////////////////////////////////////////
1936 /// Modified Struve Function of Order 1.
1937 /// By Kirill Filimonov.
1938 
1940 {
1941  const Double_t pi=TMath::Pi();
1942  Double_t a1,sl1,bi1,s;
1943  Double_t r=1.0;
1944  Int_t km,i;
1945 
1946  if (x<=20.) {
1947  s=0.0;
1948  for (i=1; i<=60;i++) {
1949  r*=x*x/(4.0*i*i-1.0);
1950  s+=r;
1951  if(TMath::Abs(r)<TMath::Abs(s)*1.e-12) break;
1952  }
1953  sl1=2.0/pi*s;
1954  } else {
1955  s=1.0;
1956  km=int(0.5*x);
1957  if(x>50.0)km=25;
1958  for (i=1; i<=km; i++) {
1959  r*=(2*i+3)*(2*i+1)/x/x;
1960  s+=r;
1961  if(TMath::Abs(r/s)<1.0e-12) break;
1962  }
1963  sl1=2.0/pi*(-1.0+1.0/(x*x)+3.0*s/(x*x*x*x));
1964  a1=TMath::Exp(x)/TMath::Sqrt(2*pi*x);
1965  r=1.0;
1966  bi1=1.0;
1967  for (i=1; i<=16; i++) {
1968  r=-0.125*r*(4.0-(2.0*i-1.0)*(2.0*i-1.0))/(i*x);
1969  bi1+=r;
1970  if(TMath::Abs(r/bi1)<1.0e-12) break;
1971  }
1972  sl1+=a1*bi1;
1973  }
1974  return sl1;
1975 }
1976 
1977 ////////////////////////////////////////////////////////////////////////////////
1978 /// Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
1979 
1981 {
1983 }
1984 
1985 ////////////////////////////////////////////////////////////////////////////////
1986 /// Continued fraction evaluation by modified Lentz's method
1987 /// used in calculation of incomplete Beta function.
1988 
1990 {
1991  Int_t itmax = 500;
1992  Double_t eps = 3.e-14;
1993  Double_t fpmin = 1.e-30;
1994 
1995  Int_t m, m2;
1996  Double_t aa, c, d, del, qab, qam, qap;
1997  Double_t h;
1998  qab = a+b;
1999  qap = a+1.0;
2000  qam = a-1.0;
2001  c = 1.0;
2002  d = 1.0 - qab*x/qap;
2003  if (TMath::Abs(d)<fpmin) d=fpmin;
2004  d=1.0/d;
2005  h=d;
2006  for (m=1; m<=itmax; m++) {
2007  m2=m*2;
2008  aa = m*(b-m)*x/((qam+ m2)*(a+m2));
2009  d = 1.0 +aa*d;
2010  if(TMath::Abs(d)<fpmin) d = fpmin;
2011  c = 1 +aa/c;
2012  if (TMath::Abs(c)<fpmin) c = fpmin;
2013  d=1.0/d;
2014  h*=d*c;
2015  aa = -(a+m)*(qab +m)*x/((a+m2)*(qap+m2));
2016  d=1.0+aa*d;
2017  if(TMath::Abs(d)<fpmin) d = fpmin;
2018  c = 1.0 +aa/c;
2019  if (TMath::Abs(c)<fpmin) c = fpmin;
2020  d=1.0/d;
2021  del = d*c;
2022  h*=del;
2023  if (TMath::Abs(del-1)<=eps) break;
2024  }
2025  if (m>itmax) {
2026  Info("TMath::BetaCf", "a or b too big, or itmax too small, a=%g, b=%g, x=%g, h=%g, itmax=%d",
2027  a,b,x,h,itmax);
2028  }
2029  return h;
2030 }
2031 
2032 ////////////////////////////////////////////////////////////////////////////////
2033 /// Computes the probability density function of the Beta distribution
2034 /// (the distribution function is computed in BetaDistI).
2035 /// The first argument is the point, where the function will be
2036 /// computed, second and third are the function parameters.
2037 /// Since the Beta distribution is bounded on both sides, it's often
2038 /// used to represent processes with natural lower and upper limits.
2039 
2041 {
2042  if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2043  Error("TMath::BetaDist", "parameter value outside allowed range");
2044  return 0;
2045  }
2046  Double_t beta = TMath::Beta(p, q);
2047  Double_t r = TMath::Power(x, p-1)*TMath::Power(1-x, q-1)/beta;
2048  return r;
2049 }
2050 
2051 ////////////////////////////////////////////////////////////////////////////////
2052 /// Computes the distribution function of the Beta distribution.
2053 /// The first argument is the point, where the function will be
2054 /// computed, second and third are the function parameters.
2055 /// Since the Beta distribution is bounded on both sides, it's often
2056 /// used to represent processes with natural lower and upper limits.
2057 
2059 {
2060  if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2061  Error("TMath::BetaDistI", "parameter value outside allowed range");
2062  return 0;
2063  }
2064  Double_t betai = TMath::BetaIncomplete(x, p, q);
2065  return betai;
2066 }
2067 
2068 ////////////////////////////////////////////////////////////////////////////////
2069 /// Calculates the incomplete Beta-function.
2070 
2072 {
2074 }
2075 
2076 ////////////////////////////////////////////////////////////////////////////////
2077 /// Calculate the binomial coefficient n over k.
2078 
2080 {
2081  if (n<0 || k<0 || n<k) return TMath::SignalingNaN();
2082  if (k==0 || n==k) return 1;
2083 
2084  Int_t k1=TMath::Min(k,n-k);
2085  Int_t k2=n-k1;
2086  Double_t fact=k2+1;
2087  for (Double_t i=k1;i>1.;--i)
2088  fact *= (k2+i)/i;
2089  return fact;
2090 }
2091 
2092 ////////////////////////////////////////////////////////////////////////////////
2093 /// Suppose an event occurs with probability _p_ per trial
2094 /// Then the probability P of its occurring _k_ or more times
2095 /// in _n_ trials is termed a cumulative binomial probability
2096 /// the formula is P = sum_from_j=k_to_n(TMath::Binomial(n, j)*
2097 /// *TMath::Power(p, j)*TMath::Power(1-p, n-j)
2098 /// For _n_ larger than 12 BetaIncomplete is a much better way
2099 /// to evaluate the sum than would be the straightforward sum calculation
2100 /// for _n_ smaller than 12 either method is acceptable
2101 /// ("Numerical Recipes")
2102 /// --implementation by Anna Kreshuk
2103 
2105 {
2106  if(k <= 0) return 1.0;
2107  if(k > n) return 0.0;
2108  if(k == n) return TMath::Power(p, n);
2109 
2110  return BetaIncomplete(p, Double_t(k), Double_t(n-k+1));
2111 }
2112 
2113 ////////////////////////////////////////////////////////////////////////////////
2114 /// Computes the density of Cauchy distribution at point x
2115 /// by default, standard Cauchy distribution is used (t=0, s=1)
2116 /// t is the location parameter
2117 /// s is the scale parameter
2118 /// The Cauchy distribution, also called Lorentzian distribution,
2119 /// is a continuous distribution describing resonance behavior
2120 /// The mean and standard deviation of the Cauchy distribution are undefined.
2121 /// The practical meaning of this is that collecting 1,000 data points gives
2122 /// no more accurate an estimate of the mean and standard deviation than
2123 /// does a single point.
2124 /// The formula was taken from "Engineering Statistics Handbook" on site
2125 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3663.htm
2126 /// Implementation by Anna Kreshuk.
2127 /// Example:
2128 /// TF1* fc = new TF1("fc", "TMath::CauchyDist(x, [0], [1])", -5, 5);
2129 /// fc->SetParameters(0, 1);
2130 /// fc->Draw();
2131 
2133 {
2134  Double_t temp= (x-t)*(x-t)/(s*s);
2135  Double_t result = 1/(s*TMath::Pi()*(1+temp));
2136  return result;
2137 }
2138 
2139 ////////////////////////////////////////////////////////////////////////////////
2140 /// Evaluate the quantiles of the chi-squared probability distribution function.
2141 /// Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
2142 /// implemented by Anna Kreshuk.
2143 /// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
2144 /// Parameters:
2145 /// p - the probability value, at which the quantile is computed
2146 /// ndf - number of degrees of freedom
2147 
2149 {
2150  Double_t c[]={0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
2151  4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
2152  84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
2153  210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
2154  462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
2155  932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
2156  2520.0, 5040.0};
2157  Double_t e = 5e-7;
2158  Double_t aa = 0.6931471806;
2159  Int_t maxit = 20;
2160  Double_t ch, p1, p2, q, t, a, b, x;
2161  Double_t s1, s2, s3, s4, s5, s6;
2162 
2163  if (ndf <= 0) return 0;
2164 
2165  Double_t g = TMath::LnGamma(0.5*ndf);
2166 
2167  Double_t xx = 0.5 * ndf;
2168  Double_t cp = xx - 1;
2169  if (ndf >= TMath::Log(p)*(-c[5])){
2170  //starting approximation for ndf less than or equal to 0.32
2171  if (ndf > c[3]) {
2172  x = TMath::NormQuantile(p);
2173  //starting approximation using Wilson and Hilferty estimate
2174  p1=c[2]/ndf;
2175  ch = ndf*TMath::Power((x*TMath::Sqrt(p1) + 1 - p1), 3);
2176  if (ch > c[6]*ndf + 6)
2177  ch = -2 * (TMath::Log(1-p) - cp * TMath::Log(0.5 * ch) + g);
2178  } else {
2179  ch = c[4];
2180  a = TMath::Log(1-p);
2181  do{
2182  q = ch;
2183  p1 = 1 + ch * (c[7]+ch);
2184  p2 = ch * (c[9] + ch * (c[8] + ch));
2185  t = -0.5 + (c[7] + 2 * ch) / p1 - (c[9] + ch * (c[10] + 3 * ch)) / p2;
2186  ch = ch - (1 - TMath::Exp(a + g + 0.5 * ch + cp * aa) *p2 / p1) / t;
2187  }while (TMath::Abs(q/ch - 1) > c[1]);
2188  }
2189  } else {
2190  ch = TMath::Power((p * xx * TMath::Exp(g + xx * aa)),(1./xx));
2191  if (ch < e) return ch;
2192  }
2193 //call to algorithm AS 239 and calculation of seven term Taylor series
2194  for (Int_t i=0; i<maxit; i++){
2195  q = ch;
2196  p1 = 0.5 * ch;
2197  p2 = p - TMath::Gamma(xx, p1);
2198 
2199  t = p2 * TMath::Exp(xx * aa + g + p1 - cp * TMath::Log(ch));
2200  b = t / ch;
2201  a = 0.5 * t - b * cp;
2202  s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] +c[11] * a))))) / c[24];
2203  s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
2204  s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
2205  s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
2206  s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
2207  s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
2208  ch = ch + t * (1 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
2209  if (TMath::Abs(q / ch - 1) > e) break;
2210  }
2211  return ch;
2212 }
2213 
2214 ////////////////////////////////////////////////////////////////////////////////
2215 /// Computes the density function of F-distribution
2216 /// (probability function, integral of density, is computed in FDistI).
2217 ///
2218 /// Parameters N and M stand for degrees of freedom of chi-squares
2219 /// mentioned above parameter F is the actual variable x of the
2220 /// density function p(x) and the point at which the density function
2221 /// is calculated.
2222 ///
2223 /// About F distribution:
2224 /// F-distribution arises in testing whether two random samples
2225 /// have the same variance. It is the ratio of two chi-square
2226 /// distributions, with N and M degrees of freedom respectively,
2227 /// where each chi-square is first divided by it's number of degrees
2228 /// of freedom.
2229 /// Implementation by Anna Kreshuk.
2230 
2232 {
2234 }
2235 
2236 ////////////////////////////////////////////////////////////////////////////////
2237 /// Calculates the cumulative distribution function of F-distribution,
2238 /// this function occurs in the statistical test of whether two observed
2239 /// samples have the same variance. For this test a certain statistic F,
2240 /// the ratio of observed dispersion of the first sample to that of the
2241 /// second sample, is calculated. N and M stand for numbers of degrees
2242 /// of freedom in the samples 1-FDistI() is the significance level at
2243 /// which the hypothesis "1 has smaller variance than 2" can be rejected.
2244 /// A small numerical value of 1 - FDistI() implies a very significant
2245 /// rejection, in turn implying high confidence in the hypothesis
2246 /// "1 has variance greater than 2".
2247 /// Implementation by Anna Kreshuk.
2248 
2250 {
2252  return fi;
2253 }
2254 
2255 ////////////////////////////////////////////////////////////////////////////////
2256 /// Computes the density function of Gamma distribution at point x.
2257 /// gamma - shape parameter
2258 /// mu - location parameter
2259 /// beta - scale parameter
2260 ///
2261 /// The definition can be found in "Engineering Statistics Handbook" on site
2262 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
2263 /// use now implementation in ROOT::Math::gamma_pdf
2264 ///
2265 /// Begin_Macro
2266 /// {
2267 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2268 ///
2269 /// c1->SetLogy();
2270 /// c1->SetGridx();
2271 /// c1->SetGridy();
2272 ///
2273 /// TF1 *gdist = new TF1("gdist", "TMath::GammaDist(x, [0], [1], [2])", 0, 10);
2274 ///
2275 /// gdist->SetParameters(0.5, 0., 1.);
2276 /// gdist->SetLineColor(2);
2277 /// TF1 *gdist1 = gdist->DrawCopy("L");
2278 /// gdist->SetParameters(1.0, 0., 1.);
2279 /// gdist->SetLineColor(3);
2280 /// TF1 *gdist2 = gdist->DrawCopy("LSAME");
2281 /// gdist->SetParameters(2.0, 0., 1.);
2282 /// gdist->SetLineColor(4);
2283 /// TF1 *gdist3 = gdist->DrawCopy("LSAME");
2284 /// gdist->SetParameters(5.0, 0., 1.);
2285 /// gdist->SetLineColor(6);
2286 /// TF1 *gdist4 = gdist->DrawCopy("LSAME");
2287 ///
2288 /// legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2289 /// legend->AddEntry(gdist1, "gamma = 0.5, mu = 0, beta = 1", "L");
2290 /// legend->AddEntry(gdist2, "gamma = 1.0, mu = 0, beta = 1", "L");
2291 /// legend->AddEntry(gdist3, "gamma = 2.0, mu = 0, beta = 1", "L");
2292 /// legend->AddEntry(gdist4, "gamma = 5.0, mu = 0, beta = 1", "L");
2293 /// legend->Draw();
2294 /// }
2295 /// End_Macro
2296 
2298 {
2299  if ((x<mu) || (gamma<=0) || (beta <=0)) {
2300  Error("TMath::GammaDist", "illegal parameter values x = %f , gamma = %f beta = %f",x,gamma,beta);
2301  return 0;
2302  }
2303  return ::ROOT::Math::gamma_pdf(x, gamma, beta, mu);
2304 }
2305 
2306 ////////////////////////////////////////////////////////////////////////////////
2307 /// Computes the probability density function of Laplace distribution
2308 /// at point x, with location parameter alpha and shape parameter beta.
2309 /// By default, alpha=0, beta=1
2310 /// This distribution is known under different names, most common is
2311 /// double exponential distribution, but it also appears as
2312 /// the two-tailed exponential or the bilateral exponential distribution
2313 
2315 {
2316  Double_t temp;
2317  temp = TMath::Exp(-TMath::Abs((x-alpha)/beta));
2318  temp /= (2.*beta);
2319  return temp;
2320 }
2321 
2322 ////////////////////////////////////////////////////////////////////////////////
2323 /// Computes the distribution function of Laplace distribution
2324 /// at point x, with location parameter alpha and shape parameter beta.
2325 /// By default, alpha=0, beta=1
2326 /// This distribution is known under different names, most common is
2327 /// double exponential distribution, but it also appears as
2328 /// the two-tailed exponential or the bilateral exponential distribution
2329 
2331 {
2332  Double_t temp;
2333  if (x<=alpha){
2334  temp = 0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2335  } else {
2336  temp = 1-0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2337  }
2338  return temp;
2339 }
2340 
2341 ////////////////////////////////////////////////////////////////////////////////
2342 /// Computes the density of LogNormal distribution at point x.
2343 /// Variable X has lognormal distribution if Y=Ln(X) has normal distribution
2344 /// - sigma is the shape parameter
2345 /// - theta is the location parameter
2346 /// - m is the scale parameter
2347 ///
2348 /// The formula was taken from "Engineering Statistics Handbook" on site
2349 /// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm
2350 /// Implementation using ROOT::Math::lognormal_pdf
2351 ///
2352 /// Begin_Macro
2353 /// {
2354 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2355 ///
2356 /// c1->SetLogy();
2357 /// c1->SetGridx();
2358 /// c1->SetGridy();
2359 ///
2360 /// TF1 *logn = new TF1("logn", "TMath::LogNormal(x, [0], [1], [2])", 0, 5);
2361 /// logn->SetMinimum(1e-3);
2362 ///
2363 /// logn->SetParameters(0.5, 0., 1.);
2364 /// logn->SetLineColor(2);
2365 /// TF1 *logn1 = logn->DrawCopy("L");
2366 /// logn->SetParameters(1.0, 0., 1.);
2367 /// logn->SetLineColor(3);
2368 /// TF1 *logn2 = logn->DrawCopy("LSAME");
2369 /// logn->SetParameters(2.0, 0., 1.);
2370 /// logn->SetLineColor(4);
2371 /// TF1 *logn3 = logn->DrawCopy("LSAME");
2372 /// logn->SetParameters(5.0, 0., 1.);
2373 /// logn->SetLineColor(6);
2374 /// TF1 *logn4 = logn->DrawCopy("LSAME");
2375 ///
2376 /// legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2377 /// legend->AddEntry(logn1, "sigma = 0.5, theta = 0, m = 1", "L");
2378 /// legend->AddEntry(logn2, "sigma = 1.0, theta = 0, m = 1", "L");
2379 /// legend->AddEntry(logn3, "sigma = 2.0, theta = 0, m = 1", "L");
2380 /// legend->AddEntry(logn4, "sigma = 5.0, theta = 0, m = 1", "L");
2381 /// legend->Draw();
2382 /// }
2383 /// End_Macro
2384 
2386 {
2387  if ((x<theta) || (sigma<=0) || (m<=0)) {
2388  Error("TMath::Lognormal", "illegal parameter values");
2389  return 0;
2390  }
2391  // lognormal_pdf uses same definition of http://en.wikipedia.org/wiki/Log-normal_distribution
2392  // where mu = log(m)
2393 
2394  return ::ROOT::Math::lognormal_pdf(x, TMath::Log(m), sigma, theta);
2395 
2396 }
2397 
2398 ////////////////////////////////////////////////////////////////////////////////
2399 /// Computes quantiles for standard normal distribution N(0, 1)
2400 /// at probability p
2401 /// ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.
2402 
2404 {
2405  if ((p<=0)||(p>=1)) {
2406  Error("TMath::NormQuantile", "probability outside (0, 1)");
2407  return 0;
2408  }
2409 
2410  Double_t a0 = 3.3871328727963666080e0;
2411  Double_t a1 = 1.3314166789178437745e+2;
2412  Double_t a2 = 1.9715909503065514427e+3;
2413  Double_t a3 = 1.3731693765509461125e+4;
2414  Double_t a4 = 4.5921953931549871457e+4;
2415  Double_t a5 = 6.7265770927008700853e+4;
2416  Double_t a6 = 3.3430575583588128105e+4;
2417  Double_t a7 = 2.5090809287301226727e+3;
2418  Double_t b1 = 4.2313330701600911252e+1;
2419  Double_t b2 = 6.8718700749205790830e+2;
2420  Double_t b3 = 5.3941960214247511077e+3;
2421  Double_t b4 = 2.1213794301586595867e+4;
2422  Double_t b5 = 3.9307895800092710610e+4;
2423  Double_t b6 = 2.8729085735721942674e+4;
2424  Double_t b7 = 5.2264952788528545610e+3;
2425  Double_t c0 = 1.42343711074968357734e0;
2426  Double_t c1 = 4.63033784615654529590e0;
2427  Double_t c2 = 5.76949722146069140550e0;
2428  Double_t c3 = 3.64784832476320460504e0;
2429  Double_t c4 = 1.27045825245236838258e0;
2430  Double_t c5 = 2.41780725177450611770e-1;
2431  Double_t c6 = 2.27238449892691845833e-2;
2432  Double_t c7 = 7.74545014278341407640e-4;
2433  Double_t d1 = 2.05319162663775882187e0;
2434  Double_t d2 = 1.67638483018380384940e0;
2435  Double_t d3 = 6.89767334985100004550e-1;
2436  Double_t d4 = 1.48103976427480074590e-1;
2437  Double_t d5 = 1.51986665636164571966e-2;
2438  Double_t d6 = 5.47593808499534494600e-4;
2439  Double_t d7 = 1.05075007164441684324e-9;
2440  Double_t e0 = 6.65790464350110377720e0;
2441  Double_t e1 = 5.46378491116411436990e0;
2442  Double_t e2 = 1.78482653991729133580e0;
2443  Double_t e3 = 2.96560571828504891230e-1;
2444  Double_t e4 = 2.65321895265761230930e-2;
2445  Double_t e5 = 1.24266094738807843860e-3;
2446  Double_t e6 = 2.71155556874348757815e-5;
2447  Double_t e7 = 2.01033439929228813265e-7;
2448  Double_t f1 = 5.99832206555887937690e-1;
2449  Double_t f2 = 1.36929880922735805310e-1;
2450  Double_t f3 = 1.48753612908506148525e-2;
2451  Double_t f4 = 7.86869131145613259100e-4;
2452  Double_t f5 = 1.84631831751005468180e-5;
2453  Double_t f6 = 1.42151175831644588870e-7;
2454  Double_t f7 = 2.04426310338993978564e-15;
2455 
2456  Double_t split1 = 0.425;
2457  Double_t split2=5.;
2458  Double_t konst1=0.180625;
2459  Double_t konst2=1.6;
2460 
2461  Double_t q, r, quantile;
2462  q=p-0.5;
2463  if (TMath::Abs(q)<split1) {
2464  r=konst1-q*q;
2465  quantile = q* (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
2466  * r + a2) * r + a1) * r + a0) /
2467  (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
2468  * r + b2) * r + b1) * r + 1.);
2469  } else {
2470  if(q<0) r=p;
2471  else r=1-p;
2472  //error case
2473  if (r<=0)
2474  quantile=0;
2475  else {
2476  r=TMath::Sqrt(-TMath::Log(r));
2477  if (r<=split2) {
2478  r=r-konst2;
2479  quantile=(((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
2480  * r + c2) * r + c1) * r + c0) /
2481  (((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
2482  * r + d2) * r + d1) * r + 1);
2483  } else{
2484  r=r-split2;
2485  quantile=(((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
2486  * r + e2) * r + e1) * r + e0) /
2487  (((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
2488  * r + f2) * r + f1) * r + 1);
2489  }
2490  if (q<0) quantile=-quantile;
2491  }
2492  }
2493  return quantile;
2494 }
2495 
2496 ////////////////////////////////////////////////////////////////////////////////
2497 /// Simple recursive algorithm to find the permutations of
2498 /// n natural numbers, not necessarily all distinct
2499 /// adapted from CERNLIB routine PERMU.
2500 /// The input array has to be initialised with a non descending
2501 /// sequence. The method returns kFALSE when
2502 /// all combinations are exhausted.
2503 
2505 {
2506  Int_t i,itmp;
2507  Int_t i1=-1;
2508 
2509  // find rightmost upward transition
2510  for(i=n-2; i>-1; i--) {
2511  if(a[i]<a[i+1]) {
2512  i1=i;
2513  break;
2514  }
2515  }
2516  // no more upward transitions, end of the story
2517  if(i1==-1) return kFALSE;
2518  else {
2519  // find lower right element higher than the lower
2520  // element of the upward transition
2521  for(i=n-1;i>i1;i--) {
2522  if(a[i] > a[i1]) {
2523  // swap the two
2524  itmp=a[i1];
2525  a[i1]=a[i];
2526  a[i]=itmp;
2527  break;
2528  }
2529  }
2530  // order the rest, in fact just invert, as there
2531  // are only downward transitions from here on
2532  for(i=0;i<(n-i1-1)/2;i++) {
2533  itmp=a[i1+i+1];
2534  a[i1+i+1]=a[n-i-1];
2535  a[n-i-1]=itmp;
2536  }
2537  }
2538  return kTRUE;
2539 }
2540 
2541 ////////////////////////////////////////////////////////////////////////////////
2542 /// Computes density function for Student's t- distribution
2543 /// (the probability function (integral of density) is computed in StudentI).
2544 ///
2545 /// First parameter stands for x - the actual variable of the
2546 /// density function p(x) and the point at which the density is calculated.
2547 /// Second parameter stands for number of degrees of freedom.
2548 ///
2549 /// About Student distribution:
2550 /// Student's t-distribution is used for many significance tests, for example,
2551 /// for the Student's t-tests for the statistical significance of difference
2552 /// between two sample means and for confidence intervals for the difference
2553 /// between two population means.
2554 ///
2555 /// Example: suppose we have a random sample of size n drawn from normal
2556 /// distribution with mean Mu and st.deviation Sigma. Then the variable
2557 ///
2558 /// t = (sample_mean - Mu)/(sample_deviation / sqrt(n))
2559 ///
2560 /// has Student's t-distribution with n-1 degrees of freedom.
2561 ///
2562 /// NOTE that this function's second argument is number of degrees of freedom,
2563 /// not the sample size.
2564 ///
2565 /// As the number of degrees of freedom grows, t-distribution approaches
2566 /// Normal(0,1) distribution.
2567 /// Implementation by Anna Kreshuk.
2568 
2570 {
2571  if (ndf < 1) {
2572  return 0;
2573  }
2574 
2575  Double_t r = ndf;
2576  Double_t rh = 0.5*r;
2577  Double_t rh1 = rh + 0.5;
2578  Double_t denom = TMath::Sqrt(r*TMath::Pi())*TMath::Gamma(rh)*TMath::Power(1+T*T/r, rh1);
2579  return TMath::Gamma(rh1)/denom;
2580 }
2581 
2582 ////////////////////////////////////////////////////////////////////////////////
2583 /// Calculates the cumulative distribution function of Student's
2584 /// t-distribution second parameter stands for number of degrees of freedom,
2585 /// not for the number of samples
2586 /// if x has Student's t-distribution, the function returns the probability of
2587 /// x being less than T.
2588 /// Implementation by Anna Kreshuk.
2589 
2591 {
2592  Double_t r = ndf;
2593 
2594  Double_t si = (T>0) ?
2595  (1 - 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5)) :
2596  0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5);
2597  return si;
2598 }
2599 
2600 ////////////////////////////////////////////////////////////////////////////////
2601 /// Computes quantiles of the Student's t-distribution
2602 /// 1st argument is the probability, at which the quantile is computed
2603 /// 2nd argument - the number of degrees of freedom of the
2604 /// Student distribution
2605 /// When the 3rd argument lower_tail is kTRUE (default)-
2606 /// the algorithm returns such x0, that
2607 /// P(x < x0)=p
2608 /// upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that
2609 /// P(x > x0)=p
2610 /// the algorithm was taken from
2611 /// G.W.Hill, "Algorithm 396, Student's t-quantiles"
2612 /// "Communications of the ACM", 13(10), October 1970
2613 
2615 {
2616  Double_t quantile;
2617  Double_t temp;
2618  Bool_t neg;
2619  Double_t q;
2620  if (ndf<1 || p>=1 || p<=0) {
2621  Error("TMath::StudentQuantile", "illegal parameter values");
2622  return 0;
2623  }
2624  if ((lower_tail && p>0.5)||(!lower_tail && p<0.5)){
2625  neg=kFALSE;
2626  q=2*(lower_tail ? (1-p) : p);
2627  } else {
2628  neg=kTRUE;
2629  q=2*(lower_tail? p : (1-p));
2630  }
2631 
2632  if ((ndf-1)<1e-8) {
2633  temp=TMath::PiOver2()*q;
2634  quantile = TMath::Cos(temp)/TMath::Sin(temp);
2635  } else {
2636  if ((ndf-2)<1e-8) {
2637  quantile = TMath::Sqrt(2./(q*(2-q))-2);
2638  } else {
2639  Double_t a=1./(ndf-0.5);
2640  Double_t b=48./(a*a);
2641  Double_t c=((20700*a/b -98)*a-16)*a+96.36;
2642  Double_t d=((94.5/(b+c)-3.)/b+1)*TMath::Sqrt(a*TMath::PiOver2())*ndf;
2643  Double_t x=q*d;
2644  Double_t y=TMath::Power(x, (2./ndf));
2645  if (y>0.05+a){
2646  //asymptotic inverse expansion about normal
2647  x=TMath::NormQuantile(q*0.5);
2648  y=x*x;
2649  if (ndf<5) c+=0.3*(ndf-4.5)*(x+0.6);
2650  c+=(((0.05*d*x-5.)*x-7.)*x-2.)*x +b;
2651  y=(((((0.4*y+6.3)*y+36.)*y+94.5)/c - y-3.)/b+1)*x;
2652  y=a*y*y;
2653  if(y>0.002) y = TMath::Exp(y)-1;
2654  else y += 0.5*y*y;
2655  } else {
2656  y=((1./(((ndf+6.)/(ndf*y)-0.089*d-0.822)*(ndf+2.)*3)+0.5/(ndf+4.))*y-1.)*
2657  (ndf+1.)/(ndf+2.)+1/y;
2658  }
2659  quantile = TMath::Sqrt(ndf*y);
2660  }
2661  }
2662  if(neg) quantile=-quantile;
2663  return quantile;
2664 }
2665 
2666 ////////////////////////////////////////////////////////////////////////////////
2667 /// Returns the value of the Vavilov density function
2668 /// Parameters: 1st - the point were the density function is evaluated
2669 /// 2nd - value of kappa (distribution parameter)
2670 /// 3rd - value of beta2 (distribution parameter)
2671 /// The algorithm was taken from the CernLib function vavden(G115)
2672 /// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2673 /// Nucl.Instr. and Meth. B47(1990), 215-224
2674 /// Accuracy: quote from the reference above:
2675 /// "The resuls of our code have been compared with the values of the Vavilov
2676 /// density function computed numerically in an accurate way: our approximation
2677 /// shows a difference of less than 3% around the peak of the density function, slowly
2678 /// increasing going towards the extreme tails to the right and to the left"
2679 ///
2680 /// Begin_Macro
2681 /// {
2682 /// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2683 ///
2684 /// c1->SetGridx();
2685 /// c1->SetGridy();
2686 ///
2687 /// TF1 *vavilov = new TF1("vavilov", "TMath::Vavilov(x, [0], [1])", -3, 11);
2688 ///
2689 /// vavilov->SetParameters(0.5, 0.);
2690 /// vavilov->SetLineColor(2);
2691 /// TF1 *vavilov1 = vavilov->DrawCopy("L");
2692 /// vavilov->SetParameters(0.3, 0.);
2693 /// vavilov->SetLineColor(3);
2694 /// TF1 *vavilov2 = vavilov->DrawCopy("LSAME");
2695 /// vavilov->SetParameters(0.2, 0.);
2696 /// vavilov->SetLineColor(4);
2697 /// TF1 *vavilov3 = vavilov->DrawCopy("LSAME");
2698 /// vavilov->SetParameters(0.1, 0.);
2699 /// vavilov->SetLineColor(6);
2700 /// TF1 *vavilov4 = vavilov->DrawCopy("LSAME");
2701 ///
2702 /// legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2703 /// legend->AddEntry(vavilov1, "kappa = 0.5, beta2 = 0", "L");
2704 /// legend->AddEntry(vavilov2, "kappa = 0.3, beta2 = 0", "L");
2705 /// legend->AddEntry(vavilov3, "kappa = 0.2, beta2 = 0", "L");
2706 /// legend->AddEntry(vavilov4, "kappa = 0.1, beta2 = 0", "L");
2707 /// legend->Draw();
2708 /// }
2709 /// End_Macro
2710 
2712 {
2713  Double_t *ac = new Double_t[14];
2714  Double_t *hc = new Double_t[9];
2715 
2716  Int_t itype;
2717  Int_t npt;
2718  TMath::VavilovSet(kappa, beta2, 0, 0, ac, hc, itype, npt);
2719  Double_t v = TMath::VavilovDenEval(x, ac, hc, itype);
2720  delete [] ac;
2721  delete [] hc;
2722  return v;
2723 }
2724 
2725 ////////////////////////////////////////////////////////////////////////////////
2726 ///Returns the value of the Vavilov distribution function
2727 ///Parameters: 1st - the point were the density function is evaluated
2728 /// 2nd - value of kappa (distribution parameter)
2729 /// 3rd - value of beta2 (distribution parameter)
2730 ///The algorithm was taken from the CernLib function vavden(G115)
2731 ///Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2732 ///Nucl.Instr. and Meth. B47(1990), 215-224
2733 ///Accuracy: quote from the reference above:
2734 ///"The resuls of our code have been compared with the values of the Vavilov
2735 ///density function computed numerically in an accurate way: our approximation
2736 ///shows a difference of less than 3% around the peak of the density function, slowly
2737 ///increasing going towards the extreme tails to the right and to the left"
2738 
2740 {
2741  Double_t *ac = new Double_t[14];
2742  Double_t *hc = new Double_t[9];
2743  Double_t *wcm = new Double_t[201];
2744  Int_t itype;
2745  Int_t npt;
2746  Int_t k;
2747  Double_t xx, v;
2748  TMath::VavilovSet(kappa, beta2, 1, wcm, ac, hc, itype, npt);
2749  if (x < ac[0]) v = 0;
2750  else if (x >=ac[8]) v = 1;
2751  else {
2752  xx = x - ac[0];
2753  k = Int_t(xx*ac[10]);
2754  v = TMath::Min(wcm[k] + (xx - k*ac[9])*(wcm[k+1]-wcm[k])*ac[10], 1.);
2755  }
2756  delete [] ac;
2757  delete [] hc;
2758  delete [] wcm;
2759  return v;
2760 }
2761 
2762 ////////////////////////////////////////////////////////////////////////////////
2763 ///Returns the value of the Landau distribution function at point x.
2764 ///The algorithm was taken from the Cernlib function dislan(G110)
2765 ///Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
2766 ///distribution", Computer Phys.Comm., 31(1984), 97-111
2767 
2769 {
2771 }
2772 
2773 
2774 ////////////////////////////////////////////////////////////////////////////////
2775 ///Internal function, called by Vavilov and VavilovI
2776 
2777 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)
2778 {
2779 
2780  Double_t BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
2781  BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
2782  BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
2783 
2784  Double_t FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
2785  FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
2786  FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
2787 
2788  Double_t FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
2789 
2790  Double_t EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
2791  0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
2792 
2793  Double_t U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
2794  0. , 0.24880692e-1, 0.47404356e-2,
2795  -0.74445130e-3, 0.73225731e-2, 0. ,
2796  0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
2797 
2798  Double_t U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
2799  0. , 0.42127077e-1, 0.73167928e-2,
2800  -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
2801  0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
2802 
2803  Double_t U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
2804  0. , 0.23834176e-1, 0.21624675e-2,
2805  -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
2806  0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
2807 
2808  Double_t U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
2809  -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
2810  0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
2811  -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
2812 
2813  Double_t U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
2814  0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
2815  0.48736023e-3, 0.34850854e-2, 0. ,
2816  -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
2817 
2818  Double_t U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
2819  -0.30188807e-2, -0.84479939e-3, 0. ,
2820  0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
2821  -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
2822 
2823  Double_t U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
2824  -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
2825  0. , 0.50505298e+0};
2826  Double_t U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
2827  -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
2828  -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
2829  0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
2830 
2831  Double_t V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
2832  0. , 0.45091424e-1, 0.80559636e-2,
2833  -0.38974523e-2, 0. , -0.30634124e-2,
2834  0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
2835 
2836  Double_t V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
2837  0. , 0.12693873e+0, 0.22999801e-1,
2838  -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
2839  0. , 0.19716857e-1, 0.32596742e-2};
2840 
2841  Double_t V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
2842  -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
2843  -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
2844  0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
2845 
2846  Double_t V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
2847  0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
2848  0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
2849  -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
2850 
2851  Double_t V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
2852  -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
2853  0.56856517e-2, -0.13438433e-1, 0. ,
2854  0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
2855 
2856  Double_t V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
2857  0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
2858  0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
2859  -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
2860 
2861  Double_t V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
2862  0. , -0.12218009e+1, -0.32324120e+0,
2863  -0.27373591e-1, 0.12173464e+0, 0. ,
2864  0. , 0.40917471e-1};
2865 
2866  Double_t V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
2867  0. , -0.18160479e+1, -0.50919193e+0,
2868  -0.51384654e-1, 0.21413992e+0, 0. ,
2869  0. , 0.66596366e-1};
2870 
2871  Double_t W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
2872  -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
2873  -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
2874  0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
2875 
2876  Double_t W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
2877  -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
2878  0. , 0.23020158e-1, 0.50574313e-2,
2879  0.94550140e-2, 0.19300232e-1};
2880 
2881  Double_t W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
2882  -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
2883  -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
2884  0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
2885 
2886  Double_t W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
2887  0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
2888  0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
2889  -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
2890 
2891  Double_t W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
2892  0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
2893  0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
2894  -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
2895 
2896  Double_t W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
2897  0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
2898  0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
2899  -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
2900 
2901  Double_t W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
2902  0. , -0.14540925e+1, -0.39529833e+0,
2903  -0.44293243e-1, 0.88741049e-1};
2904 
2905  itype = 0;
2906  if (rkappa <0.01 || rkappa >12) {
2907  Error("Vavilov distribution", "illegal value of kappa");
2908  return;
2909  }
2910 
2911  Double_t DRK[6];
2912  Double_t DSIGM[6];
2913  Double_t ALFA[8];
2914  Int_t j;
2915  Double_t x, y, xx, yy, x2, x3, y2, y3, xy, p2, p3, q2, q3, pq;
2916  if (rkappa >= 0.29) {
2917  itype = 1;
2918  npt = 100;
2919  Double_t wk = 1./TMath::Sqrt(rkappa);
2920 
2921  AC[0] = (-0.032227*beta2-0.074275)*rkappa + (0.24533*beta2+0.070152)*wk + (-0.55610*beta2-3.1579);
2922  AC[8] = (-0.013483*beta2-0.048801)*rkappa + (-1.6921*beta2+8.3656)*wk + (-0.73275*beta2-3.5226);
2923  DRK[1] = wk*wk;
2924  DSIGM[1] = TMath::Sqrt(rkappa/(1-0.5*beta2));
2925  for (j=1; j<=4; j++) {
2926  DRK[j+1] = DRK[1]*DRK[j];
2927  DSIGM[j+1] = DSIGM[1]*DSIGM[j];
2928  ALFA[j+1] = (FNINV[j]-beta2*FNINV[j+1])*DRK[j];
2929  }
2930  HC[0]=TMath::Log(rkappa)+beta2+0.42278434;
2931  HC[1]=DSIGM[1];
2932  HC[2]=ALFA[3]*DSIGM[3];
2933  HC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
2934  HC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*HC[2];
2935  HC[5]=HC[2]*HC[2];
2936  HC[6]=HC[2]*HC[3];
2937  HC[7]=HC[2]*HC[5];
2938  for (j=2; j<=7; j++)
2939  HC[j]*=EDGEC[j];
2940  HC[8]=0.39894228*HC[1];
2941  }
2942  else if (rkappa >=0.22) {
2943  itype = 2;
2944  npt = 150;
2945  x = 1+(rkappa-BKMXX3)*FBKX3;
2946  y = 1+(TMath::Sqrt(beta2)-BKMXY3)*FBKY3;
2947  xx = 2*x;
2948  yy = 2*y;
2949  x2 = xx*x-1;
2950  x3 = xx*x2-x;
2951  y2 = yy*y-1;
2952  y3 = yy*y2-y;
2953  xy = x*y;
2954  p2 = x2*y;
2955  p3 = x3*y;
2956  q2 = y2*x;
2957  q3 = y3*x;
2958  pq = x2*y2;
2959  AC[1] = W1[1] + W1[2]*x + W1[4]*x3 + W1[5]*y + W1[6]*y2 + W1[7]*y3 +
2960  W1[8]*xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
2961  AC[2] = W2[1] + W2[2]*x + W2[3]*x2 + W2[4]*x3 + W2[5]*y + W2[6]*y2 +
2962  W2[8]*xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
2963  AC[3] = W3[1] + W3[3]*x2 + W3[4]*x3 + W3[5]*y + W3[6]*y2 + W3[7]*y3 +
2964  W3[8]*xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
2965  AC[4] = W4[1] + W4[2]*x + W4[3]*x2 + W4[4]*x3 + W4[5]*y + W4[6]*y2 + W4[7]*y3 +
2966  W4[8]*xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
2967  AC[5] = W5[1] + W5[2]*x + W5[4]*x3 + W5[5]*y + W5[6]*y2 + W5[7]*y3 +
2968  W5[8]*xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
2969  AC[6] = W6[1] + W6[2]*x + W6[3]*x2 + W6[4]*x3 + W6[5]*y + W6[6]*y2 + W6[7]*y3 +
2970  W6[8]*xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
2971  AC[8] = W8[1] + W8[2]*x + W8[3]*x2 + W8[5]*y + W8[6]*y2 + W8[7]*y3 + W8[8]*xy;
2972  AC[0] = -3.05;
2973  } else if (rkappa >= 0.12) {
2974  itype = 3;
2975  npt = 200;
2976  x = 1 + (rkappa-BKMXX2)*FBKX2;
2977  y = 1 + (TMath::Sqrt(beta2)-BKMXY2)*FBKY2;
2978  xx = 2*x;
2979  yy = 2*y;
2980  x2 = xx*x-1;
2981  x3 = xx*x2-x;
2982  y2 = yy*y-1;
2983  y3 = yy*y2-y;
2984  xy = x*y;
2985  p2 = x2*y;
2986  p3 = x3*y;
2987  q2 = y2*x;
2988  q3 = y3*x;
2989  pq = x2*y2;
2990  AC[1] = V1[1] + V1[2]*x + V1[3]*x2 + V1[5]*y + V1[6]*y2 + V1[7]*y3 +
2991  V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
2992  AC[2] = V2[1] + V2[2]*x + V2[3]*x2 + V2[5]*y + V2[6]*y2 + V2[7]*y3 +
2993  V2[8]*xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
2994  AC[3] = V3[1] + V3[2]*x + V3[3]*x2 + V3[4]*x3 + V3[5]*y + V3[6]*y2 + V3[7]*y3 +
2995  V3[8]*xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
2996  AC[4] = V4[1] + V4[2]*x + V4[3]*x2 + V4[4]*x3 + V4[5]*y + V4[6]*y2 + V4[7]*y3 +
2997  V4[8]*xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
2998  AC[5] = V5[1] + V5[2]*x + V5[3]*x2 + V5[4]*x3 + V5[5]*y + V5[6]*y2 + V5[7]*y3 +
2999  V5[8]*xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
3000  AC[6] = V6[1] + V6[2]*x + V6[3]*x2 + V6[4]*x3 + V6[5]*y + V6[6]*y2 + V6[7]*y3 +
3001  V6[8]*xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
3002  AC[7] = V7[1] + V7[2]*x + V7[3]*x2 + V7[5]*y + V7[6]*y2 + V7[7]*y3 +
3003  V7[8]*xy + V7[11]*q2;
3004  AC[8] = V8[1] + V8[2]*x + V8[3]*x2 + V8[5]*y + V8[6]*y2 + V8[7]*y3 +
3005  V8[8]*xy + V8[11]*q2;
3006  AC[0] = -3.04;
3007  } else {
3008  itype = 4;
3009  if (rkappa >=0.02) itype = 3;
3010  npt = 200;
3011  x = 1+(rkappa-BKMXX1)*FBKX1;
3012  y = 1+(TMath::Sqrt(beta2)-BKMXY1)*FBKY1;
3013  xx = 2*x;
3014  yy = 2*y;
3015  x2 = xx*x-1;
3016  x3 = xx*x2-x;
3017  y2 = yy*y-1;
3018  y3 = yy*y2-y;
3019  xy = x*y;
3020  p2 = x2*y;
3021  p3 = x3*y;
3022  q2 = y2*x;
3023  q3 = y3*x;
3024  pq = x2*y2;
3025  if (itype==3){
3026  AC[1] = U1[1] + U1[2]*x + U1[3]*x2 + U1[5]*y + U1[6]*y2 + U1[7]*y3 +
3027  U1[8]*xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
3028  AC[2] = U2[1] + U2[2]*x + U2[3]*x2 + U2[5]*y + U2[6]*y2 + U2[7]*y3 +
3029  U2[8]*xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
3030  AC[3] = U3[1] + U3[2]*x + U3[3]*x2 + U3[5]*y + U3[6]*y2 + U3[7]*y3 +
3031  U3[8]*xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
3032  AC[4] = U4[1] + U4[2]*x + U4[3]*x2 + U4[4]*x3 + U4[5]*y + U4[6]*y2 + U4[7]*y3 +
3033  U4[8]*xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
3034  AC[5] = U5[1] + U5[2]*x + U5[3]*x2 + U5[4]*x3 + U5[5]*y + U5[6]*y2 + U5[7]*y3 +
3035  U5[8]*xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
3036  AC[6] = U6[1] + U6[2]*x + U6[3]*x2 + U6[4]*x3 + U6[5]*y + U6[7]*y3 +
3037  U6[8]*xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
3038  AC[7] = U7[1] + U7[2]*x + U7[3]*x2 + U7[4]*x3 + U7[5]*y + U7[6]*y2 + U7[8]*xy;
3039  }
3040  AC[8] = U8[1] + U8[2]*x + U8[3]*x2 + U8[4]*x3 + U8[5]*y + U8[6]*y2 + U8[7]*y3 +
3041  U8[8]*xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
3042  AC[0] = -3.03;
3043  }
3044 
3045  AC[9] = (AC[8] - AC[0])/npt;
3046  AC[10] = 1./AC[9];
3047  if (itype == 3) {
3048  x = (AC[7]-AC[8])/(AC[7]*AC[8]);
3049  y = 1./TMath::Log(AC[8]/AC[7]);
3050  p2 = AC[7]*AC[7];
3051  AC[11] = p2*(AC[1]*TMath::Exp(-AC[2]*(AC[7]+AC[5]*p2)-
3052  AC[3]*TMath::Exp(-AC[4]*(AC[7]+AC[6]*p2)))-0.045*y/AC[7])/(1+x*y*AC[7]);
3053  AC[12] = (0.045+x*AC[11])*y;
3054  }
3055  if (itype == 4) AC[13] = 0.995/LandauI(AC[8]);
3056 
3057  if (mode==0) return;
3058  //
3059  x = AC[0];
3060  WCM[0] = 0;
3061  Double_t fl, fu;
3062  Int_t k;
3063  fl = TMath::VavilovDenEval(x, AC, HC, itype);
3064  for (k=1; k<=npt; k++) {
3065  x += AC[9];
3066  fu = TMath::VavilovDenEval(x, AC, HC, itype);
3067  WCM[k] = WCM[k-1] + fl + fu;
3068  fl = fu;
3069  }
3070  x = 0.5*AC[9];
3071  for (k=1; k<=npt; k++)
3072  WCM[k]*=x;
3073 }
3074 
3075 ////////////////////////////////////////////////////////////////////////////////
3076 ///Internal function, called by Vavilov and VavilovSet
3077 
3079 {
3080  Double_t v = 0;
3081  if (rlam < AC[0] || rlam > AC[8])
3082  return 0;
3083  Int_t k;
3084  Double_t x, fn, s;
3085  Double_t h[10];
3086  if (itype ==1 ) {
3087  fn = 1;
3088  x = (rlam + HC[0])*HC[1];
3089  h[1] = x;
3090  h[2] = x*x -1;
3091  for (k=2; k<=8; k++) {
3092  fn++;
3093  h[k+1] = x*h[k]-fn*h[k-1];
3094  }
3095  s = 1 + HC[7]*h[9];
3096  for (k=2; k<=6; k++)
3097  s+=HC[k]*h[k+1];
3098  v = HC[8]*TMath::Exp(-0.5*x*x)*TMath::Max(s, 0.);
3099  }
3100  else if (itype == 2) {
3101  x = rlam*rlam;
3102  v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x) - AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3103  }
3104  else if (itype == 3) {
3105  if (rlam < AC[7]) {
3106  x = rlam*rlam;
3107  v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x)-AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3108  } else {
3109  x = 1./rlam;
3110  v = (AC[11]*x + AC[12])*x;
3111  }
3112  }
3113  else if (itype == 4) {
3114  v = AC[13]*TMath::Landau(rlam);
3115  }
3116  return v;
3117 }
3118 
3119 
3120 //explicity instantiate template functions from VecCore
3121 #ifdef R__HAS_VECCORE
3122 #include <Math/Types.h>
3128 // missing in veccore
3129 // template ROOT::Double_v vecCore::math::ACos(const ROOT::Double_v & x);
3130 // template ROOT::Double_v vecCore::math::Sinh(const ROOT::Double_v & x);
3131 // template ROOT::Double_v vecCore::math::Cosh(const ROOT::Double_v & x);
3132 // template ROOT::Double_v vecCore::math::Tanh(const ROOT::Double_v & x);
3133 // template ROOT::Double_v vecCore::math::Cbrt(const ROOT::Double_v & x);
3134 #endif
double lognormal_pdf(double x, double m, double s, double x0=0)
Probability density function of the lognormal distribution.
Double_t ACosH(Double_t)
Definition: TMath.cxx: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:1560
constexpr Double_t HC()
Definition: TMath.h:168
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:2231
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:2173
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:2249
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 constexpr double pi
auto * m
Definition: textangle.C:8
static double p3(double t, double a, double b, double c, double d)
T1 Sign(T1 a, T2 b)
Definition: TMathBase.h:153
Double_t Log(Double_t x)
Definition: TMath.h:648
float Float_t
Definition: RtypesCore.h:53
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:2330
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:2614
Double_t NormQuantile(Double_t p)
Computes quantiles for standard normal distribution N(0, 1) at probability p ALGORITHM AS241 APPL...
Definition: TMath.cxx:2403
TH1 * h
Definition: legend2.C:5
static constexpr double km
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:2040
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:2314
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1112
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:96
#define N
static constexpr double bar
Basic string class.
Definition: TString.h:125
Double_t StruveH1(Double_t x)
Struve Functions of Order 1.
Definition: TMath.cxx:1816
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
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:1284
Int_t FloorNint(Double_t x)
Definition: TMath.h:602
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:1639
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
Double_t GamCf(Double_t a, Double_t x)
Computation of the incomplete gamma function P(a,x) via its continued fraction representation.
Definition: TMath.cxx: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:2590
UInt_t Hash(ECaseCompare cmp=kExact) const
Return hash value.
Definition: TString.cxx:616
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:627
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:1709
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition: TMath.cxx:2148
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:1178
Double_t ASinH(Double_t)
Definition: TMath.cxx:67
Double_t Log10(Double_t x)
Definition: TMath.h:651
Double_t BesselI1(Double_t x)
Compute the modified Bessel function I_1(x) for any real x.
Definition: TMath.cxx:1464
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:3078
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_t ATan2(Double_t, Double_t)
Definition: TMath.h:580
constexpr Double_t Pi()
Definition: TMath.h:40
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:1323
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
static constexpr double second
double landau_cdf(double x, double xi=1, double x0=0)
Cumulative distribution function of the Landau distribution (lower tail).
Bool_t Even(Long_t a)
Definition: TMathBase.h:100
Double_t DiLog(Double_t x)
The DiLogarithm function Code translated by R.Brun from CERNLIB DILOG function C332.
Definition: TMath.cxx:113
Double_t Double_v
Definition: Types.h:50
#define F(x, y, z)
double gamma(double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
ROOT::R::TRInterface & r
Definition: Object.C:4
Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2)
Returns the value of the Vavilov distribution function Parameters: 1st - the point were the density f...
Definition: TMath.cxx:2739
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:1675
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:187
auto * a
Definition: textangle.C:12
Double_t BetaDistI(Double_t x, Double_t p, Double_t q)
Computes the distribution function of the Beta distribution.
Definition: TMath.cxx:2058
Double_t StruveL1(Double_t x)
Modified Struve Function of Order 1.
Definition: TMath.cxx:1939
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:2777
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:2504
Double_t StruveL0(Double_t x)
Modified Struve Function of Order 0.
Definition: TMath.cxx:1893
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
#define NamespaceImp(name)
Definition: Rtypes.h:374
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2079
Double_t LandauI(Double_t x)
Returns the value of the Landau distribution function at point x.
Definition: TMath.cxx:2768
Double_t ACos(Double_t)
Definition: TMath.h:571
static constexpr double m2
Double_t BesselJ0(Double_t x)
Returns the Bessel function J0(x) for any real x.
Definition: TMath.cxx:1604
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:790
ULong_t Hash(const void *txt, Int_t ntxt)
Calculates hash index from any char string.
Definition: TMath.cxx:1378
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:2569
Double_t Beta(Double_t p, Double_t q)
Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
Definition: TMath.cxx:1980
Double_t Cos(Double_t)
Definition: TMath.h:550
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
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t BesselI0(Double_t x)
Compute the modified Bessel function I_0(x) for any real x.
Definition: TMath.cxx:1396
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:1531
long Long_t
Definition: RtypesCore.h:50
static constexpr double pi2
Double_t Exp(Double_t x)
Definition: TMath.h:621
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:1747
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:570
static constexpr double s
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:60
Double_t BesselK0(Double_t x)
Compute the modified Bessel function K_0(x) for positive real x.
Definition: TMath.cxx:1430
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:2104
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Double_t BetaCf(Double_t x, Double_t a, Double_t b)
Continued fraction evaluation by modified Lentz&#39;s method used in calculation of incomplete Beta funct...
Definition: TMath.cxx:1989
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Definition: TMath.h:1315
double chisquared_cdf_c(double x, double r, double x0=0)
Complement of the cumulative distribution function of the distribution with degrees of freedom (upp...
Double_t GammaDist(Double_t x, Double_t gamma, Double_t mu=0, Double_t beta=1)
Computes the density function of Gamma distribution at point x.
Definition: TMath.cxx:2297
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
Double_t Factorial(Int_t i)
Compute factorial(n).
Definition: TMath.cxx:250
Double_t Sin(Double_t)
Definition: TMath.h:547
TF1 * f1
Definition: legend1.C:11
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t ASin(Double_t)
Definition: TMath.h:565
Double_t BesselK1(Double_t x)
Compute the modified Bessel function K_1(x) for positive real x.
Definition: TMath.cxx:1499
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:2385
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)
Definition: TMath.cxx:1082
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
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:2071
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:2132
double inc_gamma(double a, double x)
Calculates the normalized (regularized) lower incomplete gamma function (lower integral) ...
static T Epsilon()
Definition: TMath.h:814
const Bool_t kTRUE
Definition: RtypesCore.h:87
Int_t Nint(T x)
Definition: TMath.h:606
return c3
Definition: legend3.C:15
const Int_t n
Definition: legend1.C:16
Double_t ATanH(Double_t)
Definition: TMath.cxx:93
constexpr Double_t PiOver2()
Definition: TMath.h:48
double log(double)
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:2711
static const double x3[11]
Double_t ATan(Double_t)
Definition: TMath.h:577
static constexpr double g