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