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