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