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