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