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] type method to compute (from 1 to 9).
1153///
1154/// #### NOTE:
1155/// When the input is not sorted, an array of integers of size n needs
1156/// to be allocated. It can be passed by the user in parameter index,
1157/// or, if not passed, it will be allocated inside the function
1158///
1159/// ### Following types are provided:
1160/// - Discontinuous:
1161/// - type=1 - inverse of the empirical distribution function
1162/// - type=2 - like type 1, but with averaging at discontinuities
1163/// - type=3 - SAS definition: nearest even order statistic
1164/// - Piecewise linear continuous:
1165/// - In this case, sample quantiles can be obtained by linear interpolation
1166/// between the k-th order statistic and p(k).
1167/// -type=4 - linear interpolation of empirical cdf, p(k)=k/n;
1168/// - type=5 - a very popular definition, p(k) = (k-0.5)/n;
1169/// - type=6 - used by Minitab and SPSS, p(k) = k/(n+1);
1170/// - type=7 - used by S-Plus and R, p(k) = (k-1)/(n-1);
1171/// - type=8 - resulting sample quantiles are approximately median unbiased
1172/// regardless of the distribution of x. p(k) = (k-1/3)/(n+1/3);
1173/// - type=9 - resulting sample quantiles are approximately unbiased, when
1174/// the sample comes from Normal distribution. p(k)=(k-3/8)/(n+1/4);
1175///
1176/// default type = 7
1177///
1178/// ### References:
1179/// 1. Hyndman, R.J and Fan, Y, (1996) "Sample quantiles in statistical packages"
1180/// American Statistician, 50, 361-365
1181/// 2. R Project documentation for the function quantile of package {stats}
1182
1183void TMath::Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted, Int_t *index, Int_t type)
1184{
1185
1186 if (type<1 || type>9){
1187 printf("illegal value of type\n");
1188 return;
1189 }
1190 Int_t *ind = 0;
1191 Bool_t isAllocated = kFALSE;
1192 if (!isSorted){
1193 if (index) ind = index;
1194 else {
1195 ind = new Int_t[n];
1196 isAllocated = kTRUE;
1197 }
1198 }
1199
1200 // re-implemented by L.M (9/11/2011) to fix bug https://savannah.cern.ch/bugs/?87251
1201 // following now exactly R implementation (available in library/stats/R/quantile.R )
1202 // which follows precisely Hyndman-Fan paper
1203 // (older implementation had a bug for type =3)
1204
1205 for (Int_t i=0; i<nprob; i++){
1206
1207 Double_t nppm = 0;
1208 Double_t gamma = 0;
1209 Int_t j = 0;
1210
1211 //Discontinuous functions
1212 // type = 1,2, or 3
1213 if (type < 4 ){
1214 if (type == 3)
1215 nppm = n*prob[i]-0.5; // use m = -0.5
1216 else
1217 nppm = n*prob[i]; // use m = 0
1218
1219 // be careful with machine precision
1220 double eps = 4 * TMath::Limits<Double_t>::Epsilon();
1221 j = TMath::FloorNint(nppm + eps);
1222
1223 // LM : fix for numerical problems if nppm is actually equal to j, but results different for numerical error
1224 // g in the paper is nppm -j
1225 if (type == 1)
1226 gamma = ( (nppm -j) > j*TMath::Limits<Double_t>::Epsilon() ) ? 1 : 0;
1227 else if (type == 2)
1228 gamma = ( (nppm -j) > j*TMath::Limits<Double_t>::Epsilon() ) ? 1 : 0.5;
1229 else if (type == 3)
1230 // gamma = 0 for g=0 and j even
1231 gamma = ( TMath::Abs(nppm -j) <= j*TMath::Limits<Double_t>::Epsilon() && TMath::Even(j) ) ? 0 : 1;
1232
1233 }
1234 else {
1235 // for continuous quantiles type 4 to 9)
1236 // define alpha and beta
1237 double a = 0;
1238 double b = 0;
1239 if (type == 4) { a = 0; b = 1; }
1240 else if (type == 5) { a = 0.5; b = 0.5; }
1241 else if (type == 6) { a = 0.; b = 0.; }
1242 else if (type == 7) { a = 1.; b = 1.; }
1243 else if (type == 8) { a = 1./3.; b = a; }
1244 else if (type == 9) { a = 3./8.; b = a; }
1245
1246 // be careful with machine precision
1247 double eps = 4 * TMath::Limits<Double_t>::Epsilon();
1248 nppm = a + prob[i] * ( n + 1 -a -b); // n * p + m
1249 j = TMath::FloorNint(nppm + eps);
1250 gamma = nppm - j;
1251 if (gamma < eps) gamma = 0;
1252 }
1253
1254 // since index j starts from 1 first is j-1 and second is j
1255 // add protection to keep indices in range [0,n-1]
1256 int first = (j > 0 && j <=n) ? j-1 : ( j <=0 ) ? 0 : n-1;
1257 int second = (j > 0 && j < n) ? j : ( j <=0 ) ? 0 : n-1;
1258
1259 Double_t xj, xjj;
1260 if (isSorted){
1261 xj = x[first];
1262 xjj = x[second];
1263 } else {
1264 xj = TMath::KOrdStat(n, x, first, ind);
1265 xjj = TMath::KOrdStat(n, x, second, ind);
1266 }
1267 quantiles[i] = (1-gamma)*xj + gamma*xjj;
1268 // printf("x[j] = %12f x[j+1] = %12f \n",xj, xjj);
1269
1270 }
1271
1272
1273
1274 if (isAllocated)
1275 delete [] ind;
1276}
1277
1278////////////////////////////////////////////////////////////////////////////////
1279/// Bubble sort variant to obtain the order of an array's elements into
1280/// an index in order to do more useful things than the standard built
1281/// in functions.
1282/// \param[in] *arr1 is unchanged;
1283/// \param[in] *arr2 is the array of indicies corresponding to the descending value
1284/// of arr1 with arr2[0] corresponding to the largest arr1 value and
1285/// arr2[Narr] the smallest.
1286///
1287/// \author Adrian Bevan (bevan@slac.stanford.edu)
1288
1289void TMath::BubbleHigh(Int_t Narr, Double_t *arr1, Int_t *arr2)
1290{
1291 if (Narr <= 0) return;
1292 double *localArr1 = new double[Narr];
1293 int *localArr2 = new int[Narr];
1294 int iEl;
1295 int iEl2;
1296
1297 for(iEl = 0; iEl < Narr; iEl++) {
1298 localArr1[iEl] = arr1[iEl];
1299 localArr2[iEl] = iEl;
1300 }
1301
1302 for (iEl = 0; iEl < Narr; iEl++) {
1303 for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1304 if (localArr1[iEl2-1] < localArr1[iEl2]) {
1305 double tmp = localArr1[iEl2-1];
1306 localArr1[iEl2-1] = localArr1[iEl2];
1307 localArr1[iEl2] = tmp;
1308
1309 int tmp2 = localArr2[iEl2-1];
1310 localArr2[iEl2-1] = localArr2[iEl2];
1311 localArr2[iEl2] = tmp2;
1312 }
1313 }
1314 }
1315
1316 for (iEl = 0; iEl < Narr; iEl++) {
1317 arr2[iEl] = localArr2[iEl];
1318 }
1319 delete [] localArr2;
1320 delete [] localArr1;
1321}
1322
1323////////////////////////////////////////////////////////////////////////////////
1324/// Opposite ordering of the array arr2[] to that of BubbleHigh.
1325///
1326/// \author Adrian Bevan (bevan@slac.stanford.edu)
1327
1328void TMath::BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
1329{
1330 if (Narr <= 0) return;
1331 double *localArr1 = new double[Narr];
1332 int *localArr2 = new int[Narr];
1333 int iEl;
1334 int iEl2;
1335
1336 for (iEl = 0; iEl < Narr; iEl++) {
1337 localArr1[iEl] = arr1[iEl];
1338 localArr2[iEl] = iEl;
1339 }
1340
1341 for (iEl = 0; iEl < Narr; iEl++) {
1342 for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1343 if (localArr1[iEl2-1] > localArr1[iEl2]) {
1344 double tmp = localArr1[iEl2-1];
1345 localArr1[iEl2-1] = localArr1[iEl2];
1346 localArr1[iEl2] = tmp;
1347
1348 int tmp2 = localArr2[iEl2-1];
1349 localArr2[iEl2-1] = localArr2[iEl2];
1350 localArr2[iEl2] = tmp2;
1351 }
1352 }
1353 }
1354
1355 for (iEl = 0; iEl < Narr; iEl++) {
1356 arr2[iEl] = localArr2[iEl];
1357 }
1358 delete [] localArr2;
1359 delete [] localArr1;
1360}
1361
1362
1363////////////////////////////////////////////////////////////////////////////////
1364/// Calculates hash index from any char string.
1365/// Based on pre-calculated table of 256 specially selected numbers.
1366/// These numbers are selected in such a way, that for string
1367/// length == 4 (integer number) the hash is unambiguous, i.e.
1368/// from hash value we can recalculate input (no degeneration).
1369///
1370/// The quality of hash method is good enough, that
1371/// "random" numbers made as R = Hash(1), Hash(2), ...Hash(N)
1372/// tested by <R>, <R*R>, <Ri*Ri+1> gives the same result
1373/// as for libc rand().
1374///
1375/// - For string: i = TMath::Hash(string,nstring);
1376/// - For int: i = TMath::Hash(&intword,sizeof(int));
1377/// - For pointer: i = TMath::Hash(&pointer,sizeof(void*));
1378///
1379/// V.Perev
1380/// This function is kept for back compatibility. The code previously in this function
1381/// has been moved to the static function TString::Hash
1382
1383ULong_t TMath::Hash(const void *txt, Int_t ntxt)
1384{
1385 return TString::Hash(txt,ntxt);
1386}
1387
1388
1389////////////////////////////////////////////////////////////////////////////////
1390
1391ULong_t TMath::Hash(const char *txt)
1392{
1393 return Hash(txt, Int_t(strlen(txt)));
1394}
1395
1396////////////////////////////////////////////////////////////////////////////////
1397/// Compute the modified Bessel function I_0(x) for any real x.
1398///
1399/// \author NvE 12-mar-2000 UU-SAP Utrecht
1400
1402{
1403 // Parameters of the polynomial approximation
1404 const Double_t p1=1.0, p2=3.5156229, p3=3.0899424,
1405 p4=1.2067492, p5=0.2659732, p6=3.60768e-2, p7=4.5813e-3;
1406
1407 const Double_t q1= 0.39894228, q2= 1.328592e-2, q3= 2.25319e-3,
1408 q4=-1.57565e-3, q5= 9.16281e-3, q6=-2.057706e-2,
1409 q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3;
1410
1411 const Double_t k1 = 3.75;
1412 Double_t ax = TMath::Abs(x);
1413
1414 Double_t y=0, result=0;
1415
1416 if (ax < k1) {
1417 Double_t xx = x/k1;
1418 y = xx*xx;
1419 result = p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))));
1420 } else {
1421 y = k1/ax;
1422 result = (TMath::Exp(ax)/TMath::Sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
1423 }
1424 return result;
1425}
1426
1427////////////////////////////////////////////////////////////////////////////////
1428/// Compute the modified Bessel function K_0(x) for positive real x.
1429///
1430/// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1431/// Applied Mathematics Series vol. 55 (1964), Washington.
1432///
1433/// \author NvE 12-mar-2000 UU-SAP Utrecht
1434
1436{
1437 // Parameters of the polynomial approximation
1438 const Double_t p1=-0.57721566, p2=0.42278420, p3=0.23069756,
1439 p4= 3.488590e-2, p5=2.62698e-3, p6=1.0750e-4, p7=7.4e-6;
1440
1441 const Double_t q1= 1.25331414, q2=-7.832358e-2, q3= 2.189568e-2,
1442 q4=-1.062446e-2, q5= 5.87872e-3, q6=-2.51540e-3, q7=5.3208e-4;
1443
1444 if (x <= 0) {
1445 Error("TMath::BesselK0", "*K0* Invalid argument x = %g\n",x);
1446 return 0;
1447 }
1448
1449 Double_t y=0, result=0;
1450
1451 if (x <= 2) {
1452 y = x*x/4;
1453 result = (-log(x/2.)*TMath::BesselI0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1454 } else {
1455 y = 2/x;
1456 result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
1457 }
1458 return result;
1459}
1460
1461////////////////////////////////////////////////////////////////////////////////
1462/// Compute the modified Bessel function I_1(x) for any real x.
1463///
1464/// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1465/// Applied Mathematics Series vol. 55 (1964), Washington.
1466///
1467/// \author NvE 12-mar-2000 UU-SAP Utrecht
1468
1470{
1471 // Parameters of the polynomial approximation
1472 const Double_t p1=0.5, p2=0.87890594, p3=0.51498869,
1473 p4=0.15084934, p5=2.658733e-2, p6=3.01532e-3, p7=3.2411e-4;
1474
1475 const Double_t q1= 0.39894228, q2=-3.988024e-2, q3=-3.62018e-3,
1476 q4= 1.63801e-3, q5=-1.031555e-2, q6= 2.282967e-2,
1477 q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3;
1478
1479 const Double_t k1 = 3.75;
1480 Double_t ax = TMath::Abs(x);
1481
1482 Double_t y=0, result=0;
1483
1484 if (ax < k1) {
1485 Double_t xx = x/k1;
1486 y = xx*xx;
1487 result = x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1488 } else {
1489 y = k1/ax;
1490 result = (exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))));
1491 if (x < 0) result = -result;
1492 }
1493 return result;
1494}
1495
1496////////////////////////////////////////////////////////////////////////////////
1497/// Compute the modified Bessel function K_1(x) for positive real x.
1498///
1499/// M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
1500/// Applied Mathematics Series vol. 55 (1964), Washington.
1501///
1502/// \author NvE 12-mar-2000 UU-SAP Utrecht
1503
1505{
1506 // Parameters of the polynomial approximation
1507 const Double_t p1= 1., p2= 0.15443144, p3=-0.67278579,
1508 p4=-0.18156897, p5=-1.919402e-2, p6=-1.10404e-3, p7=-4.686e-5;
1509
1510 const Double_t q1= 1.25331414, q2= 0.23498619, q3=-3.655620e-2,
1511 q4= 1.504268e-2, q5=-7.80353e-3, q6= 3.25614e-3, q7=-6.8245e-4;
1512
1513 if (x <= 0) {
1514 Error("TMath::BesselK1", "*K1* Invalid argument x = %g\n",x);
1515 return 0;
1516 }
1517
1518 Double_t y=0,result=0;
1519
1520 if (x <= 2) {
1521 y = x*x/4;
1522 result = (log(x/2.)*TMath::BesselI1(x))+(1./x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
1523 } else {
1524 y = 2/x;
1525 result = (exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7))))));
1526 }
1527 return result;
1528}
1529
1530////////////////////////////////////////////////////////////////////////////////
1531/// Compute the Integer Order Modified Bessel function K_n(x)
1532/// for n=0,1,2,... and positive real x.
1533///
1534/// \author NvE 12-mar-2000 UU-SAP Utrecht
1535
1537{
1538 if (x <= 0 || n < 0) {
1539 Error("TMath::BesselK", "*K* Invalid argument(s) (n,x) = (%d, %g)\n",n,x);
1540 return 0;
1541 }
1542
1543 if (n==0) return TMath::BesselK0(x);
1544 if (n==1) return TMath::BesselK1(x);
1545
1546 // Perform upward recurrence for all x
1547 Double_t tox = 2/x;
1548 Double_t bkm = TMath::BesselK0(x);
1550 Double_t bkp = 0;
1551 for (Int_t j=1; j<n; j++) {
1552 bkp = bkm+Double_t(j)*tox*bk;
1553 bkm = bk;
1554 bk = bkp;
1555 }
1556 return bk;
1557}
1558
1559////////////////////////////////////////////////////////////////////////////////
1560/// Compute the Integer Order Modified Bessel function I_n(x)
1561/// for n=0,1,2,... and any real x.
1562///
1563/// \author NvE 12-mar-2000 UU-SAP Utrecht
1564
1566{
1567 Int_t iacc = 40; // Increase to enhance accuracy
1568 const Double_t kBigPositive = 1.e10;
1569 const Double_t kBigNegative = 1.e-10;
1570
1571 if (n < 0) {
1572 Error("TMath::BesselI", "*I* Invalid argument (n,x) = (%d, %g)\n",n,x);
1573 return 0;
1574 }
1575
1576 if (n==0) return TMath::BesselI0(x);
1577 if (n==1) return TMath::BesselI1(x);
1578
1579 if (x == 0) return 0;
1580 if (TMath::Abs(x) > kBigPositive) return 0;
1581
1582 Double_t tox = 2/TMath::Abs(x);
1583 Double_t bip = 0, bim = 0;
1584 Double_t bi = 1;
1585 Double_t result = 0;
1586 Int_t m = 2*((n+Int_t(sqrt(Float_t(iacc*n)))));
1587 for (Int_t j=m; j>=1; j--) {
1588 bim = bip+Double_t(j)*tox*bi;
1589 bip = bi;
1590 bi = bim;
1591 // Renormalise to prevent overflows
1592 if (TMath::Abs(bi) > kBigPositive) {
1593 result *= kBigNegative;
1594 bi *= kBigNegative;
1595 bip *= kBigNegative;
1596 }
1597 if (j==n) result=bip;
1598 }
1599
1600 result *= TMath::BesselI0(x)/bi; // Normalise with BesselI0(x)
1601 if ((x < 0) && (n%2 == 1)) result = -result;
1602
1603 return result;
1604}
1605
1606////////////////////////////////////////////////////////////////////////////////
1607/// Returns the Bessel function J0(x) for any real x.
1608
1610{
1611 Double_t ax,z;
1612 Double_t xx,y,result,result1,result2;
1613 const Double_t p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
1614 const Double_t p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
1615 const Double_t p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
1616 const Double_t p10 = 59272.64853, p11 = 267.8532712;
1617
1618 const Double_t q1 = 0.785398164;
1619 const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1620 const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1621 const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1622 const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1623 const Double_t q10 = 0.934935152e-7, q11 = 0.636619772;
1624
1625 if ((ax=fabs(x)) < 8) {
1626 y=x*x;
1627 result1 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6))));
1628 result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1629 result = result1/result2;
1630 } else {
1631 z = 8/ax;
1632 y = z*z;
1633 xx = ax-q1;
1634 result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1635 result2 = q6 + y*(q7 + y*(q8 + y*(q9 - y*q10)));
1636 result = sqrt(q11/ax)*(cos(xx)*result1-z*sin(xx)*result2);
1637 }
1638 return result;
1639}
1640
1641////////////////////////////////////////////////////////////////////////////////
1642/// Returns the Bessel function J1(x) for any real x.
1643
1645{
1646 Double_t ax,z;
1647 Double_t xx,y,result,result1,result2;
1648 const Double_t p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
1649 const Double_t p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
1650 const Double_t p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
1651 const Double_t p10 = 99447.43394, p11 = 376.9991397;
1652
1653 const Double_t q1 = 2.356194491;
1654 const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1655 const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1656 const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1657 const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1658 const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1659
1660 if ((ax=fabs(x)) < 8) {
1661 y=x*x;
1662 result1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6)))));
1663 result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1664 result = result1/result2;
1665 } else {
1666 z = 8/ax;
1667 y = z*z;
1668 xx = ax-q1;
1669 result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1670 result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1671 result = sqrt(q11/ax)*(cos(xx)*result1-z*sin(xx)*result2);
1672 if (x < 0) result = -result;
1673 }
1674 return result;
1675}
1676
1677////////////////////////////////////////////////////////////////////////////////
1678/// Returns the Bessel function Y0(x) for positive x.
1679
1681{
1682 Double_t z,xx,y,result,result1,result2;
1683 const Double_t p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
1684 const Double_t p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
1685 const Double_t p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
1686 const Double_t p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
1687
1688 const Double_t q1 = 0.785398164;
1689 const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1690 const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1691 const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1692 const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1693 const Double_t q10 = -0.934945152e-7, q11 = 0.636619772;
1694
1695 if (x < 8) {
1696 y = x*x;
1697 result1 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6))));
1698 result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y))));
1699 result = (result1/result2) + p12*TMath::BesselJ0(x)*log(x);
1700 } else {
1701 z = 8/x;
1702 y = z*z;
1703 xx = x-q1;
1704 result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1705 result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1706 result = sqrt(q11/x)*(sin(xx)*result1+z*cos(xx)*result2);
1707 }
1708 return result;
1709}
1710
1711////////////////////////////////////////////////////////////////////////////////
1712/// Returns the Bessel function Y1(x) for positive x.
1713
1715{
1716 Double_t z,xx,y,result,result1,result2;
1717 const Double_t p1 = -0.4900604943e13, p2 = 0.1275274390e13;
1718 const Double_t p3 = -0.5153438139e11, p4 = 0.7349264551e9;
1719 const Double_t p5 = -0.4237922726e7, p6 = 0.8511937935e4;
1720 const Double_t p7 = 0.2499580570e14, p8 = 0.4244419664e12;
1721 const Double_t p9 = 0.3733650367e10, p10 = 0.2245904002e8;
1722 const Double_t p11 = 0.1020426050e6, p12 = 0.3549632885e3;
1723 const Double_t p13 = 0.636619772;
1724 const Double_t q1 = 2.356194491;
1725 const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1726 const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1727 const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1728 const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1729 const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1730
1731 if (x < 8) {
1732 y=x*x;
1733 result1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*p6)))));
1734 result2 = p7 + y*(p8 + y*(p9 + y*(p10 + y*(p11 + y*(p12+y)))));
1735 result = (result1/result2) + p13*(TMath::BesselJ1(x)*log(x)-1/x);
1736 } else {
1737 z = 8/x;
1738 y = z*z;
1739 xx = x-q1;
1740 result1 = 1 + y*(q2 + y*(q3 + y*(q4 + y*q5)));
1741 result2 = q6 + y*(q7 + y*(q8 + y*(q9 + y*q10)));
1742 result = sqrt(q11/x)*(sin(xx)*result1+z*cos(xx)*result2);
1743 }
1744 return result;
1745}
1746
1747////////////////////////////////////////////////////////////////////////////////
1748/// Struve Functions of Order 0
1749///
1750/// Converted from CERNLIB M342 by Rene Brun.
1751
1753{
1754 const Int_t n1 = 15;
1755 const Int_t n2 = 25;
1756 const Double_t c1[16] = { 1.00215845609911981, -1.63969292681309147,
1757 1.50236939618292819, -.72485115302121872,
1758 .18955327371093136, -.03067052022988,
1759 .00337561447375194, -2.6965014312602e-4,
1760 1.637461692612e-5, -7.8244408508e-7,
1761 3.021593188e-8, -9.6326645e-10,
1762 2.579337e-11, -5.8854e-13,
1763 1.158e-14, -2e-16 };
1764 const Double_t c2[26] = { .99283727576423943, -.00696891281138625,
1765 1.8205103787037e-4, -1.063258252844e-5,
1766 9.8198294287e-7, -1.2250645445e-7,
1767 1.894083312e-8, -3.44358226e-9,
1768 7.1119102e-10, -1.6288744e-10,
1769 4.065681e-11, -1.091505e-11,
1770 3.12005e-12, -9.4202e-13,
1771 2.9848e-13, -9.872e-14,
1772 3.394e-14, -1.208e-14,
1773 4.44e-15, -1.68e-15,
1774 6.5e-16, -2.6e-16,
1775 1.1e-16, -4e-17,
1776 2e-17, -1e-17 };
1777
1778 const Double_t c0 = 2/TMath::Pi();
1779
1780 Int_t i;
1781 Double_t alfa, h, r, y, b0, b1, b2;
1783
1784 v = TMath::Abs(x);
1785 if (v < 8) {
1786 y = v/8;
1787 h = 2*y*y -1;
1788 alfa = h + h;
1789 b0 = 0;
1790 b1 = 0;
1791 b2 = 0;
1792 for (i = n1; i >= 0; --i) {
1793 b0 = c1[i] + alfa*b1 - b2;
1794 b2 = b1;
1795 b1 = b0;
1796 }
1797 h = y*(b0 - h*b2);
1798 } else {
1799 r = 1/v;
1800 h = 128*r*r -1;
1801 alfa = h + h;
1802 b0 = 0;
1803 b1 = 0;
1804 b2 = 0;
1805 for (i = n2; i >= 0; --i) {
1806 b0 = c2[i] + alfa*b1 - b2;
1807 b2 = b1;
1808 b1 = b0;
1809 }
1810 h = TMath::BesselY0(v) + r*c0*(b0 - h*b2);
1811 }
1812 if (x < 0) h = -h;
1813 return h;
1814}
1815
1816////////////////////////////////////////////////////////////////////////////////
1817/// Struve Functions of Order 1
1818///
1819/// Converted from CERNLIB M342 by Rene Brun.
1820
1822{
1823 const Int_t n3 = 16;
1824 const Int_t n4 = 22;
1825 const Double_t c3[17] = { .5578891446481605, -.11188325726569816,
1826 -.16337958125200939, .32256932072405902,
1827 -.14581632367244242, .03292677399374035,
1828 -.00460372142093573, 4.434706163314e-4,
1829 -3.142099529341e-5, 1.7123719938e-6,
1830 -7.416987005e-8, 2.61837671e-9,
1831 -7.685839e-11, 1.9067e-12,
1832 -4.052e-14, 7.5e-16,
1833 -1e-17 };
1834 const Double_t c4[23] = { 1.00757647293865641, .00750316051248257,
1835 -7.043933264519e-5, 2.66205393382e-6,
1836 -1.8841157753e-7, 1.949014958e-8,
1837 -2.6126199e-9, 4.236269e-10,
1838 -7.955156e-11, 1.679973e-11,
1839 -3.9072e-12, 9.8543e-13,
1840 -2.6636e-13, 7.645e-14,
1841 -2.313e-14, 7.33e-15,
1842 -2.42e-15, 8.3e-16,
1843 -3e-16, 1.1e-16,
1844 -4e-17, 2e-17,-1e-17 };
1845
1846 const Double_t c0 = 2/TMath::Pi();
1847 const Double_t cc = 2/(3*TMath::Pi());
1848
1849 Int_t i, i1;
1850 Double_t alfa, h, r, y, b0, b1, b2;
1852
1853 if (v == 0) {
1854 h = 0;
1855 } else if (v <= 0.3) {
1856 y = v*v;
1857 r = 1;
1858 h = 1;
1859 i1 = (Int_t)(-8. / TMath::Log10(v));
1860 for (i = 1; i <= i1; ++i) {
1861 h = -h*y / ((2*i+ 1)*(2*i + 3));
1862 r += h;
1863 }
1864 h = cc*y*r;
1865 } else if (v < 8) {
1866 h = v*v/32 -1;
1867 alfa = h + h;
1868 b0 = 0;
1869 b1 = 0;
1870 b2 = 0;
1871 for (i = n3; i >= 0; --i) {
1872 b0 = c3[i] + alfa*b1 - b2;
1873 b2 = b1;
1874 b1 = b0;
1875 }
1876 h = b0 - h*b2;
1877 } else {
1878 h = 128/(v*v) -1;
1879 alfa = h + h;
1880 b0 = 0;
1881 b1 = 0;
1882 b2 = 0;
1883 for (i = n4; i >= 0; --i) {
1884 b0 = c4[i] + alfa*b1 - b2;
1885 b2 = b1;
1886 b1 = b0;
1887 }
1888 h = TMath::BesselY1(v) + c0*(b0 - h*b2);
1889 }
1890 return h;
1891}
1892
1893////////////////////////////////////////////////////////////////////////////////
1894/// Modified Struve Function of Order 0.
1895/// By Kirill Filimonov.
1896
1898{
1899 const Double_t pi=TMath::Pi();
1900
1901 Double_t s=1.0;
1902 Double_t r=1.0;
1903
1904 Double_t a0,sl0,a1,bi0;
1905
1906 Int_t km,i;
1907
1908 if (x<=20.) {
1909 a0=2.0*x/pi;
1910 for (i=1; i<=60;i++) {
1911 r*=(x/(2*i+1))*(x/(2*i+1));
1912 s+=r;
1913 if(TMath::Abs(r/s)<1.e-12) break;
1914 }
1915 sl0=a0*s;
1916 } else {
1917 km=int(5*(x+1.0));
1918 if(x>=50.0)km=25;
1919 for (i=1; i<=km; i++) {
1920 r*=(2*i-1)*(2*i-1)/x/x;
1921 s+=r;
1922 if(TMath::Abs(r/s)<1.0e-12) break;
1923 }
1924 a1=TMath::Exp(x)/TMath::Sqrt(2*pi*x);
1925 r=1.0;
1926 bi0=1.0;
1927 for (i=1; i<=16; i++) {
1928 r=0.125*r*(2.0*i-1.0)*(2.0*i-1.0)/(i*x);
1929 bi0+=r;
1930 if(TMath::Abs(r/bi0)<1.0e-12) break;
1931 }
1932
1933 bi0=a1*bi0;
1934 sl0=-2.0/(pi*x)*s+bi0;
1935 }
1936 return sl0;
1937}
1938
1939////////////////////////////////////////////////////////////////////////////////
1940/// Modified Struve Function of Order 1.
1941/// By Kirill Filimonov.
1942
1944{
1945 const Double_t pi=TMath::Pi();
1946 Double_t a1,sl1,bi1,s;
1947 Double_t r=1.0;
1948 Int_t km,i;
1949
1950 if (x<=20.) {
1951 s=0.0;
1952 for (i=1; i<=60;i++) {
1953 r*=x*x/(4.0*i*i-1.0);
1954 s+=r;
1955 if(TMath::Abs(r)<TMath::Abs(s)*1.e-12) break;
1956 }
1957 sl1=2.0/pi*s;
1958 } else {
1959 s=1.0;
1960 km=int(0.5*x);
1961 if(x>50.0)km=25;
1962 for (i=1; i<=km; i++) {
1963 r*=(2*i+3)*(2*i+1)/x/x;
1964 s+=r;
1965 if(TMath::Abs(r/s)<1.0e-12) break;
1966 }
1967 sl1=2.0/pi*(-1.0+1.0/(x*x)+3.0*s/(x*x*x*x));
1968 a1=TMath::Exp(x)/TMath::Sqrt(2*pi*x);
1969 r=1.0;
1970 bi1=1.0;
1971 for (i=1; i<=16; i++) {
1972 r=-0.125*r*(4.0-(2.0*i-1.0)*(2.0*i-1.0))/(i*x);
1973 bi1+=r;
1974 if(TMath::Abs(r/bi1)<1.0e-12) break;
1975 }
1976 sl1+=a1*bi1;
1977 }
1978 return sl1;
1979}
1980
1981////////////////////////////////////////////////////////////////////////////////
1982/// Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
1983
1985{
1987}
1988
1989////////////////////////////////////////////////////////////////////////////////
1990/// Continued fraction evaluation by modified Lentz's method
1991/// used in calculation of incomplete Beta function.
1992
1994{
1995 Int_t itmax = 500;
1996 Double_t eps = 3.e-14;
1997 Double_t fpmin = 1.e-30;
1998
1999 Int_t m, m2;
2000 Double_t aa, c, d, del, qab, qam, qap;
2001 Double_t h;
2002 qab = a+b;
2003 qap = a+1.0;
2004 qam = a-1.0;
2005 c = 1.0;
2006 d = 1.0 - qab*x/qap;
2007 if (TMath::Abs(d)<fpmin) d=fpmin;
2008 d=1.0/d;
2009 h=d;
2010 for (m=1; m<=itmax; m++) {
2011 m2=m*2;
2012 aa = m*(b-m)*x/((qam+ m2)*(a+m2));
2013 d = 1.0 +aa*d;
2014 if(TMath::Abs(d)<fpmin) d = fpmin;
2015 c = 1 +aa/c;
2016 if (TMath::Abs(c)<fpmin) c = fpmin;
2017 d=1.0/d;
2018 h*=d*c;
2019 aa = -(a+m)*(qab +m)*x/((a+m2)*(qap+m2));
2020 d=1.0+aa*d;
2021 if(TMath::Abs(d)<fpmin) d = fpmin;
2022 c = 1.0 +aa/c;
2023 if (TMath::Abs(c)<fpmin) c = fpmin;
2024 d=1.0/d;
2025 del = d*c;
2026 h*=del;
2027 if (TMath::Abs(del-1)<=eps) break;
2028 }
2029 if (m>itmax) {
2030 Info("TMath::BetaCf", "a or b too big, or itmax too small, a=%g, b=%g, x=%g, h=%g, itmax=%d",
2031 a,b,x,h,itmax);
2032 }
2033 return h;
2034}
2035
2036////////////////////////////////////////////////////////////////////////////////
2037/// Computes the probability density function of the Beta distribution
2038/// (the distribution function is computed in BetaDistI).
2039/// The first argument is the point, where the function will be
2040/// computed, second and third are the function parameters.
2041/// Since the Beta distribution is bounded on both sides, it's often
2042/// used to represent processes with natural lower and upper limits.
2043
2045{
2046 if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2047 Error("TMath::BetaDist", "parameter value outside allowed range");
2048 return 0;
2049 }
2050 Double_t beta = TMath::Beta(p, q);
2051 Double_t r = TMath::Power(x, p-1)*TMath::Power(1-x, q-1)/beta;
2052 return r;
2053}
2054
2055////////////////////////////////////////////////////////////////////////////////
2056/// Computes the distribution function of the Beta distribution.
2057/// The first argument is the point, where the function will be
2058/// computed, second and third are the function parameters.
2059/// Since the Beta distribution is bounded on both sides, it's often
2060/// used to represent processes with natural lower and upper limits.
2061
2063{
2064 if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2065 Error("TMath::BetaDistI", "parameter value outside allowed range");
2066 return 0;
2067 }
2068 Double_t betai = TMath::BetaIncomplete(x, p, q);
2069 return betai;
2070}
2071
2072////////////////////////////////////////////////////////////////////////////////
2073/// Calculates the incomplete Beta-function.
2074
2076{
2078}
2079
2080////////////////////////////////////////////////////////////////////////////////
2081/// Calculate the binomial coefficient n over k.
2082
2084{
2085 if (n<0 || k<0 || n<k) return TMath::SignalingNaN();
2086 if (k==0 || n==k) return 1;
2087
2088 Int_t k1=TMath::Min(k,n-k);
2089 Int_t k2=n-k1;
2090 Double_t fact=k2+1;
2091 for (Double_t i=k1;i>1.;--i)
2092 fact *= (k2+i)/i;
2093 return fact;
2094}
2095
2096////////////////////////////////////////////////////////////////////////////////
2097/// Suppose an event occurs with probability _p_ per trial
2098/// Then the probability P of its occurring _k_ or more times
2099/// in _n_ trials is termed a cumulative binomial probability
2100/// the formula is P = sum_from_j=k_to_n(TMath::Binomial(n, j)*
2101/// *TMath::Power(p, j)*TMath::Power(1-p, n-j)
2102/// For _n_ larger than 12 BetaIncomplete is a much better way
2103/// to evaluate the sum than would be the straightforward sum calculation
2104/// for _n_ smaller than 12 either method is acceptable
2105/// ("Numerical Recipes")
2106/// --implementation by Anna Kreshuk
2107
2109{
2110 if(k <= 0) return 1.0;
2111 if(k > n) return 0.0;
2112 if(k == n) return TMath::Power(p, n);
2113
2114 return BetaIncomplete(p, Double_t(k), Double_t(n-k+1));
2115}
2116
2117////////////////////////////////////////////////////////////////////////////////
2118/// Computes the density of Cauchy distribution at point x
2119/// by default, standard Cauchy distribution is used (t=0, s=1)
2120/// - t is the location parameter
2121/// - s is the scale parameter
2122///
2123/// The Cauchy distribution, also called Lorentzian distribution,
2124/// is a continuous distribution describing resonance behavior
2125/// The mean and standard deviation of the Cauchy distribution are undefined.
2126/// The practical meaning of this is that collecting 1,000 data points gives
2127/// no more accurate an estimate of the mean and standard deviation than
2128/// does a single point.
2129/// The formula was taken from "Engineering Statistics Handbook" on site
2130/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3663.htm
2131/// Implementation by Anna Kreshuk.
2132///
2133/// Example:
2134///
2135/// ~~~ {cpp}
2136/// TF1* fc = new TF1("fc", "TMath::CauchyDist(x, [0], [1])", -5, 5);
2137/// fc->SetParameters(0, 1);
2138/// fc->Draw();
2139/// ~~~
2140
2142{
2143 Double_t temp= (x-t)*(x-t)/(s*s);
2144 Double_t result = 1/(s*TMath::Pi()*(1+temp));
2145 return result;
2146}
2147
2148////////////////////////////////////////////////////////////////////////////////
2149/// Evaluate the quantiles of the chi-squared probability distribution function.
2150/// Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
2151/// implemented by Anna Kreshuk.
2152/// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
2153///
2154/// \param[in] p the probability value, at which the quantile is computed
2155/// \param[in] ndf number of degrees of freedom
2156
2158{
2159 Double_t c[]={0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
2160 4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
2161 84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
2162 210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
2163 462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
2164 932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
2165 2520.0, 5040.0};
2166 Double_t e = 5e-7;
2167 Double_t aa = 0.6931471806;
2168 Int_t maxit = 20;
2169 Double_t ch, p1, p2, q, t, a, b, x;
2170 Double_t s1, s2, s3, s4, s5, s6;
2171
2172 if (ndf <= 0) return 0;
2173
2174 Double_t g = TMath::LnGamma(0.5*ndf);
2175
2176 Double_t xx = 0.5 * ndf;
2177 Double_t cp = xx - 1;
2178 if (ndf >= TMath::Log(p)*(-c[5])){
2179 //starting approximation for ndf less than or equal to 0.32
2180 if (ndf > c[3]) {
2182 //starting approximation using Wilson and Hilferty estimate
2183 p1=c[2]/ndf;
2184 ch = ndf*TMath::Power((x*TMath::Sqrt(p1) + 1 - p1), 3);
2185 if (ch > c[6]*ndf + 6)
2186 ch = -2 * (TMath::Log(1-p) - cp * TMath::Log(0.5 * ch) + g);
2187 } else {
2188 ch = c[4];
2189 a = TMath::Log(1-p);
2190 do{
2191 q = ch;
2192 p1 = 1 + ch * (c[7]+ch);
2193 p2 = ch * (c[9] + ch * (c[8] + ch));
2194 t = -0.5 + (c[7] + 2 * ch) / p1 - (c[9] + ch * (c[10] + 3 * ch)) / p2;
2195 ch = ch - (1 - TMath::Exp(a + g + 0.5 * ch + cp * aa) *p2 / p1) / t;
2196 }while (TMath::Abs(q/ch - 1) > c[1]);
2197 }
2198 } else {
2199 ch = TMath::Power((p * xx * TMath::Exp(g + xx * aa)),(1./xx));
2200 if (ch < e) return ch;
2201 }
2202//call to algorithm AS 239 and calculation of seven term Taylor series
2203 for (Int_t i=0; i<maxit; i++){
2204 q = ch;
2205 p1 = 0.5 * ch;
2206 p2 = p - TMath::Gamma(xx, p1);
2207
2208 t = p2 * TMath::Exp(xx * aa + g + p1 - cp * TMath::Log(ch));
2209 b = t / ch;
2210 a = 0.5 * t - b * cp;
2211 s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] +c[11] * a))))) / c[24];
2212 s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
2213 s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
2214 s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
2215 s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
2216 s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
2217 ch = ch + t * (1 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
2218 if (TMath::Abs(q / ch - 1) > e) break;
2219 }
2220 return ch;
2221}
2222
2223////////////////////////////////////////////////////////////////////////////////
2224/// Computes the density function of F-distribution
2225/// (probability function, integral of density, is computed in FDistI).
2226///
2227/// Parameters N and M stand for degrees of freedom of chi-squares
2228/// mentioned above parameter F is the actual variable x of the
2229/// density function p(x) and the point at which the density function
2230/// is calculated.
2231///
2232/// ### About F distribution:
2233/// F-distribution arises in testing whether two random samples
2234/// have the same variance. It is the ratio of two chi-square
2235/// distributions, with N and M degrees of freedom respectively,
2236/// where each chi-square is first divided by it's number of degrees
2237/// of freedom.
2238/// Implementation by Anna Kreshuk.
2239
2241{
2243}
2244
2245////////////////////////////////////////////////////////////////////////////////
2246/// Calculates the cumulative distribution function of F-distribution,
2247/// this function occurs in the statistical test of whether two observed
2248/// samples have the same variance. For this test a certain statistic F,
2249/// the ratio of observed dispersion of the first sample to that of the
2250/// second sample, is calculated. N and M stand for numbers of degrees
2251/// of freedom in the samples 1-FDistI() is the significance level at
2252/// which the hypothesis "1 has smaller variance than 2" can be rejected.
2253/// A small numerical value of 1 - FDistI() implies a very significant
2254/// rejection, in turn implying high confidence in the hypothesis
2255/// "1 has variance greater than 2".
2256///
2257/// Implementation by Anna Kreshuk.
2258
2260{
2262 return fi;
2263}
2264
2265////////////////////////////////////////////////////////////////////////////////
2266/// Computes the density function of Gamma distribution at point x.
2267///
2268/// \param[in] gamma shape parameter
2269/// \param[in] mu location parameter
2270/// \param[in] beta scale parameter
2271///
2272/// The definition can be found in "Engineering Statistics Handbook" on site
2273/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
2274/// use now implementation in ROOT::Math::gamma_pdf
2275///
2276/// Begin_Macro
2277/// {
2278/// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2279///
2280/// c1->SetLogy();
2281/// c1->SetGridx();
2282/// c1->SetGridy();
2283///
2284/// TF1 *gdist = new TF1("gdist", "TMath::GammaDist(x, [0], [1], [2])", 0, 10);
2285///
2286/// gdist->SetParameters(0.5, 0., 1.);
2287/// gdist->SetLineColor(2);
2288/// TF1 *gdist1 = gdist->DrawCopy("L");
2289/// gdist->SetParameters(1.0, 0., 1.);
2290/// gdist->SetLineColor(3);
2291/// TF1 *gdist2 = gdist->DrawCopy("LSAME");
2292/// gdist->SetParameters(2.0, 0., 1.);
2293/// gdist->SetLineColor(4);
2294/// TF1 *gdist3 = gdist->DrawCopy("LSAME");
2295/// gdist->SetParameters(5.0, 0., 1.);
2296/// gdist->SetLineColor(6);
2297/// TF1 *gdist4 = gdist->DrawCopy("LSAME");
2298///
2299/// legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2300/// legend->AddEntry(gdist1, "gamma = 0.5, mu = 0, beta = 1", "L");
2301/// legend->AddEntry(gdist2, "gamma = 1.0, mu = 0, beta = 1", "L");
2302/// legend->AddEntry(gdist3, "gamma = 2.0, mu = 0, beta = 1", "L");
2303/// legend->AddEntry(gdist4, "gamma = 5.0, mu = 0, beta = 1", "L");
2304/// legend->Draw();
2305/// }
2306/// End_Macro
2307
2309{
2310 if ((x<mu) || (gamma<=0) || (beta <=0)) {
2311 Error("TMath::GammaDist", "illegal parameter values x = %f , gamma = %f beta = %f",x,gamma,beta);
2312 return 0;
2313 }
2315}
2316
2317////////////////////////////////////////////////////////////////////////////////
2318/// Computes the probability density function of Laplace distribution
2319/// at point x, with location parameter alpha and shape parameter beta.
2320/// By default, alpha=0, beta=1
2321/// This distribution is known under different names, most common is
2322/// double exponential distribution, but it also appears as
2323/// the two-tailed exponential or the bilateral exponential distribution
2324
2326{
2327 Double_t temp;
2328 temp = TMath::Exp(-TMath::Abs((x-alpha)/beta));
2329 temp /= (2.*beta);
2330 return temp;
2331}
2332
2333////////////////////////////////////////////////////////////////////////////////
2334/// Computes the distribution function of Laplace distribution
2335/// at point x, with location parameter alpha and shape parameter beta.
2336/// By default, alpha=0, beta=1
2337/// This distribution is known under different names, most common is
2338/// double exponential distribution, but it also appears as
2339/// the two-tailed exponential or the bilateral exponential distribution
2340
2342{
2343 Double_t temp;
2344 if (x<=alpha){
2345 temp = 0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2346 } else {
2347 temp = 1-0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2348 }
2349 return temp;
2350}
2351
2352////////////////////////////////////////////////////////////////////////////////
2353/// Computes the density of LogNormal distribution at point x.
2354/// Variable X has lognormal distribution if Y=Ln(X) has normal distribution
2355///
2356/// \param[in] sigma is the shape parameter
2357/// \param[in] theta is the location parameter
2358/// \param[in] m is the scale parameter
2359///
2360/// The formula was taken from "Engineering Statistics Handbook" on site
2361/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm
2362/// Implementation using ROOT::Math::lognormal_pdf
2363///
2364/// Begin_Macro
2365/// {
2366/// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2367///
2368/// c1->SetLogy();
2369/// c1->SetGridx();
2370/// c1->SetGridy();
2371///
2372/// TF1 *logn = new TF1("logn", "TMath::LogNormal(x, [0], [1], [2])", 0, 5);
2373/// logn->SetMinimum(1e-3);
2374///
2375/// logn->SetParameters(0.5, 0., 1.);
2376/// logn->SetLineColor(2);
2377/// TF1 *logn1 = logn->DrawCopy("L");
2378/// logn->SetParameters(1.0, 0., 1.);
2379/// logn->SetLineColor(3);
2380/// TF1 *logn2 = logn->DrawCopy("LSAME");
2381/// logn->SetParameters(2.0, 0., 1.);
2382/// logn->SetLineColor(4);
2383/// TF1 *logn3 = logn->DrawCopy("LSAME");
2384/// logn->SetParameters(5.0, 0., 1.);
2385/// logn->SetLineColor(6);
2386/// TF1 *logn4 = logn->DrawCopy("LSAME");
2387///
2388/// legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2389/// legend->AddEntry(logn1, "sigma = 0.5, theta = 0, m = 1", "L");
2390/// legend->AddEntry(logn2, "sigma = 1.0, theta = 0, m = 1", "L");
2391/// legend->AddEntry(logn3, "sigma = 2.0, theta = 0, m = 1", "L");
2392/// legend->AddEntry(logn4, "sigma = 5.0, theta = 0, m = 1", "L");
2393/// legend->Draw();
2394/// }
2395/// End_Macro
2396
2398{
2399 if ((x<theta) || (sigma<=0) || (m<=0)) {
2400 Error("TMath::Lognormal", "illegal parameter values");
2401 return 0;
2402 }
2403 // lognormal_pdf uses same definition of http://en.wikipedia.org/wiki/Log-normal_distribution
2404 // where mu = log(m)
2405
2407
2408}
2409
2410////////////////////////////////////////////////////////////////////////////////
2411/// Computes quantiles for standard normal distribution N(0, 1)
2412/// at probability p
2413///
2414/// ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.
2415
2417{
2418 if ((p<=0)||(p>=1)) {
2419 Error("TMath::NormQuantile", "probability outside (0, 1)");
2420 return 0;
2421 }
2422
2423 Double_t a0 = 3.3871328727963666080e0;
2424 Double_t a1 = 1.3314166789178437745e+2;
2425 Double_t a2 = 1.9715909503065514427e+3;
2426 Double_t a3 = 1.3731693765509461125e+4;
2427 Double_t a4 = 4.5921953931549871457e+4;
2428 Double_t a5 = 6.7265770927008700853e+4;
2429 Double_t a6 = 3.3430575583588128105e+4;
2430 Double_t a7 = 2.5090809287301226727e+3;
2431 Double_t b1 = 4.2313330701600911252e+1;
2432 Double_t b2 = 6.8718700749205790830e+2;
2433 Double_t b3 = 5.3941960214247511077e+3;
2434 Double_t b4 = 2.1213794301586595867e+4;
2435 Double_t b5 = 3.9307895800092710610e+4;
2436 Double_t b6 = 2.8729085735721942674e+4;
2437 Double_t b7 = 5.2264952788528545610e+3;
2438 Double_t c0 = 1.42343711074968357734e0;
2439 Double_t c1 = 4.63033784615654529590e0;
2440 Double_t c2 = 5.76949722146069140550e0;
2441 Double_t c3 = 3.64784832476320460504e0;
2442 Double_t c4 = 1.27045825245236838258e0;
2443 Double_t c5 = 2.41780725177450611770e-1;
2444 Double_t c6 = 2.27238449892691845833e-2;
2445 Double_t c7 = 7.74545014278341407640e-4;
2446 Double_t d1 = 2.05319162663775882187e0;
2447 Double_t d2 = 1.67638483018380384940e0;
2448 Double_t d3 = 6.89767334985100004550e-1;
2449 Double_t d4 = 1.48103976427480074590e-1;
2450 Double_t d5 = 1.51986665636164571966e-2;
2451 Double_t d6 = 5.47593808499534494600e-4;
2452 Double_t d7 = 1.05075007164441684324e-9;
2453 Double_t e0 = 6.65790464350110377720e0;
2454 Double_t e1 = 5.46378491116411436990e0;
2455 Double_t e2 = 1.78482653991729133580e0;
2456 Double_t e3 = 2.96560571828504891230e-1;
2457 Double_t e4 = 2.65321895265761230930e-2;
2458 Double_t e5 = 1.24266094738807843860e-3;
2459 Double_t e6 = 2.71155556874348757815e-5;
2460 Double_t e7 = 2.01033439929228813265e-7;
2461 Double_t f1 = 5.99832206555887937690e-1;
2462 Double_t f2 = 1.36929880922735805310e-1;
2463 Double_t f3 = 1.48753612908506148525e-2;
2464 Double_t f4 = 7.86869131145613259100e-4;
2465 Double_t f5 = 1.84631831751005468180e-5;
2466 Double_t f6 = 1.42151175831644588870e-7;
2467 Double_t f7 = 2.04426310338993978564e-15;
2468
2469 Double_t split1 = 0.425;
2470 Double_t split2=5.;
2471 Double_t konst1=0.180625;
2472 Double_t konst2=1.6;
2473
2474 Double_t q, r, quantile;
2475 q=p-0.5;
2476 if (TMath::Abs(q)<split1) {
2477 r=konst1-q*q;
2478 quantile = q* (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
2479 * r + a2) * r + a1) * r + a0) /
2480 (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
2481 * r + b2) * r + b1) * r + 1.);
2482 } else {
2483 if(q<0) r=p;
2484 else r=1-p;
2485 //error case
2486 if (r<=0)
2487 quantile=0;
2488 else {
2490 if (r<=split2) {
2491 r=r-konst2;
2492 quantile=(((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
2493 * r + c2) * r + c1) * r + c0) /
2494 (((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
2495 * r + d2) * r + d1) * r + 1);
2496 } else{
2497 r=r-split2;
2498 quantile=(((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
2499 * r + e2) * r + e1) * r + e0) /
2500 (((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
2501 * r + f2) * r + f1) * r + 1);
2502 }
2503 if (q<0) quantile=-quantile;
2504 }
2505 }
2506 return quantile;
2507}
2508
2509////////////////////////////////////////////////////////////////////////////////
2510/// Simple recursive algorithm to find the permutations of
2511/// n natural numbers, not necessarily all distinct
2512/// adapted from CERNLIB routine PERMU.
2513/// The input array has to be initialised with a non descending
2514/// sequence. The method returns kFALSE when
2515/// all combinations are exhausted.
2516
2518{
2519 Int_t i,itmp;
2520 Int_t i1=-1;
2521
2522 // find rightmost upward transition
2523 for(i=n-2; i>-1; i--) {
2524 if(a[i]<a[i+1]) {
2525 i1=i;
2526 break;
2527 }
2528 }
2529 // no more upward transitions, end of the story
2530 if(i1==-1) return kFALSE;
2531 else {
2532 // find lower right element higher than the lower
2533 // element of the upward transition
2534 for(i=n-1;i>i1;i--) {
2535 if(a[i] > a[i1]) {
2536 // swap the two
2537 itmp=a[i1];
2538 a[i1]=a[i];
2539 a[i]=itmp;
2540 break;
2541 }
2542 }
2543 // order the rest, in fact just invert, as there
2544 // are only downward transitions from here on
2545 for(i=0;i<(n-i1-1)/2;i++) {
2546 itmp=a[i1+i+1];
2547 a[i1+i+1]=a[n-i-1];
2548 a[n-i-1]=itmp;
2549 }
2550 }
2551 return kTRUE;
2552}
2553
2554////////////////////////////////////////////////////////////////////////////////
2555/// Computes density function for Student's t- distribution
2556/// (the probability function (integral of density) is computed in StudentI).
2557///
2558/// First parameter stands for x - the actual variable of the
2559/// density function p(x) and the point at which the density is calculated.
2560/// Second parameter stands for number of degrees of freedom.
2561///
2562/// About Student distribution:
2563/// Student's t-distribution is used for many significance tests, for example,
2564/// for the Student's t-tests for the statistical significance of difference
2565/// between two sample means and for confidence intervals for the difference
2566/// between two population means.
2567///
2568/// Example: suppose we have a random sample of size n drawn from normal
2569/// distribution with mean Mu and st.deviation Sigma. Then the variable
2570///
2571/// t = (sample_mean - Mu)/(sample_deviation / sqrt(n))
2572///
2573/// has Student's t-distribution with n-1 degrees of freedom.
2574///
2575/// NOTE that this function's second argument is number of degrees of freedom,
2576/// not the sample size.
2577///
2578/// As the number of degrees of freedom grows, t-distribution approaches
2579/// Normal(0,1) distribution.
2580///
2581/// Implementation by Anna Kreshuk.
2582
2584{
2585 if (ndf < 1) {
2586 return 0;
2587 }
2588
2589 Double_t r = ndf;
2590 Double_t rh = 0.5*r;
2591 Double_t rh1 = rh + 0.5;
2592 Double_t denom = TMath::Sqrt(r*TMath::Pi())*TMath::Gamma(rh)*TMath::Power(1+T*T/r, rh1);
2593 return TMath::Gamma(rh1)/denom;
2594}
2595
2596////////////////////////////////////////////////////////////////////////////////
2597/// Calculates the cumulative distribution function of Student's
2598/// t-distribution second parameter stands for number of degrees of freedom,
2599/// not for the number of samples
2600/// if x has Student's t-distribution, the function returns the probability of
2601/// x being less than T.
2602///
2603/// Implementation by Anna Kreshuk.
2604
2606{
2607 Double_t r = ndf;
2608
2609 Double_t si = (T>0) ?
2610 (1 - 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5)) :
2611 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5);
2612 return si;
2613}
2614
2615////////////////////////////////////////////////////////////////////////////////
2616/// Computes quantiles of the Student's t-distribution
2617/// 1st argument is the probability, at which the quantile is computed
2618/// 2nd argument - the number of degrees of freedom of the
2619/// Student distribution
2620/// When the 3rd argument lower_tail is kTRUE (default)-
2621/// the algorithm returns such x0, that
2622///
2623/// P(x < x0)=p
2624///
2625/// upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that
2626///
2627/// P(x > x0)=p
2628///
2629/// the algorithm was taken from:
2630/// G.W.Hill, "Algorithm 396, Student's t-quantiles"
2631/// "Communications of the ACM", 13(10), October 1970
2632
2634{
2635 Double_t quantile;
2636 Double_t temp;
2637 Bool_t neg;
2638 Double_t q;
2639 if (ndf<1 || p>=1 || p<=0) {
2640 Error("TMath::StudentQuantile", "illegal parameter values");
2641 return 0;
2642 }
2643 if ((lower_tail && p>0.5)||(!lower_tail && p<0.5)){
2644 neg=kFALSE;
2645 q=2*(lower_tail ? (1-p) : p);
2646 } else {
2647 neg=kTRUE;
2648 q=2*(lower_tail? p : (1-p));
2649 }
2650
2651 if ((ndf-1)<1e-8) {
2652 temp=TMath::PiOver2()*q;
2653 quantile = TMath::Cos(temp)/TMath::Sin(temp);
2654 } else {
2655 if ((ndf-2)<1e-8) {
2656 quantile = TMath::Sqrt(2./(q*(2-q))-2);
2657 } else {
2658 Double_t a=1./(ndf-0.5);
2659 Double_t b=48./(a*a);
2660 Double_t c=((20700*a/b -98)*a-16)*a+96.36;
2661 Double_t d=((94.5/(b+c)-3.)/b+1)*TMath::Sqrt(a*TMath::PiOver2())*ndf;
2662 Double_t x=q*d;
2663 Double_t y=TMath::Power(x, (2./ndf));
2664 if (y>0.05+a){
2665 //asymptotic inverse expansion about normal
2666 x=TMath::NormQuantile(q*0.5);
2667 y=x*x;
2668 if (ndf<5) c+=0.3*(ndf-4.5)*(x+0.6);
2669 c+=(((0.05*d*x-5.)*x-7.)*x-2.)*x +b;
2670 y=(((((0.4*y+6.3)*y+36.)*y+94.5)/c - y-3.)/b+1)*x;
2671 y=a*y*y;
2672 if(y>0.002) y = TMath::Exp(y)-1;
2673 else y += 0.5*y*y;
2674 } else {
2675 y=((1./(((ndf+6.)/(ndf*y)-0.089*d-0.822)*(ndf+2.)*3)+0.5/(ndf+4.))*y-1.)*
2676 (ndf+1.)/(ndf+2.)+1/y;
2677 }
2678 quantile = TMath::Sqrt(ndf*y);
2679 }
2680 }
2681 if(neg) quantile=-quantile;
2682 return quantile;
2683}
2684
2685////////////////////////////////////////////////////////////////////////////////
2686/// Returns the value of the Vavilov density function
2687///
2688/// \param[in] x the point were the density function is evaluated
2689/// \param[in] kappa value of kappa (distribution parameter)
2690/// \param[in] beta2 value of beta2 (distribution parameter)
2691///
2692/// The algorithm was taken from the CernLib function vavden(G115)
2693/// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2694/// Nucl.Instr. and Meth. B47(1990), 215-224
2695///
2696/// Accuracy: quote from the reference above:
2697///
2698/// "The results of our code have been compared with the values of the Vavilov
2699/// density function computed numerically in an accurate way: our approximation
2700/// shows a difference of less than 3% around the peak of the density function, slowly
2701/// increasing going towards the extreme tails to the right and to the left"
2702///
2703/// Begin_Macro
2704/// {
2705/// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2706///
2707/// c1->SetGridx();
2708/// c1->SetGridy();
2709///
2710/// TF1 *vavilov = new TF1("vavilov", "TMath::Vavilov(x, [0], [1])", -3, 11);
2711///
2712/// vavilov->SetParameters(0.5, 0.);
2713/// vavilov->SetLineColor(2);
2714/// TF1 *vavilov1 = vavilov->DrawCopy("L");
2715/// vavilov->SetParameters(0.3, 0.);
2716/// vavilov->SetLineColor(3);
2717/// TF1 *vavilov2 = vavilov->DrawCopy("LSAME");
2718/// vavilov->SetParameters(0.2, 0.);
2719/// vavilov->SetLineColor(4);
2720/// TF1 *vavilov3 = vavilov->DrawCopy("LSAME");
2721/// vavilov->SetParameters(0.1, 0.);
2722/// vavilov->SetLineColor(6);
2723/// TF1 *vavilov4 = vavilov->DrawCopy("LSAME");
2724///
2725/// legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2726/// legend->AddEntry(vavilov1, "kappa = 0.5, beta2 = 0", "L");
2727/// legend->AddEntry(vavilov2, "kappa = 0.3, beta2 = 0", "L");
2728/// legend->AddEntry(vavilov3, "kappa = 0.2, beta2 = 0", "L");
2729/// legend->AddEntry(vavilov4, "kappa = 0.1, beta2 = 0", "L");
2730/// legend->Draw();
2731/// }
2732/// End_Macro
2733
2735{
2736 Double_t *ac = new Double_t[14];
2737 Double_t *hc = new Double_t[9];
2738
2739 Int_t itype;
2740 Int_t npt;
2741 TMath::VavilovSet(kappa, beta2, 0, 0, ac, hc, itype, npt);
2742 Double_t v = TMath::VavilovDenEval(x, ac, hc, itype);
2743 delete [] ac;
2744 delete [] hc;
2745 return v;
2746}
2747
2748////////////////////////////////////////////////////////////////////////////////
2749/// Returns the value of the Vavilov distribution function
2750///
2751/// \param[in] x the point were the density function is evaluated
2752/// \param[in] kappa value of kappa (distribution parameter)
2753/// \param[in] beta2 value of beta2 (distribution parameter)
2754///
2755/// The algorithm was taken from the CernLib function vavden(G115)
2756///
2757/// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2758/// Nucl.Instr. and Meth. B47(1990), 215-224
2759///
2760/// Accuracy: quote from the reference above:
2761///
2762/// "The results of our code have been compared with the values of the Vavilov
2763/// density function computed numerically in an accurate way: our approximation
2764/// shows a difference of less than 3% around the peak of the density function, slowly
2765/// increasing going towards the extreme tails to the right and to the left"
2766
2768{
2769 Double_t *ac = new Double_t[14];
2770 Double_t *hc = new Double_t[9];
2771 Double_t *wcm = new Double_t[201];
2772 Int_t itype;
2773 Int_t npt;
2774 Int_t k;
2775 Double_t xx, v;
2776 TMath::VavilovSet(kappa, beta2, 1, wcm, ac, hc, itype, npt);
2777 if (x < ac[0]) v = 0;
2778 else if (x >=ac[8]) v = 1;
2779 else {
2780 xx = x - ac[0];
2781 k = Int_t(xx*ac[10]);
2782 v = TMath::Min(wcm[k] + (xx - k*ac[9])*(wcm[k+1]-wcm[k])*ac[10], 1.);
2783 }
2784 delete [] ac;
2785 delete [] hc;
2786 delete [] wcm;
2787 return v;
2788}
2789
2790////////////////////////////////////////////////////////////////////////////////
2791/// Returns the value of the Landau distribution function at point x.
2792/// The algorithm was taken from the Cernlib function dislan(G110)
2793/// Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
2794/// distribution", Computer Phys.Comm., 31(1984), 97-111
2795
2797{
2799}
2800
2801
2802////////////////////////////////////////////////////////////////////////////////
2803/// Internal function, called by Vavilov and VavilovI
2804
2805void TMath::VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt)
2806{
2807
2808 Double_t BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
2809 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
2810 BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
2811
2812 Double_t FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
2813 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
2814 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
2815
2816 Double_t FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
2817
2818 Double_t EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
2819 0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
2820
2821 Double_t U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
2822 0. , 0.24880692e-1, 0.47404356e-2,
2823 -0.74445130e-3, 0.73225731e-2, 0. ,
2824 0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
2825
2826 Double_t U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
2827 0. , 0.42127077e-1, 0.73167928e-2,
2828 -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
2829 0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
2830
2831 Double_t U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
2832 0. , 0.23834176e-1, 0.21624675e-2,
2833 -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
2834 0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
2835
2836 Double_t U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
2837 -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
2838 0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
2839 -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
2840
2841 Double_t U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
2842 0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
2843 0.48736023e-3, 0.34850854e-2, 0. ,
2844 -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
2845
2846 Double_t U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
2847 -0.30188807e-2, -0.84479939e-3, 0. ,
2848 0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
2849 -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
2850
2851 Double_t U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
2852 -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
2853 0. , 0.50505298e+0};
2854 Double_t U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
2855 -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
2856 -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
2857 0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
2858
2859 Double_t V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
2860 0. , 0.45091424e-1, 0.80559636e-2,
2861 -0.38974523e-2, 0. , -0.30634124e-2,
2862 0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
2863
2864 Double_t V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
2865 0. , 0.12693873e+0, 0.22999801e-1,
2866 -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
2867 0. , 0.19716857e-1, 0.32596742e-2};
2868
2869 Double_t V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
2870 -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
2871 -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
2872 0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
2873
2874 Double_t V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
2875 0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
2876 0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
2877 -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
2878
2879 Double_t V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
2880 -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
2881 0.56856517e-2, -0.13438433e-1, 0. ,
2882 0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
2883
2884 Double_t V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
2885 0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
2886 0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
2887 -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
2888
2889 Double_t V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
2890 0. , -0.12218009e+1, -0.32324120e+0,
2891 -0.27373591e-1, 0.12173464e+0, 0. ,
2892 0. , 0.40917471e-1};
2893
2894 Double_t V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
2895 0. , -0.18160479e+1, -0.50919193e+0,
2896 -0.51384654e-1, 0.21413992e+0, 0. ,
2897 0. , 0.66596366e-1};
2898
2899 Double_t W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
2900 -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
2901 -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
2902 0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
2903
2904 Double_t W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
2905 -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
2906 0. , 0.23020158e-1, 0.50574313e-2,
2907 0.94550140e-2, 0.19300232e-1};
2908
2909 Double_t W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
2910 -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
2911 -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
2912 0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
2913
2914 Double_t W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
2915 0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
2916 0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
2917 -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
2918
2919 Double_t W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
2920 0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
2921 0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
2922 -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
2923
2924 Double_t W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
2925 0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
2926 0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
2927 -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
2928
2929 Double_t W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
2930 0. , -0.14540925e+1, -0.39529833e+0,
2931 -0.44293243e-1, 0.88741049e-1};
2932
2933 itype = 0;
2934 if (rkappa <0.01 || rkappa >12) {
2935 Error("Vavilov distribution", "illegal value of kappa");
2936 return;
2937 }
2938
2939 Double_t DRK[6];
2940 Double_t DSIGM[6];
2941 Double_t ALFA[8];
2942 Int_t j;
2943 Double_t x, y, xx, yy, x2, x3, y2, y3, xy, p2, p3, q2, q3, pq;
2944 if (rkappa >= 0.29) {
2945 itype = 1;
2946 npt = 100;
2947 Double_t wk = 1./TMath::Sqrt(rkappa);
2948
2949 AC[0] = (-0.032227*beta2-0.074275)*rkappa + (0.24533*beta2+0.070152)*wk + (-0.55610*beta2-3.1579);
2950 AC[8] = (-0.013483*beta2-0.048801)*rkappa + (-1.6921*beta2+8.3656)*wk + (-0.73275*beta2-3.5226);
2951 DRK[1] = wk*wk;
2952 DSIGM[1] = TMath::Sqrt(rkappa/(1-0.5*beta2));
2953 for (j=1; j<=4; j++) {
2954 DRK[j+1] = DRK[1]*DRK[j];
2955 DSIGM[j+1] = DSIGM[1]*DSIGM[j];
2956 ALFA[j+1] = (FNINV[j]-beta2*FNINV[j+1])*DRK[j];
2957 }
2958 HC[0]=TMath::Log(rkappa)+beta2+0.42278434;
2959 HC[1]=DSIGM[1];
2960 HC[2]=ALFA[3]*DSIGM[3];
2961 HC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
2962 HC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*HC[2];
2963 HC[5]=HC[2]*HC[2];
2964 HC[6]=HC[2]*HC[3];
2965 HC[7]=HC[2]*HC[5];
2966 for (j=2; j<=7; j++)
2967 HC[j]*=EDGEC[j];
2968 HC[8]=0.39894228*HC[1];
2969 }
2970 else if (rkappa >=0.22) {
2971 itype = 2;
2972 npt = 150;
2973 x = 1+(rkappa-BKMXX3)*FBKX3;
2974 y = 1+(TMath::Sqrt(beta2)-BKMXY3)*FBKY3;
2975 xx = 2*x;
2976 yy = 2*y;
2977 x2 = xx*x-1;
2978 x3 = xx*x2-x;
2979 y2 = yy*y-1;
2980 y3 = yy*y2-y;
2981 xy = x*y;
2982 p2 = x2*y;
2983 p3 = x3*y;
2984 q2 = y2*x;
2985 q3 = y3*x;
2986 pq = x2*y2;
2987 AC[1] = W1[1] + W1[2]*x + W1[4]*x3 + W1[5]*y + W1[6]*y2 + W1[7]*y3 +
2988 W1[8]*xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
2989 AC[2] = W2[1] + W2[2]*x + W2[3]*x2 + W2[4]*x3 + W2[5]*y + W2[6]*y2 +
2990 W2[8]*xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
2991 AC[3] = W3[1] + W3[3]*x2 + W3[4]*x3 + W3[5]*y + W3[6]*y2 + W3[7]*y3 +
2992 W3[8]*xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
2993 AC[4] = W4[1] + W4[2]*x + W4[3]*x2 + W4[4]*x3 + W4[5]*y + W4[6]*y2 + W4[7]*y3 +
2994 W4[8]*xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
2995 AC[5] = W5[1] + W5[2]*x + W5[4]*x3 + W5[5]*y + W5[6]*y2 + W5[7]*y3 +
2996 W5[8]*xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
2997 AC[6] = W6[1] + W6[2]*x + W6[3]*x2 + W6[4]*x3 + W6[5]*y + W6[6]*y2 + W6[7]*y3 +
2998 W6[8]*xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
2999 AC[8] = W8[1] + W8[2]*x + W8[3]*x2 + W8[5]*y + W8[6]*y2 + W8[7]*y3 + W8[8]*xy;
3000 AC[0] = -3.05;
3001 } else if (rkappa >= 0.12) {
3002 itype = 3;
3003 npt = 200;
3004 x = 1 + (rkappa-BKMXX2)*FBKX2;
3005 y = 1 + (TMath::Sqrt(beta2)-BKMXY2)*FBKY2;
3006 xx = 2*x;
3007 yy = 2*y;
3008 x2 = xx*x-1;
3009 x3 = xx*x2-x;
3010 y2 = yy*y-1;
3011 y3 = yy*y2-y;
3012 xy = x*y;
3013 p2 = x2*y;
3014 p3 = x3*y;
3015 q2 = y2*x;
3016 q3 = y3*x;
3017 pq = x2*y2;
3018 AC[1] = V1[1] + V1[2]*x + V1[3]*x2 + V1[5]*y + V1[6]*y2 + V1[7]*y3 +
3019 V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
3020 AC[2] = V2[1] + V2[2]*x + V2[3]*x2 + V2[5]*y + V2[6]*y2 + V2[7]*y3 +
3021 V2[8]*xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
3022 AC[3] = V3[1] + V3[2]*x + V3[3]*x2 + V3[4]*x3 + V3[5]*y + V3[6]*y2 + V3[7]*y3 +
3023 V3[8]*xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
3024 AC[4] = V4[1] + V4[2]*x + V4[3]*x2 + V4[4]*x3 + V4[5]*y + V4[6]*y2 + V4[7]*y3 +
3025 V4[8]*xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
3026 AC[5] = V5[1] + V5[2]*x + V5[3]*x2 + V5[4]*x3 + V5[5]*y + V5[6]*y2 + V5[7]*y3 +
3027 V5[8]*xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
3028 AC[6] = V6[1] + V6[2]*x + V6[3]*x2 + V6[4]*x3 + V6[5]*y + V6[6]*y2 + V6[7]*y3 +
3029 V6[8]*xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
3030 AC[7] = V7[1] + V7[2]*x + V7[3]*x2 + V7[5]*y + V7[6]*y2 + V7[7]*y3 +
3031 V7[8]*xy + V7[11]*q2;
3032 AC[8] = V8[1] + V8[2]*x + V8[3]*x2 + V8[5]*y + V8[6]*y2 + V8[7]*y3 +
3033 V8[8]*xy + V8[11]*q2;
3034 AC[0] = -3.04;
3035 } else {
3036 itype = 4;
3037 if (rkappa >=0.02) itype = 3;
3038 npt = 200;
3039 x = 1+(rkappa-BKMXX1)*FBKX1;
3040 y = 1+(TMath::Sqrt(beta2)-BKMXY1)*FBKY1;
3041 xx = 2*x;
3042 yy = 2*y;
3043 x2 = xx*x-1;
3044 x3 = xx*x2-x;
3045 y2 = yy*y-1;
3046 y3 = yy*y2-y;
3047 xy = x*y;
3048 p2 = x2*y;
3049 p3 = x3*y;
3050 q2 = y2*x;
3051 q3 = y3*x;
3052 pq = x2*y2;
3053 if (itype==3){
3054 AC[1] = U1[1] + U1[2]*x + U1[3]*x2 + U1[5]*y + U1[6]*y2 + U1[7]*y3 +
3055 U1[8]*xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
3056 AC[2] = U2[1] + U2[2]*x + U2[3]*x2 + U2[5]*y + U2[6]*y2 + U2[7]*y3 +
3057 U2[8]*xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
3058 AC[3] = U3[1] + U3[2]*x + U3[3]*x2 + U3[5]*y + U3[6]*y2 + U3[7]*y3 +
3059 U3[8]*xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
3060 AC[4] = U4[1] + U4[2]*x + U4[3]*x2 + U4[4]*x3 + U4[5]*y + U4[6]*y2 + U4[7]*y3 +
3061 U4[8]*xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
3062 AC[5] = U5[1] + U5[2]*x + U5[3]*x2 + U5[4]*x3 + U5[5]*y + U5[6]*y2 + U5[7]*y3 +
3063 U5[8]*xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
3064 AC[6] = U6[1] + U6[2]*x + U6[3]*x2 + U6[4]*x3 + U6[5]*y + U6[7]*y3 +
3065 U6[8]*xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
3066 AC[7] = U7[1] + U7[2]*x + U7[3]*x2 + U7[4]*x3 + U7[5]*y + U7[6]*y2 + U7[8]*xy;
3067 }
3068 AC[8] = U8[1] + U8[2]*x + U8[3]*x2 + U8[4]*x3 + U8[5]*y + U8[6]*y2 + U8[7]*y3 +
3069 U8[8]*xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
3070 AC[0] = -3.03;
3071 }
3072
3073 AC[9] = (AC[8] - AC[0])/npt;
3074 AC[10] = 1./AC[9];
3075 if (itype == 3) {
3076 x = (AC[7]-AC[8])/(AC[7]*AC[8]);
3077 y = 1./TMath::Log(AC[8]/AC[7]);
3078 p2 = AC[7]*AC[7];
3079 AC[11] = p2*(AC[1]*TMath::Exp(-AC[2]*(AC[7]+AC[5]*p2)-
3080 AC[3]*TMath::Exp(-AC[4]*(AC[7]+AC[6]*p2)))-0.045*y/AC[7])/(1+x*y*AC[7]);
3081 AC[12] = (0.045+x*AC[11])*y;
3082 }
3083 if (itype == 4) AC[13] = 0.995/LandauI(AC[8]);
3084
3085 if (mode==0) return;
3086 //
3087 x = AC[0];
3088 WCM[0] = 0;
3089 Double_t fl, fu;
3090 Int_t k;
3091 fl = TMath::VavilovDenEval(x, AC, HC, itype);
3092 for (k=1; k<=npt; k++) {
3093 x += AC[9];
3094 fu = TMath::VavilovDenEval(x, AC, HC, itype);
3095 WCM[k] = WCM[k-1] + fl + fu;
3096 fl = fu;
3097 }
3098 x = 0.5*AC[9];
3099 for (k=1; k<=npt; k++)
3100 WCM[k]*=x;
3101}
3102
3103////////////////////////////////////////////////////////////////////////////////
3104/// Internal function, called by Vavilov and VavilovSet
3105
3107{
3108 Double_t v = 0;
3109 if (rlam < AC[0] || rlam > AC[8])
3110 return 0;
3111 Int_t k;
3112 Double_t x, fn, s;
3113 Double_t h[10];
3114 if (itype ==1 ) {
3115 fn = 1;
3116 x = (rlam + HC[0])*HC[1];
3117 h[1] = x;
3118 h[2] = x*x -1;
3119 for (k=2; k<=8; k++) {
3120 fn++;
3121 h[k+1] = x*h[k]-fn*h[k-1];
3122 }
3123 s = 1 + HC[7]*h[9];
3124 for (k=2; k<=6; k++)
3125 s+=HC[k]*h[k+1];
3126 v = HC[8]*TMath::Exp(-0.5*x*x)*TMath::Max(s, 0.);
3127 }
3128 else if (itype == 2) {
3129 x = rlam*rlam;
3130 v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x) - AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3131 }
3132 else if (itype == 3) {
3133 if (rlam < AC[7]) {
3134 x = rlam*rlam;
3135 v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x)-AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3136 } else {
3137 x = 1./rlam;
3138 v = (AC[11]*x + AC[12])*x;
3139 }
3140 }
3141 else if (itype == 4) {
3142 v = AC[13]*TMath::Landau(rlam);
3143 }
3144 return v;
3145}
3146
3147
3148//explicitly instantiate template functions from VecCore
3149#ifdef R__HAS_VECCORE
3150#include <Math/Types.h>
3156// missing in veccore
3157// template ROOT::Double_v vecCore::math::ACos(const ROOT::Double_v & x);
3158// template ROOT::Double_v vecCore::math::Sinh(const ROOT::Double_v & x);
3159// template ROOT::Double_v vecCore::math::Cosh(const ROOT::Double_v & x);
3160// template ROOT::Double_v vecCore::math::Tanh(const ROOT::Double_v & x);
3161// template ROOT::Double_v vecCore::math::Cbrt(const ROOT::Double_v & x);
3162#endif
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
double cos(double)
double sqrt(double)
double sin(double)
double exp(double)
double log(double)
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< 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:2259
Double_t LogNormal(Double_t x, Double_t sigma, Double_t theta=0, Double_t m=1)
Computes the density of LogNormal distribution at point x.
Definition: TMath.cxx:2397
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:2044
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:3106
Double_t ACos(Double_t)
Definition: TMath.h:669
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:1565
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:1333
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:713
Double_t BinomialI(Double_t p, Int_t n, Int_t k)
Suppose an event occurs with probability p per trial Then the probability P of its occurring k or mor...
Definition: TMath.cxx:2108
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2)
Returns the value of the Vavilov density function.
Definition: TMath.cxx:2734
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2083
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:663
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:1504
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:1289
Double_t Exp(Double_t x)
Definition: TMath.h:727
Double_t BesselI1(Double_t x)
modified Bessel function K_0(x)
Definition: TMath.cxx:1469
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:2517
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:901
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:2141
Double_t ATan(Double_t)
Definition: TMath.h:675
Double_t StruveL1(Double_t x)
Modified Struve functions of order 0.
Definition: TMath.cxx:1943
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:2341
ULong_t Hash(const void *txt, Int_t ntxt)
Calculates hash index from any char string.
Definition: TMath.cxx:1383
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculate a Breit Wigner function with mean and gamma.
Definition: TMath.cxx:437
T1 Sign(T1 a, T2 b)
Definition: TMathBase.h:165
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
Definition: TMath.cxx:469
Double_t Voigt(Double_t x, Double_t sigma, Double_t lg, Int_t r=4)
Computation of Voigt function (normalised).
Definition: TMath.cxx: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:2583
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:2062
Int_t FloorNint(Double_t x)
Definition: TMath.h:707
Double_t ACosH(Double_t)
Definition: TMath.cxx:77
Double_t BesselK0(Double_t x)
modified Bessel function I_0(x)
Definition: TMath.cxx:1435
Double_t BesselY0(Double_t x)
Bessel function J1(x) for any real x.
Definition: TMath.cxx:1680
Double_t ATan2(Double_t y, Double_t x)
Definition: TMath.h:679
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:1993
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:2325
Double_t Log(Double_t x)
Definition: TMath.h:760
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:2767
Double_t Beta(Double_t p, Double_t q)
Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
Definition: TMath.cxx:1984
Double_t 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:691
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:735
Double_t BesselJ0(Double_t x)
modified Bessel function K_1(x)
Definition: TMath.cxx:1609
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:348
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Double_t StruveL0(Double_t x)
Struve functions of order 1.
Definition: TMath.cxx:1897
Double_t NormQuantile(Double_t p)
Computes quantiles for standard normal distribution N(0, 1) at probability p.
Definition: TMath.cxx:2416
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:643
constexpr Double_t Pi()
Definition: TMath.h:37
Double_t StruveH0(Double_t x)
Bessel function Y1(x) for positive x.
Definition: TMath.cxx:1752
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:2157
Double_t Sin(Double_t)
Definition: TMath.h:639
Double_t FDist(Double_t F, Double_t N, Double_t M)
Computes the density function of F-distribution (probability function, integral of density,...
Definition: TMath.cxx:2240
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition: TMath.h:908
void BubbleLow(Int_t Narr, Double_t *arr1, Int_t *arr2)
Opposite ordering of the array arr2[] to that of BubbleHigh.
Definition: TMath.cxx:1328
Double_t BesselK(Int_t n, Double_t x)
integer order modified Bessel function I_n(x)
Definition: TMath.cxx:1536
Double_t BesselJ1(Double_t x)
Bessel function J0(x) for any real x.
Definition: TMath.cxx:1644
Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b)
Calculates the incomplete Beta-function.
Definition: TMath.cxx:2075
Double_t StruveH1(Double_t x)
Struve functions of order 0.
Definition: TMath.cxx:1821
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:2796
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:1183
Double_t BesselI0(Double_t x)
integer order modified Bessel function K_n(x)
Definition: TMath.cxx:1401
void VavilovSet(Double_t rkappa, Double_t beta2, Bool_t mode, Double_t *WCM, Double_t *AC, Double_t *HC, Int_t &itype, Int_t &npt)
Internal function, called by Vavilov and VavilovI.
Definition: TMath.cxx:2805
Double_t Log10(Double_t x)
Definition: TMath.h:764
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:2605
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:2633
Double_t BesselY1(Double_t x)
Bessel function Y0(x) for positive x.
Definition: TMath.cxx:1714
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:2308
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:935
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345