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