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