Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
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
25#include <Math/SpecFuncMathCore.h>
26#include <Math/PdfFuncMathCore.h>
27#include <Math/ProbFuncMathCore.h>
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)
36NamespaceImp(TMath)
37#endif
38
39namespace TMath {
40
41 Double_t GamCf(Double_t a,Double_t x);
42 Double_t GamSer(Double_t a,Double_t x);
43 Double_t VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t itype);
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
51Long_t TMath::Hypot(Long_t x, Long_t y)
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
59Double_t TMath::Hypot(Double_t x, Double_t y)
60{
61 return hypot(x, y);
62}
63
64////////////////////////////////////////////////////////////////////////////////
65/// Returns the area hyperbolic sine of `x`.
66
67Double_t TMath::ASinH(Double_t x)
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
81Double_t TMath::ACosH(Double_t x)
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
95Double_t TMath::ATanH(Double_t x)
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
107Double_t TMath::Log2(Double_t x)
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
116Double_t TMath::DiLog(Double_t x)
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
190Double_t TMath::Erf(Double_t x)
191{
192 return ::ROOT::Math::erf(x);
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
199Double_t TMath::Erfc(Double_t x)
200{
201 return ::ROOT::Math::erfc(x);
202}
203
204////////////////////////////////////////////////////////////////////////////////
205/// Returns the inverse error function.
206/// x must be <-1<x<1
207
208Double_t TMath::ErfInverse(Double_t x)
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
242Double_t TMath::ErfcInverse(Double_t x)
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
252Double_t TMath::Factorial(Int_t n)
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
270Double_t TMath::Freq(Double_t x)
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
353Double_t TMath::Gamma(Double_t z)
354{
355 return ::ROOT::Math::tgamma(z);
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
369Double_t TMath::Gamma(Double_t a,Double_t x)
370{
371 return ::ROOT::Math::inc_gamma(a, x);
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
380Double_t TMath::GamCf(Double_t a,Double_t x)
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
417Double_t TMath::GamSer(Double_t a,Double_t x)
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
442Double_t TMath::BreitWigner(Double_t x, Double_t mean, Double_t gamma)
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
452Double_t TMath::BreitWignerRelativistic(Double_t x, Double_t median, Double_t gamma)
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
471Double_t TMath::Gaus(Double_t x, Double_t mean, Double_t sigma, Bool_t norm)
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
492Double_t TMath::Landau(Double_t x, Double_t mu, Double_t sigma, Bool_t norm)
493{
494 if (sigma <= 0) return 0;
495 Double_t den = ::ROOT::Math::landau_pdf( (x-mu)/sigma );
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
509Double_t TMath::LnGamma(Double_t z)
510{
511 return ::ROOT::Math::lgamma(z);
512}
513
514////////////////////////////////////////////////////////////////////////////////
515/// Normalize a vector v in place.
516/// Returns the norm of the original vector.
517
518Float_t TMath::Normalize(Float_t v[3])
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
535Double_t TMath::Normalize(Double_t v[3])
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
587Double_t TMath::Poisson(Double_t x, Double_t par)
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
615Double_t TMath::PoissonI(Double_t x, Double_t par)
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
637Double_t TMath::Prob(Double_t chi2,Int_t ndf)
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
646 return ::ROOT::Math::chisquared_cdf_c(chi2,ndf);
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
679Double_t TMath::KolmogorovProb(Double_t z)
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
805Double_t TMath::KolmogorovTest(Int_t na, const Double_t *a, Int_t nb, const Double_t *b, Option_t *option)
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 algorithm 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 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
898Double_t TMath::Voigt(Double_t xx, Double_t sigma, Double_t lg, Int_t r)
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
1107Bool_t TMath::RootsCubic(const Double_t coef[4],Double_t &a, Double_t &b, Double_t &c)
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 = nullptr;
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
1426Double_t TMath::BesselI0(Double_t x)
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
1460Double_t TMath::BesselK0(Double_t x)
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
1494Double_t TMath::BesselI1(Double_t x)
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
1529Double_t TMath::BesselK1(Double_t x)
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
1561Double_t TMath::BesselK(Int_t n,Double_t x)
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);
1574 Double_t bk = TMath::BesselK1(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
1590Double_t TMath::BesselI(Int_t n,Double_t x)
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
1634Double_t TMath::BesselJ0(Double_t x)
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
1669Double_t TMath::BesselJ1(Double_t x)
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
1705Double_t TMath::BesselY0(Double_t x)
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
1739Double_t TMath::BesselY1(Double_t x)
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
1777Double_t TMath::StruveH0(Double_t x)
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;
1807 Double_t v = TMath::Abs(x);
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
1846Double_t TMath::StruveH1(Double_t x)
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;
1876 Double_t v = TMath::Abs(x);
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
1923Double_t TMath::StruveL0(Double_t x)
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
1970Double_t TMath::StruveL1(Double_t x)
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
2011Double_t TMath::Beta(Double_t p, Double_t q)
2012{
2013 return ::ROOT::Math::beta(p, q);
2014}
2015
2016////////////////////////////////////////////////////////////////////////////////
2017/// Continued fraction evaluation by modified Lentz's method
2018/// used in calculation of incomplete Beta function.
2019
2020Double_t TMath::BetaCf(Double_t x, Double_t a, Double_t b)
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 cumulative 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
2071Double_t TMath::BetaDist(Double_t x, Double_t p, Double_t q)
2072{
2073 if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2074 Error("TMath::BetaDist", "parameter value outside allowed range");
2075 return 0;
2076 }
2077 Double_t beta = TMath::Beta(p, q);
2078 Double_t r = TMath::Power(x, p-1)*TMath::Power(1-x, q-1)/beta;
2079 return r;
2080}
2081
2082////////////////////////////////////////////////////////////////////////////////
2083/// Computes the cumulative distribution function of the Beta distribution,
2084/// i.e. the lower tail integral of TMath::BetaDist
2085/// The first argument is the point, where the function will be
2086/// computed, second and third are the function parameters.
2087/// Since the Beta distribution is bounded on both sides, it's often
2088/// used to represent processes with natural lower and upper limits.
2089
2090Double_t TMath::BetaDistI(Double_t x, Double_t p, Double_t q)
2091{
2092 if ((x<0) || (x>1) || (p<=0) || (q<=0)){
2093 Error("TMath::BetaDistI", "parameter value outside allowed range");
2094 return 0;
2095 }
2096 Double_t betai = TMath::BetaIncomplete(x, p, q);
2097 return betai;
2098}
2099
2100////////////////////////////////////////////////////////////////////////////////
2101/// Calculates the incomplete Beta-function.
2102
2103Double_t TMath::BetaIncomplete(Double_t x, Double_t a, Double_t b)
2104{
2105 return ::ROOT::Math::inc_beta(x, a, b);
2106}
2107
2108////////////////////////////////////////////////////////////////////////////////
2109/// Calculates the binomial coefficient n over k.
2110
2111Double_t TMath::Binomial(Int_t n,Int_t k)
2112{
2113 if (n<0 || k<0 || n<k) return TMath::SignalingNaN();
2114 if (k==0 || n==k) return 1;
2115
2116 Int_t k1=TMath::Min(k,n-k);
2117 Int_t k2=n-k1;
2118 Double_t fact=k2+1;
2119 for (Double_t i=k1;i>1.;--i)
2120 fact *= (k2+i)/i;
2121 return fact;
2122}
2123
2124////////////////////////////////////////////////////////////////////////////////
2125/// Suppose an event occurs with probability _p_ per trial
2126/// Then the probability P of its occurring _k_ or more times
2127/// in _n_ trials is termed a cumulative binomial probability
2128/// the formula is:
2129/// ~~~ {cpp}
2130/// P = sum_from_j=k_to_n(TMath::Binomial (n, j)**TMath::Power (p, j)*TMath::Power (1-p, n-j)
2131/// ~~~
2132/// For _n_ larger than 12 BetaIncomplete is a much better way
2133/// to evaluate the sum than would be the straightforward sum calculation
2134/// for _n_ smaller than 12 either method is acceptable ("Numerical Recipes")
2135///
2136/// Note this function is not exactly implementing the cumulative or the complement of the cumulative of the
2137/// Binomial distrinution. It is equivalent to ROOT::Math::binomial_cdf_c(k-1,p,n)
2138///
2139/// \author Anna Kreshuk
2140
2141Double_t TMath::BinomialI(Double_t p, Int_t n, Int_t k)
2142{
2143 if(k <= 0) return 1.0;
2144 if(k > n) return 0.0;
2145 if(k == n) return TMath::Power(p, n);
2146
2147 return BetaIncomplete(p, Double_t(k), Double_t(n-k+1));
2148}
2149
2150////////////////////////////////////////////////////////////////////////////////
2151/// Computes the density of Cauchy distribution at point x
2152/// by default, standard Cauchy distribution is used (t=0, s=1)
2153/// - t is the location parameter
2154/// - s is the scale parameter
2155///
2156/// The Cauchy distribution, also called Lorentzian distribution,
2157/// is a continuous distribution describing resonance behavior
2158/// The mean and standard deviation of the Cauchy distribution are undefined.
2159/// The practical meaning of this is that collecting 1,000 data points gives
2160/// no more accurate an estimate of the mean and standard deviation than
2161/// does a single point.
2162/// The formula was taken from "Engineering Statistics Handbook" on site
2163/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3663.htm
2164///
2165/// Example:
2166///
2167/// ~~~ {cpp}
2168/// TF1* fc = new TF1("fc", "TMath::CauchyDist(x, [0], [1])", -5, 5);
2169/// fc->SetParameters(0, 1);
2170/// fc->Draw();
2171/// ~~~
2172///
2173/// \author Anna Kreshuk
2174
2175Double_t TMath::CauchyDist(Double_t x, Double_t t, Double_t s)
2176{
2177 Double_t temp= (x-t)*(x-t)/(s*s);
2178 Double_t result = 1/(s*TMath::Pi()*(1+temp));
2179 return result;
2180}
2181
2182////////////////////////////////////////////////////////////////////////////////
2183/// Evaluate the quantiles of the chi-squared probability distribution function.
2184/// Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
2185/// .
2186/// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
2187///
2188/// \param[in] p the probability value, at which the quantile is computed
2189/// \param[in] ndf number of degrees of freedom
2190///
2191/// \author Anna Kreshuk
2192
2193Double_t TMath::ChisquareQuantile(Double_t p, Double_t ndf)
2194{
2195 Double_t c[]={0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
2196 4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
2197 84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
2198 210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
2199 462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
2200 932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
2201 2520.0, 5040.0};
2202 Double_t e = 5e-7;
2203 Double_t aa = 0.6931471806;
2204 Int_t maxit = 20;
2205 Double_t ch, p1, p2, q, t, a, b, x;
2206 Double_t s1, s2, s3, s4, s5, s6;
2207
2208 if (ndf <= 0) return 0;
2209
2210 Double_t g = TMath::LnGamma(0.5*ndf);
2211
2212 Double_t xx = 0.5 * ndf;
2213 Double_t cp = xx - 1;
2214 if (ndf >= TMath::Log(p)*(-c[5])){
2215 //starting approximation for ndf less than or equal to 0.32
2216 if (ndf > c[3]) {
2217 x = TMath::NormQuantile(p);
2218 //starting approximation using Wilson and Hilferty estimate
2219 p1=c[2]/ndf;
2220 ch = ndf*TMath::Power((x*TMath::Sqrt(p1) + 1 - p1), 3);
2221 if (ch > c[6]*ndf + 6)
2222 ch = -2 * (TMath::Log(1-p) - cp * TMath::Log(0.5 * ch) + g);
2223 } else {
2224 ch = c[4];
2225 a = TMath::Log(1-p);
2226 do{
2227 q = ch;
2228 p1 = 1 + ch * (c[7]+ch);
2229 p2 = ch * (c[9] + ch * (c[8] + ch));
2230 t = -0.5 + (c[7] + 2 * ch) / p1 - (c[9] + ch * (c[10] + 3 * ch)) / p2;
2231 ch = ch - (1 - TMath::Exp(a + g + 0.5 * ch + cp * aa) *p2 / p1) / t;
2232 }while (TMath::Abs(q/ch - 1) > c[1]);
2233 }
2234 } else {
2235 ch = TMath::Power((p * xx * TMath::Exp(g + xx * aa)),(1./xx));
2236 if (ch < e) return ch;
2237 }
2238//call to algorithm AS 239 and calculation of seven term Taylor series
2239 for (Int_t i=0; i<maxit; i++){
2240 q = ch;
2241 p1 = 0.5 * ch;
2242 p2 = p - TMath::Gamma(xx, p1);
2243
2244 t = p2 * TMath::Exp(xx * aa + g + p1 - cp * TMath::Log(ch));
2245 b = t / ch;
2246 a = 0.5 * t - b * cp;
2247 s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] +c[11] * a))))) / c[24];
2248 s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
2249 s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
2250 s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
2251 s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
2252 s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
2253 ch = ch + t * (1 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
2254 if (TMath::Abs(q / ch - 1) > e) break;
2255 }
2256 return ch;
2257}
2258
2259////////////////////////////////////////////////////////////////////////////////
2260/// Computes the density function of F-distribution
2261/// (probability function, integral of density, is computed in FDistI).
2262///
2263/// Parameters N and M stand for degrees of freedom of chi-squares
2264/// mentioned above parameter F is the actual variable x of the
2265/// density function p(x) and the point at which the density function
2266/// is calculated.
2267///
2268/// ### About F distribution:
2269/// F-distribution arises in testing whether two random samples
2270/// have the same variance. It is the ratio of two chi-square
2271/// distributions, with N and M degrees of freedom respectively,
2272/// where each chi-square is first divided by it's number of degrees
2273/// of freedom.
2274///
2275/// \author Anna Kreshuk
2276
2277Double_t TMath::FDist(Double_t F, Double_t N, Double_t M)
2278{
2279 return ::ROOT::Math::fdistribution_pdf(F,N,M);
2280}
2281
2282////////////////////////////////////////////////////////////////////////////////
2283/// Calculates the cumulative distribution function of F-distribution
2284/// (see ROOT::Math::fdistribution_cdf).
2285/// This function occurs in the statistical test of whether two observed
2286/// samples have the same variance. For this test a certain statistic F,
2287/// the ratio of observed dispersion of the first sample to that of the
2288/// second sample, is calculated. N and M stand for numbers of degrees
2289/// of freedom in the samples 1-FDistI() is the significance level at
2290/// which the hypothesis "1 has smaller variance than 2" can be rejected.
2291/// A small numerical value of 1 - FDistI() implies a very significant
2292/// rejection, in turn implying high confidence in the hypothesis
2293/// "1 has variance greater than 2".
2294///
2295/// \author Anna Kreshuk
2296
2302
2303////////////////////////////////////////////////////////////////////////////////
2304/// Computes the density function of Gamma distribution at point x.
2305///
2306/// \param[in] x evaluation point
2307/// \param[in] gamma shape parameter
2308/// \param[in] mu location parameter
2309/// \param[in] beta scale parameter
2310///
2311/// The definition can be found in "Engineering Statistics Handbook" on site
2312/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
2313/// use now implementation in ROOT::Math::gamma_pdf
2314///
2315/// Begin_Macro
2316/// {
2317/// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2318///
2319/// c1->SetLogy();
2320/// c1->SetGridx();
2321/// c1->SetGridy();
2322///
2323/// TF1 *gdist = new TF1("gdist", "TMath::GammaDist(x, [0], [1], [2])", 0, 10);
2324///
2325/// gdist->SetParameters(0.5, 0., 1.);
2326/// gdist->SetLineColor(2);
2327/// TF1 *gdist1 = gdist->DrawCopy("L");
2328/// gdist->SetParameters(1.0, 0., 1.);
2329/// gdist->SetLineColor(3);
2330/// TF1 *gdist2 = gdist->DrawCopy("LSAME");
2331/// gdist->SetParameters(2.0, 0., 1.);
2332/// gdist->SetLineColor(4);
2333/// TF1 *gdist3 = gdist->DrawCopy("LSAME");
2334/// gdist->SetParameters(5.0, 0., 1.);
2335/// gdist->SetLineColor(6);
2336/// TF1 *gdist4 = gdist->DrawCopy("LSAME");
2337///
2338/// auto legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2339/// legend->AddEntry(gdist1, "gamma = 0.5, mu = 0, beta = 1", "L");
2340/// legend->AddEntry(gdist2, "gamma = 1.0, mu = 0, beta = 1", "L");
2341/// legend->AddEntry(gdist3, "gamma = 2.0, mu = 0, beta = 1", "L");
2342/// legend->AddEntry(gdist4, "gamma = 5.0, mu = 0, beta = 1", "L");
2343/// legend->Draw();
2344/// }
2345/// End_Macro
2346
2347Double_t TMath::GammaDist(Double_t x, Double_t gamma, Double_t mu, Double_t beta)
2348{
2349 if ((x<mu) || (gamma<=0) || (beta <=0)) {
2350 Error("TMath::GammaDist", "illegal parameter values x = %f , gamma = %f beta = %f",x,gamma,beta);
2351 return 0;
2352 }
2353 return ::ROOT::Math::gamma_pdf(x, gamma, beta, mu);
2354}
2355
2356////////////////////////////////////////////////////////////////////////////////
2357/// Computes the probability density function of Laplace distribution
2358/// at point x, with location parameter alpha and shape parameter beta.
2359/// By default, alpha=0, beta=1
2360/// This distribution is known under different names, most common is
2361/// double exponential distribution, but it also appears as
2362/// the two-tailed exponential or the bilateral exponential distribution
2363
2364Double_t TMath::LaplaceDist(Double_t x, Double_t alpha, Double_t beta)
2365{
2366 Double_t temp;
2367 temp = TMath::Exp(-TMath::Abs((x-alpha)/beta));
2368 temp /= (2.*beta);
2369 return temp;
2370}
2371
2372////////////////////////////////////////////////////////////////////////////////
2373/// Computes the cumulative distribution function (lower tail integral)
2374/// of Laplace distribution at point x, with location parameter alpha and shape parameter beta.
2375/// By default, alpha=0, beta=1
2376/// This distribution is known under different names, most common is
2377/// double exponential distribution, but it also appears as
2378/// the two-tailed exponential or the bilateral exponential distribution
2379
2380Double_t TMath::LaplaceDistI(Double_t x, Double_t alpha, Double_t beta)
2381{
2382 Double_t temp;
2383 if (x<=alpha){
2384 temp = 0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2385 } else {
2386 temp = 1-0.5*TMath::Exp(-TMath::Abs((x-alpha)/beta));
2387 }
2388 return temp;
2389}
2390
2391////////////////////////////////////////////////////////////////////////////////
2392/// Computes the density of LogNormal distribution at point x.
2393/// Variable X has lognormal distribution if Y=Ln(X) has normal distribution
2394///
2395/// \param[in] x is the evaluation point
2396/// \param[in] sigma is the shape parameter
2397/// \param[in] theta is the location parameter
2398/// \param[in] m is the scale parameter
2399///
2400/// The formula was taken from "Engineering Statistics Handbook" on site
2401/// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm
2402/// Implementation using ROOT::Math::lognormal_pdf
2403///
2404/// Begin_Macro
2405/// {
2406/// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2407///
2408/// c1->SetLogy();
2409/// c1->SetGridx();
2410/// c1->SetGridy();
2411///
2412/// TF1 *logn = new TF1("logn", "TMath::LogNormal(x, [0], [1], [2])", 0, 5);
2413/// logn->SetMinimum(1e-3);
2414///
2415/// logn->SetParameters(0.5, 0., 1.);
2416/// logn->SetLineColor(2);
2417/// TF1 *logn1 = logn->DrawCopy("L");
2418/// logn->SetParameters(1.0, 0., 1.);
2419/// logn->SetLineColor(3);
2420/// TF1 *logn2 = logn->DrawCopy("LSAME");
2421/// logn->SetParameters(2.0, 0., 1.);
2422/// logn->SetLineColor(4);
2423/// TF1 *logn3 = logn->DrawCopy("LSAME");
2424/// logn->SetParameters(5.0, 0., 1.);
2425/// logn->SetLineColor(6);
2426/// TF1 *logn4 = logn->DrawCopy("LSAME");
2427///
2428/// auto legend = new TLegend(0.15, 0.15, 0.5, 0.35);
2429/// legend->AddEntry(logn1, "sigma = 0.5, theta = 0, m = 1", "L");
2430/// legend->AddEntry(logn2, "sigma = 1.0, theta = 0, m = 1", "L");
2431/// legend->AddEntry(logn3, "sigma = 2.0, theta = 0, m = 1", "L");
2432/// legend->AddEntry(logn4, "sigma = 5.0, theta = 0, m = 1", "L");
2433/// legend->Draw();
2434/// }
2435/// End_Macro
2436
2437Double_t TMath::LogNormal(Double_t x, Double_t sigma, Double_t theta, Double_t m)
2438{
2439 if ((x<theta) || (sigma<=0) || (m<=0)) {
2440 Error("TMath::Lognormal", "illegal parameter values");
2441 return 0;
2442 }
2443 // lognormal_pdf uses same definition of http://en.wikipedia.org/wiki/Log-normal_distribution
2444 // where mu = log(m)
2445
2446 return ::ROOT::Math::lognormal_pdf(x, TMath::Log(m), sigma, theta);
2447
2448}
2449
2450////////////////////////////////////////////////////////////////////////////////
2451/// Computes quantiles for standard normal distribution N(0, 1)
2452/// at probability p
2453///
2454/// ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.
2455
2456Double_t TMath::NormQuantile(Double_t p)
2457{
2458 if ((p<=0)||(p>=1)) {
2459 Error("TMath::NormQuantile", "probability outside (0, 1)");
2460 return 0;
2461 }
2462
2463 Double_t a0 = 3.3871328727963666080e0;
2464 Double_t a1 = 1.3314166789178437745e+2;
2465 Double_t a2 = 1.9715909503065514427e+3;
2466 Double_t a3 = 1.3731693765509461125e+4;
2467 Double_t a4 = 4.5921953931549871457e+4;
2468 Double_t a5 = 6.7265770927008700853e+4;
2469 Double_t a6 = 3.3430575583588128105e+4;
2470 Double_t a7 = 2.5090809287301226727e+3;
2471 Double_t b1 = 4.2313330701600911252e+1;
2472 Double_t b2 = 6.8718700749205790830e+2;
2473 Double_t b3 = 5.3941960214247511077e+3;
2474 Double_t b4 = 2.1213794301586595867e+4;
2475 Double_t b5 = 3.9307895800092710610e+4;
2476 Double_t b6 = 2.8729085735721942674e+4;
2477 Double_t b7 = 5.2264952788528545610e+3;
2478 Double_t c0 = 1.42343711074968357734e0;
2479 Double_t c1 = 4.63033784615654529590e0;
2480 Double_t c2 = 5.76949722146069140550e0;
2481 Double_t c3 = 3.64784832476320460504e0;
2482 Double_t c4 = 1.27045825245236838258e0;
2483 Double_t c5 = 2.41780725177450611770e-1;
2484 Double_t c6 = 2.27238449892691845833e-2;
2485 Double_t c7 = 7.74545014278341407640e-4;
2486 Double_t d1 = 2.05319162663775882187e0;
2487 Double_t d2 = 1.67638483018380384940e0;
2488 Double_t d3 = 6.89767334985100004550e-1;
2489 Double_t d4 = 1.48103976427480074590e-1;
2490 Double_t d5 = 1.51986665636164571966e-2;
2491 Double_t d6 = 5.47593808499534494600e-4;
2492 Double_t d7 = 1.05075007164441684324e-9;
2493 Double_t e0 = 6.65790464350110377720e0;
2494 Double_t e1 = 5.46378491116411436990e0;
2495 Double_t e2 = 1.78482653991729133580e0;
2496 Double_t e3 = 2.96560571828504891230e-1;
2497 Double_t e4 = 2.65321895265761230930e-2;
2498 Double_t e5 = 1.24266094738807843860e-3;
2499 Double_t e6 = 2.71155556874348757815e-5;
2500 Double_t e7 = 2.01033439929228813265e-7;
2501 Double_t f1 = 5.99832206555887937690e-1;
2502 Double_t f2 = 1.36929880922735805310e-1;
2503 Double_t f3 = 1.48753612908506148525e-2;
2504 Double_t f4 = 7.86869131145613259100e-4;
2505 Double_t f5 = 1.84631831751005468180e-5;
2506 Double_t f6 = 1.42151175831644588870e-7;
2507 Double_t f7 = 2.04426310338993978564e-15;
2508
2509 Double_t split1 = 0.425;
2510 Double_t split2=5.;
2511 Double_t konst1=0.180625;
2512 Double_t konst2=1.6;
2513
2514 Double_t q, r, quantile;
2515 q=p-0.5;
2516 if (TMath::Abs(q)<split1) {
2517 r=konst1-q*q;
2518 quantile = q* (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
2519 * r + a2) * r + a1) * r + a0) /
2520 (((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
2521 * r + b2) * r + b1) * r + 1.);
2522 } else {
2523 if(q<0) r=p;
2524 else r=1-p;
2525 //error case
2526 if (r<=0)
2527 quantile=0;
2528 else {
2529 r=TMath::Sqrt(-TMath::Log(r));
2530 if (r<=split2) {
2531 r=r-konst2;
2532 quantile=(((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
2533 * r + c2) * r + c1) * r + c0) /
2534 (((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
2535 * r + d2) * r + d1) * r + 1);
2536 } else{
2537 r=r-split2;
2538 quantile=(((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
2539 * r + e2) * r + e1) * r + e0) /
2540 (((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
2541 * r + f2) * r + f1) * r + 1);
2542 }
2543 if (q<0) quantile=-quantile;
2544 }
2545 }
2546 return quantile;
2547}
2548
2549////////////////////////////////////////////////////////////////////////////////
2550/// Simple recursive algorithm to find the permutations of
2551/// n natural numbers, not necessarily all distinct
2552/// adapted from CERNLIB routine PERMU.
2553/// The input array has to be initialised with a non descending
2554/// sequence. The method returns kFALSE when
2555/// all combinations are exhausted.
2556
2557Bool_t TMath::Permute(Int_t n, Int_t *a)
2558{
2559 Int_t i,itmp;
2560 Int_t i1=-1;
2561
2562 // find rightmost upward transition
2563 for(i=n-2; i>-1; i--) {
2564 if(a[i]<a[i+1]) {
2565 i1=i;
2566 break;
2567 }
2568 }
2569 // no more upward transitions, end of the story
2570 if(i1==-1) return kFALSE;
2571 else {
2572 // find lower right element higher than the lower
2573 // element of the upward transition
2574 for(i=n-1;i>i1;i--) {
2575 if(a[i] > a[i1]) {
2576 // swap the two
2577 itmp=a[i1];
2578 a[i1]=a[i];
2579 a[i]=itmp;
2580 break;
2581 }
2582 }
2583 // order the rest, in fact just invert, as there
2584 // are only downward transitions from here on
2585 for(i=0;i<(n-i1-1)/2;i++) {
2586 itmp=a[i1+i+1];
2587 a[i1+i+1]=a[n-i-1];
2588 a[n-i-1]=itmp;
2589 }
2590 }
2591 return kTRUE;
2592}
2593
2594////////////////////////////////////////////////////////////////////////////////
2595/// Computes density function for Student's t- distribution
2596/// (the probability function (integral of density) is computed in StudentI).
2597///
2598/// First parameter stands for x - the actual variable of the
2599/// density function p(x) and the point at which the density is calculated.
2600/// Second parameter stands for number of degrees of freedom.
2601///
2602/// About Student distribution:
2603/// Student's t-distribution is used for many significance tests, for example,
2604/// for the Student's t-tests for the statistical significance of difference
2605/// between two sample means and for confidence intervals for the difference
2606/// between two population means.
2607///
2608/// Example: suppose we have a random sample of size n drawn from normal
2609/// distribution with mean Mu and st.deviation Sigma. Then the variable
2610///
2611/// t = (sample_mean - Mu)/(sample_deviation / sqrt(n))
2612///
2613/// has Student's t-distribution with n-1 degrees of freedom.
2614///
2615/// NOTE that this function's second argument is number of degrees of freedom,
2616/// not the sample size.
2617///
2618/// As the number of degrees of freedom grows, t-distribution approaches
2619/// Normal(0,1) distribution.
2620///
2621/// \author Anna Kreshuk
2622
2623Double_t TMath::Student(Double_t T, Double_t ndf)
2624{
2625 if (ndf < 1) {
2626 return 0;
2627 }
2628
2629 Double_t r = ndf;
2630 Double_t rh = 0.5*r;
2631 Double_t rh1 = rh + 0.5;
2632 Double_t denom = TMath::Sqrt(r*TMath::Pi())*TMath::Gamma(rh)*TMath::Power(1+T*T/r, rh1);
2633 return TMath::Gamma(rh1)/denom;
2634}
2635
2636////////////////////////////////////////////////////////////////////////////////
2637/// Calculates the cumulative distribution function of Student's
2638/// t-distribution second parameter stands for number of degrees of freedom,
2639/// not for the number of samples
2640/// if x has Student's t-distribution, the function returns the probability of
2641/// x being less than T.
2642/// This is equivalent to ROOT::Math::tdistribution_cdf(T,ndf)
2643///
2644/// \author Anna Kreshuk
2645
2646Double_t TMath::StudentI(Double_t T, Double_t ndf)
2647{
2648 Double_t r = ndf;
2649
2650 Double_t si = (T>0) ?
2651 (1 - 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5)) :
2652 0.5*BetaIncomplete((r/(r + T*T)), r*0.5, 0.5);
2653 return si;
2654}
2655
2656////////////////////////////////////////////////////////////////////////////////
2657/// Computes quantiles of the Student's t-distribution
2658/// 1st argument is the probability, at which the quantile is computed
2659/// 2nd argument - the number of degrees of freedom of the
2660/// Student distribution
2661/// When the 3rd argument lower_tail is kTRUE (default)-
2662/// the algorithm returns such x0, that
2663///
2664/// P(x < x0)=p
2665///
2666/// upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that
2667///
2668/// P(x > x0)=p
2669///
2670/// the algorithm was taken from:
2671/// G.W.Hill, "Algorithm 396, Student's t-quantiles"
2672/// "Communications of the ACM", 13(10), October 1970
2673
2674Double_t TMath::StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail)
2675{
2676 Double_t quantile;
2677 Double_t temp;
2678 Bool_t neg;
2679 Double_t q;
2680 if (ndf<1 || p>=1 || p<=0) {
2681 Error("TMath::StudentQuantile", "illegal parameter values");
2682 return 0;
2683 }
2684 if ((lower_tail && p>0.5)||(!lower_tail && p<0.5)){
2685 neg=kFALSE;
2686 q=2*(lower_tail ? (1-p) : p);
2687 } else {
2688 neg=kTRUE;
2689 q=2*(lower_tail? p : (1-p));
2690 }
2691
2692 if ((ndf-1)<1e-8) {
2693 temp=TMath::PiOver2()*q;
2694 quantile = TMath::Cos(temp)/TMath::Sin(temp);
2695 } else {
2696 if ((ndf-2)<1e-8) {
2697 quantile = TMath::Sqrt(2./(q*(2-q))-2);
2698 } else {
2699 Double_t a=1./(ndf-0.5);
2700 Double_t b=48./(a*a);
2701 Double_t c=((20700*a/b -98)*a-16)*a+96.36;
2702 Double_t d=((94.5/(b+c)-3.)/b+1)*TMath::Sqrt(a*TMath::PiOver2())*ndf;
2703 Double_t x=q*d;
2704 Double_t y=TMath::Power(x, (2./ndf));
2705 if (y>0.05+a){
2706 //asymptotic inverse expansion about normal
2707 x=TMath::NormQuantile(q*0.5);
2708 y=x*x;
2709 if (ndf<5) c+=0.3*(ndf-4.5)*(x+0.6);
2710 c+=(((0.05*d*x-5.)*x-7.)*x-2.)*x +b;
2711 y=(((((0.4*y+6.3)*y+36.)*y+94.5)/c - y-3.)/b+1)*x;
2712 y=a*y*y;
2713 if(y>0.002) y = TMath::Exp(y)-1;
2714 else y += 0.5*y*y;
2715 } else {
2716 y=((1./(((ndf+6.)/(ndf*y)-0.089*d-0.822)*(ndf+2.)*3)+0.5/(ndf+4.))*y-1.)*
2717 (ndf+1.)/(ndf+2.)+1/y;
2718 }
2719 quantile = TMath::Sqrt(ndf*y);
2720 }
2721 }
2722 if(neg) quantile=-quantile;
2723 return quantile;
2724}
2725
2726////////////////////////////////////////////////////////////////////////////////
2727/// Returns the value of the Vavilov probability density function
2728///
2729/// \param[in] x the point were the density function is evaluated
2730/// \param[in] kappa value of kappa (distribution parameter)
2731/// \param[in] beta2 value of beta2 (distribution parameter)
2732///
2733/// The algorithm was taken from the CernLib function vavden(G115)
2734/// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2735/// Nucl.Instr. and Meth. B47(1990), 215-224
2736///
2737/// Accuracy: quote from the reference above:
2738///
2739/// "The results of our code have been compared with the values of the Vavilov
2740/// density function computed numerically in an accurate way: our approximation
2741/// shows a difference of less than 3% around the peak of the density function, slowly
2742/// increasing going towards the extreme tails to the right and to the left"
2743///
2744/// For a more accurate implementation see the documentation in the Vavilov class and
2745/// ROOT::Math::vavilov_accurate_pdf
2746///
2747/// Begin_Macro
2748/// {
2749/// TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
2750///
2751/// c1->SetGridx();
2752/// c1->SetGridy();
2753///
2754/// TF1 *vavilov = new TF1("vavilov", "TMath::Vavilov(x, [0], [1])", -3, 11);
2755///
2756/// vavilov->SetParameters(0.5, 0.);
2757/// vavilov->SetLineColor(2);
2758/// TF1 *vavilov1 = vavilov->DrawCopy("L");
2759/// vavilov->SetParameters(0.3, 0.);
2760/// vavilov->SetLineColor(3);
2761/// TF1 *vavilov2 = vavilov->DrawCopy("LSAME");
2762/// vavilov->SetParameters(0.2, 0.);
2763/// vavilov->SetLineColor(4);
2764/// TF1 *vavilov3 = vavilov->DrawCopy("LSAME");
2765/// vavilov->SetParameters(0.1, 0.);
2766/// vavilov->SetLineColor(6);
2767/// TF1 *vavilov4 = vavilov->DrawCopy("LSAME");
2768///
2769/// auto legend = new TLegend(0.5, 0.65, 0.85, 0.85);
2770/// legend->AddEntry(vavilov1, "kappa = 0.5, beta2 = 0", "L");
2771/// legend->AddEntry(vavilov2, "kappa = 0.3, beta2 = 0", "L");
2772/// legend->AddEntry(vavilov3, "kappa = 0.2, beta2 = 0", "L");
2773/// legend->AddEntry(vavilov4, "kappa = 0.1, beta2 = 0", "L");
2774/// legend->Draw();
2775/// }
2776/// End_Macro
2777
2778Double_t TMath::Vavilov(Double_t x, Double_t kappa, Double_t beta2)
2779{
2780 Double_t *ac = new Double_t[14];
2781 Double_t *hc = new Double_t[9];
2782
2783 Int_t itype;
2784 Int_t npt;
2785 TMath::VavilovSet(kappa, beta2, false, nullptr, ac, hc, itype, npt);
2786 Double_t v = TMath::VavilovDenEval(x, ac, hc, itype);
2787 delete [] ac;
2788 delete [] hc;
2789 return v;
2790}
2791
2792////////////////////////////////////////////////////////////////////////////////
2793/// Returns the value of the Vavilov cumulative distribution function
2794/// (lower tail integral of the probability distribution function)
2795///
2796/// \param[in] x the point were the density function is evaluated
2797/// \param[in] kappa value of kappa (distribution parameter)
2798/// \param[in] beta2 value of beta2 (distribution parameter)
2799///
2800/// The algorithm was taken from the CernLib function vavden(G115)
2801///
2802/// Reference: A.Rotondi and P.Montagna, Fast Calculation of Vavilov distribution
2803/// Nucl.Instr. and Meth. B47(1990), 215-224
2804///
2805/// Accuracy: quote from the reference above:
2806///
2807/// "The results of our code have been compared with the values of the Vavilov
2808/// density function computed numerically in an accurate way: our approximation
2809/// shows a difference of less than 3% around the peak of the density function, slowly
2810/// increasing going towards the extreme tails to the right and to the left"
2811///
2812/// For a more accurate implementation see the documentation of the Vavilov class and the cumulative
2813/// ROOT::Math::vavilov_accurate_cdf
2814
2815Double_t TMath::VavilovI(Double_t x, Double_t kappa, Double_t beta2)
2816{
2817 Double_t *ac = new Double_t[14];
2818 Double_t *hc = new Double_t[9];
2819 Double_t *wcm = new Double_t[201];
2820 Int_t itype;
2821 Int_t npt;
2822 Int_t k;
2823 Double_t xx, v;
2824 TMath::VavilovSet(kappa, beta2, true, wcm, ac, hc, itype, npt);
2825 if (x < ac[0]) v = 0;
2826 else if (x >=ac[8]) v = 1;
2827 else {
2828 xx = x - ac[0];
2829 k = Int_t(xx*ac[10]);
2830 v = TMath::Min(wcm[k] + (xx - k*ac[9])*(wcm[k+1]-wcm[k])*ac[10], 1.);
2831 }
2832 delete [] ac;
2833 delete [] hc;
2834 delete [] wcm;
2835 return v;
2836}
2837
2838////////////////////////////////////////////////////////////////////////////////
2839/// Returns the cumulative (lower tail integral) of the Landau distribution function at point x.
2840/// (see ROOT::Math::landau_cdf)
2841/// The algorithm was taken from the Cernlib function dislan(G110)
2842/// Reference: K.S.Kolbig and B.Schorr, "A program package for the Landau
2843/// distribution", Computer Phys.Comm., 31(1984), 97-111
2844
2845Double_t TMath::LandauI(Double_t x)
2846{
2847 return ::ROOT::Math::landau_cdf(x);
2848}
2849
2850
2851////////////////////////////////////////////////////////////////////////////////
2852/// Internal function, called by Vavilov and VavilovI
2853
2854void 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)
2855{
2856
2857 Double_t BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
2858 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
2859 BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
2860
2861 Double_t FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
2862 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
2863 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
2864
2865 Double_t FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
2866
2867 Double_t EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
2868 0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
2869
2870 Double_t U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
2871 0. , 0.24880692e-1, 0.47404356e-2,
2872 -0.74445130e-3, 0.73225731e-2, 0. ,
2873 0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
2874
2875 Double_t U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
2876 0. , 0.42127077e-1, 0.73167928e-2,
2877 -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
2878 0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
2879
2880 Double_t U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
2881 0. , 0.23834176e-1, 0.21624675e-2,
2882 -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
2883 0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
2884
2885 Double_t U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
2886 -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
2887 0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
2888 -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
2889
2890 Double_t U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
2891 0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
2892 0.48736023e-3, 0.34850854e-2, 0. ,
2893 -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
2894
2895 Double_t U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
2896 -0.30188807e-2, -0.84479939e-3, 0. ,
2897 0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
2898 -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
2899
2900 Double_t U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
2901 -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
2902 0. , 0.50505298e+0};
2903 Double_t U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
2904 -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
2905 -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
2906 0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
2907
2908 Double_t V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
2909 0. , 0.45091424e-1, 0.80559636e-2,
2910 -0.38974523e-2, 0. , -0.30634124e-2,
2911 0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
2912
2913 Double_t V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
2914 0. , 0.12693873e+0, 0.22999801e-1,
2915 -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
2916 0. , 0.19716857e-1, 0.32596742e-2};
2917
2918 Double_t V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
2919 -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
2920 -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
2921 0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
2922
2923 Double_t V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
2924 0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
2925 0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
2926 -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
2927
2928 Double_t V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
2929 -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
2930 0.56856517e-2, -0.13438433e-1, 0. ,
2931 0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
2932
2933 Double_t V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
2934 0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
2935 0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
2936 -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
2937
2938 Double_t V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
2939 0. , -0.12218009e+1, -0.32324120e+0,
2940 -0.27373591e-1, 0.12173464e+0, 0. ,
2941 0. , 0.40917471e-1};
2942
2943 Double_t V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
2944 0. , -0.18160479e+1, -0.50919193e+0,
2945 -0.51384654e-1, 0.21413992e+0, 0. ,
2946 0. , 0.66596366e-1};
2947
2948 Double_t W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
2949 -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
2950 -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
2951 0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
2952
2953 Double_t W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
2954 -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
2955 0. , 0.23020158e-1, 0.50574313e-2,
2956 0.94550140e-2, 0.19300232e-1};
2957
2958 Double_t W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
2959 -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
2960 -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
2961 0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
2962
2963 Double_t W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
2964 0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
2965 0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
2966 -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
2967
2968 Double_t W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
2969 0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
2970 0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
2971 -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
2972
2973 Double_t W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
2974 0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
2975 0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
2976 -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
2977
2978 Double_t W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
2979 0. , -0.14540925e+1, -0.39529833e+0,
2980 -0.44293243e-1, 0.88741049e-1};
2981
2982 itype = 0;
2983 if (rkappa <0.01 || rkappa >12) {
2984 Error("Vavilov distribution", "illegal value of kappa");
2985 return;
2986 }
2987
2988 Double_t DRK[6];
2989 Double_t DSIGM[6];
2990 Double_t ALFA[8];
2991 Int_t j;
2992 Double_t x, y, xx, yy, x2, x3, y2, y3, xy, p2, p3, q2, q3, pq;
2993 if (rkappa >= 0.29) {
2994 itype = 1;
2995 npt = 100;
2996 Double_t wk = 1./TMath::Sqrt(rkappa);
2997
2998 AC[0] = (-0.032227*beta2-0.074275)*rkappa + (0.24533*beta2+0.070152)*wk + (-0.55610*beta2-3.1579);
2999 AC[8] = (-0.013483*beta2-0.048801)*rkappa + (-1.6921*beta2+8.3656)*wk + (-0.73275*beta2-3.5226);
3000 DRK[1] = wk*wk;
3001 DSIGM[1] = TMath::Sqrt(rkappa/(1-0.5*beta2));
3002 for (j=1; j<=4; j++) {
3003 DRK[j+1] = DRK[1]*DRK[j];
3004 DSIGM[j+1] = DSIGM[1]*DSIGM[j];
3005 ALFA[j+1] = (FNINV[j]-beta2*FNINV[j+1])*DRK[j];
3006 }
3007 HC[0]=TMath::Log(rkappa)+beta2+0.42278434;
3008 HC[1]=DSIGM[1];
3009 HC[2]=ALFA[3]*DSIGM[3];
3010 HC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
3011 HC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*HC[2];
3012 HC[5]=HC[2]*HC[2];
3013 HC[6]=HC[2]*HC[3];
3014 HC[7]=HC[2]*HC[5];
3015 for (j=2; j<=7; j++)
3016 HC[j]*=EDGEC[j];
3017 HC[8]=0.39894228*HC[1];
3018 }
3019 else if (rkappa >=0.22) {
3020 itype = 2;
3021 npt = 150;
3022 x = 1+(rkappa-BKMXX3)*FBKX3;
3023 y = 1+(TMath::Sqrt(beta2)-BKMXY3)*FBKY3;
3024 xx = 2*x;
3025 yy = 2*y;
3026 x2 = xx*x-1;
3027 x3 = xx*x2-x;
3028 y2 = yy*y-1;
3029 y3 = yy*y2-y;
3030 xy = x*y;
3031 p2 = x2*y;
3032 p3 = x3*y;
3033 q2 = y2*x;
3034 q3 = y3*x;
3035 pq = x2*y2;
3036 AC[1] = W1[1] + W1[2]*x + W1[4]*x3 + W1[5]*y + W1[6]*y2 + W1[7]*y3 +
3037 W1[8]*xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
3038 AC[2] = W2[1] + W2[2]*x + W2[3]*x2 + W2[4]*x3 + W2[5]*y + W2[6]*y2 +
3039 W2[8]*xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
3040 AC[3] = W3[1] + W3[3]*x2 + W3[4]*x3 + W3[5]*y + W3[6]*y2 + W3[7]*y3 +
3041 W3[8]*xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
3042 AC[4] = W4[1] + W4[2]*x + W4[3]*x2 + W4[4]*x3 + W4[5]*y + W4[6]*y2 + W4[7]*y3 +
3043 W4[8]*xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
3044 AC[5] = W5[1] + W5[2]*x + W5[4]*x3 + W5[5]*y + W5[6]*y2 + W5[7]*y3 +
3045 W5[8]*xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
3046 AC[6] = W6[1] + W6[2]*x + W6[3]*x2 + W6[4]*x3 + W6[5]*y + W6[6]*y2 + W6[7]*y3 +
3047 W6[8]*xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
3048 AC[8] = W8[1] + W8[2]*x + W8[3]*x2 + W8[5]*y + W8[6]*y2 + W8[7]*y3 + W8[8]*xy;
3049 AC[0] = -3.05;
3050 } else if (rkappa >= 0.12) {
3051 itype = 3;
3052 npt = 200;
3053 x = 1 + (rkappa-BKMXX2)*FBKX2;
3054 y = 1 + (TMath::Sqrt(beta2)-BKMXY2)*FBKY2;
3055 xx = 2*x;
3056 yy = 2*y;
3057 x2 = xx*x-1;
3058 x3 = xx*x2-x;
3059 y2 = yy*y-1;
3060 y3 = yy*y2-y;
3061 xy = x*y;
3062 p2 = x2*y;
3063 p3 = x3*y;
3064 q2 = y2*x;
3065 q3 = y3*x;
3066 pq = x2*y2;
3067 AC[1] = V1[1] + V1[2]*x + V1[3]*x2 + V1[5]*y + V1[6]*y2 + V1[7]*y3 +
3068 V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
3069 AC[2] = V2[1] + V2[2]*x + V2[3]*x2 + V2[5]*y + V2[6]*y2 + V2[7]*y3 +
3070 V2[8]*xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
3071 AC[3] = V3[1] + V3[2]*x + V3[3]*x2 + V3[4]*x3 + V3[5]*y + V3[6]*y2 + V3[7]*y3 +
3072 V3[8]*xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
3073 AC[4] = V4[1] + V4[2]*x + V4[3]*x2 + V4[4]*x3 + V4[5]*y + V4[6]*y2 + V4[7]*y3 +
3074 V4[8]*xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
3075 AC[5] = V5[1] + V5[2]*x + V5[3]*x2 + V5[4]*x3 + V5[5]*y + V5[6]*y2 + V5[7]*y3 +
3076 V5[8]*xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
3077 AC[6] = V6[1] + V6[2]*x + V6[3]*x2 + V6[4]*x3 + V6[5]*y + V6[6]*y2 + V6[7]*y3 +
3078 V6[8]*xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
3079 AC[7] = V7[1] + V7[2]*x + V7[3]*x2 + V7[5]*y + V7[6]*y2 + V7[7]*y3 +
3080 V7[8]*xy + V7[11]*q2;
3081 AC[8] = V8[1] + V8[2]*x + V8[3]*x2 + V8[5]*y + V8[6]*y2 + V8[7]*y3 +
3082 V8[8]*xy + V8[11]*q2;
3083 AC[0] = -3.04;
3084 } else {
3085 itype = 4;
3086 if (rkappa >=0.02) itype = 3;
3087 npt = 200;
3088 x = 1+(rkappa-BKMXX1)*FBKX1;
3089 y = 1+(TMath::Sqrt(beta2)-BKMXY1)*FBKY1;
3090 xx = 2*x;
3091 yy = 2*y;
3092 x2 = xx*x-1;
3093 x3 = xx*x2-x;
3094 y2 = yy*y-1;
3095 y3 = yy*y2-y;
3096 xy = x*y;
3097 p2 = x2*y;
3098 p3 = x3*y;
3099 q2 = y2*x;
3100 q3 = y3*x;
3101 pq = x2*y2;
3102 if (itype==3){
3103 AC[1] = U1[1] + U1[2]*x + U1[3]*x2 + U1[5]*y + U1[6]*y2 + U1[7]*y3 +
3104 U1[8]*xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
3105 AC[2] = U2[1] + U2[2]*x + U2[3]*x2 + U2[5]*y + U2[6]*y2 + U2[7]*y3 +
3106 U2[8]*xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
3107 AC[3] = U3[1] + U3[2]*x + U3[3]*x2 + U3[5]*y + U3[6]*y2 + U3[7]*y3 +
3108 U3[8]*xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
3109 AC[4] = U4[1] + U4[2]*x + U4[3]*x2 + U4[4]*x3 + U4[5]*y + U4[6]*y2 + U4[7]*y3 +
3110 U4[8]*xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
3111 AC[5] = U5[1] + U5[2]*x + U5[3]*x2 + U5[4]*x3 + U5[5]*y + U5[6]*y2 + U5[7]*y3 +
3112 U5[8]*xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
3113 AC[6] = U6[1] + U6[2]*x + U6[3]*x2 + U6[4]*x3 + U6[5]*y + U6[7]*y3 +
3114 U6[8]*xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
3115 AC[7] = U7[1] + U7[2]*x + U7[3]*x2 + U7[4]*x3 + U7[5]*y + U7[6]*y2 + U7[8]*xy;
3116 }
3117 AC[8] = U8[1] + U8[2]*x + U8[3]*x2 + U8[4]*x3 + U8[5]*y + U8[6]*y2 + U8[7]*y3 +
3118 U8[8]*xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
3119 AC[0] = -3.03;
3120 }
3121
3122 AC[9] = (AC[8] - AC[0])/npt;
3123 AC[10] = 1./AC[9];
3124 if (itype == 3) {
3125 x = (AC[7]-AC[8])/(AC[7]*AC[8]);
3126 y = 1./TMath::Log(AC[8]/AC[7]);
3127 p2 = AC[7]*AC[7];
3128 AC[11] = p2*(AC[1]*TMath::Exp(-AC[2]*(AC[7]+AC[5]*p2)-
3129 AC[3]*TMath::Exp(-AC[4]*(AC[7]+AC[6]*p2)))-0.045*y/AC[7])/(1+x*y*AC[7]);
3130 AC[12] = (0.045+x*AC[11])*y;
3131 }
3132 if (itype == 4) AC[13] = 0.995/LandauI(AC[8]);
3133
3134 if (mode==0) return;
3135 //
3136 x = AC[0];
3137 WCM[0] = 0;
3138 Double_t fl, fu;
3139 Int_t k;
3140 fl = TMath::VavilovDenEval(x, AC, HC, itype);
3141 for (k=1; k<=npt; k++) {
3142 x += AC[9];
3143 fu = TMath::VavilovDenEval(x, AC, HC, itype);
3144 WCM[k] = WCM[k-1] + fl + fu;
3145 fl = fu;
3146 }
3147 x = 0.5*AC[9];
3148 for (k=1; k<=npt; k++)
3149 WCM[k]*=x;
3150}
3151
3152////////////////////////////////////////////////////////////////////////////////
3153/// Internal function, called by Vavilov and VavilovSet
3154
3155Double_t TMath::VavilovDenEval(Double_t rlam, Double_t *AC, Double_t *HC, Int_t itype)
3156{
3157 Double_t v = 0;
3158 if (rlam < AC[0] || rlam > AC[8])
3159 return 0;
3160 Int_t k;
3161 Double_t x, fn, s;
3162 Double_t h[10];
3163 if (itype ==1 ) {
3164 fn = 1;
3165 x = (rlam + HC[0])*HC[1];
3166 h[1] = x;
3167 h[2] = x*x -1;
3168 for (k=2; k<=8; k++) {
3169 fn++;
3170 h[k+1] = x*h[k]-fn*h[k-1];
3171 }
3172 s = 1 + HC[7]*h[9];
3173 for (k=2; k<=6; k++)
3174 s+=HC[k]*h[k+1];
3175 v = HC[8]*TMath::Exp(-0.5*x*x)*TMath::Max(s, 0.);
3176 }
3177 else if (itype == 2) {
3178 x = rlam*rlam;
3179 v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x) - AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3180 }
3181 else if (itype == 3) {
3182 if (rlam < AC[7]) {
3183 x = rlam*rlam;
3184 v = AC[1]*TMath::Exp(-AC[2]*(rlam+AC[5]*x)-AC[3]*TMath::Exp(-AC[4]*(rlam+AC[6]*x)));
3185 } else {
3186 x = 1./rlam;
3187 v = (AC[11]*x + AC[12])*x;
3188 }
3189 }
3190 else if (itype == 4) {
3191 v = AC[13]*TMath::Landau(rlam);
3192 }
3193 return v;
3194}
3195
3196
3197//explicitly instantiate template functions from VecCore
3198#ifdef R__HAS_VECCORE
3199#include <Math/Types.h>
3200template ROOT::Double_v vecCore::math::Sin(const ROOT::Double_v & x);
3201template ROOT::Double_v vecCore::math::Cos(const ROOT::Double_v & x);
3202template ROOT::Double_v vecCore::math::ASin(const ROOT::Double_v & x);
3203template ROOT::Double_v vecCore::math::ATan(const ROOT::Double_v & x);
3204template ROOT::Double_v vecCore::math::ATan2(const ROOT::Double_v & x,const ROOT::Double_v & y);
3205// missing in veccore
3206// template ROOT::Double_v vecCore::math::ACos(const ROOT::Double_v & x);
3207// template ROOT::Double_v vecCore::math::Sinh(const ROOT::Double_v & x);
3208// template ROOT::Double_v vecCore::math::Cosh(const ROOT::Double_v & x);
3209// template ROOT::Double_v vecCore::math::Tanh(const ROOT::Double_v & x);
3210// template ROOT::Double_v vecCore::math::Cbrt(const ROOT::Double_v & x);
3211#endif
d
#define d(i)
Definition RSha256.hxx:102
b
#define b(i)
Definition RSha256.hxx:100
c
#define c(i)
Definition RSha256.hxx:101
g
#define g(i)
Definition RSha256.hxx:105
a
#define a(i)
Definition RSha256.hxx:99
s1
#define s1(x)
Definition RSha256.hxx:91
h
#define h(i)
Definition RSha256.hxx:106
e
#define e(i)
Definition RSha256.hxx:103
Bool_t
bool Bool_t
Definition RtypesCore.h:63
Int_t
int Int_t
Definition RtypesCore.h:45
ULong_t
unsigned long ULong_t
Definition RtypesCore.h:55
Long_t
long Long_t
Definition RtypesCore.h:54
Float_t
float Float_t
Definition RtypesCore.h:57
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
Double_t
double Double_t
Definition RtypesCore.h:59
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
Option_t
const char Option_t
Definition RtypesCore.h:66
NamespaceImp
#define NamespaceImp(name)
Definition Rtypes.h:393
R__ASSERT
#define R__ASSERT(e)
Definition TError.h:118
Info
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:218
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
N
#define N
option
Option_t Option_t option
Definition TGWin32VirtualXProxy.cxx:44
w
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void w
Definition TGWin32VirtualXProxy.cxx:138
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 Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t del
Definition TGWin32VirtualXProxy.cxx:188
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 r
Definition TGWin32VirtualXProxy.cxx:168
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 Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Definition TGWin32VirtualXProxy.cxx:174
index
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
Definition TGWin32VirtualXProxy.cxx:168
x2
Option_t Option_t TPoint TPoint const char x2
Definition TGWin32VirtualXProxy.cxx:70
xy
Option_t Option_t TPoint xy
Definition TGWin32VirtualXProxy.cxx:62
mode
Option_t Option_t TPoint TPoint const char mode
Definition TGWin32VirtualXProxy.cxx:68
y2
Option_t Option_t TPoint TPoint const char y2
Definition TGWin32VirtualXProxy.cxx:70
type
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
Definition TGWin32VirtualXProxy.cxx:249
y1
Option_t Option_t TPoint TPoint const char y1
Definition TGWin32VirtualXProxy.cxx:70
i
Int_t i
Definition TGraphAsymmErrors.cxx:1306
q
float * q
Definition THbookFile.cxx:89
Hash
UInt_t Hash(const TString &s)
Definition TString.h:494
TString
Basic string class.
Definition TString.h:139
TString::ToUpper
void ToUpper()
Change string to upper case.
Definition TString.cxx:1195
TString::Hash
UInt_t Hash(ECaseCompare cmp=kExact) const
Return hash value.
Definition TString.cxx:677
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
ROOT::Math::landau_pdf
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution:
Definition PdfFuncMathCore.cxx:21
ROOT::Math::fdistribution_cdf
double fdistribution_cdf(double x, double n, double m, double x0=0)
Cumulative distribution function of the F-distribution (lower tail).
Definition ProbFuncMathCore.cxx:183
sigma
const Double_t sigma
Definition h1analysisProxy.h:11
y
Double_t y[n]
Definition legend1.C:17
c1
return c1
Definition legend1.C:41
x
Double_t x[n]
Definition legend1.C:17
n
const Int_t n
Definition legend1.C:16
f1
TF1 * f1
Definition legend1.C:11
c2
return c2
Definition legend2.C:14
c3
return c3
Definition legend3.C:15
F
#define F(x, y, z)
ROOT::Double_v
Double_t Double_v
Definition Types.h:55
TMath
TMath.
Definition TMathBase.h:35
TMath::FDistI
Double_t FDistI(Double_t F, Double_t N, Double_t M)
Calculates the cumulative distribution function of F-distribution (see ROOT::Math::fdistribution_cdf)...
Definition TMath.cxx:2297
TMath::LogNormal
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:2437
TMath::DiLog
Double_t DiLog(Double_t x)
Modified Struve functions of order 1.
Definition TMath.cxx:116
TMath::BetaDist
Double_t BetaDist(Double_t x, Double_t p, Double_t q)
Computes the probability density function of the Beta distribution (the cumulative distribution funct...
Definition TMath.cxx:2071
TMath::GamSer
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
TMath::VavilovDenEval
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:3155
TMath::ACos
Double_t ACos(Double_t)
Returns the principal value of the arc cosine of x, expressed in radians.
Definition TMath.h:632
TMath::BesselI
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
TMath::KOrdStat
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:1359
TMath::Gaus
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
TMath::Factorial
Double_t Factorial(Int_t i)
Computes factorial(n).
Definition TMath.cxx:252
TMath::KolmogorovTest
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
TMath::Nint
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition TMath.h:693
TMath::BinomialI
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:2141
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
TMath::Vavilov
Double_t Vavilov(Double_t x, Double_t kappa, Double_t beta2)
Returns the value of the Vavilov probability density function.
Definition TMath.cxx:2778
TMath::Binomial
Double_t Binomial(Int_t n, Int_t k)
Calculates the binomial coefficient n over k.
Definition TMath.cxx:2111
TMath::Normalize
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition TMath.cxx:518
TMath::Prob
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
TMath::Log2
Double_t Log2(Double_t x)
Returns the binary (base-2) logarithm of x.
Definition TMath.cxx:107
TMath::BesselK1
Double_t BesselK1(Double_t x)
Modified Bessel function I_1(x)
Definition TMath.cxx:1529
TMath::BubbleHigh
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
TMath::Exp
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:709
TMath::BesselI1
Double_t BesselI1(Double_t x)
Modified Bessel function K_0(x)
Definition TMath.cxx:1494
TMath::Erf
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition TMath.cxx:190
TMath::Permute
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:2557
TMath::QuietNaN
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition TMath.h:902
TMath::PoissonI
Double_t PoissonI(Double_t x, Double_t par)
Computes the Discrete Poisson distribution function for (x,par).
Definition TMath.cxx:615
TMath::CauchyDist
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:2175
TMath::StruveL1
Double_t StruveL1(Double_t x)
Modified Struve functions of order 0.
Definition TMath.cxx:1970
TMath::ASinH
Double_t ASinH(Double_t)
Returns the area hyperbolic sine of x.
Definition TMath.cxx:67
TMath::LaplaceDistI
Double_t LaplaceDistI(Double_t x, Double_t alpha=0, Double_t beta=1)
Computes the cumulative distribution function (lower tail integral) of Laplace distribution at point ...
Definition TMath.cxx:2380
TMath::Hash
ULong_t Hash(const void *txt, Int_t ntxt)
Calculates hash index from any char string.
Definition TMath.cxx:1408
TMath::BreitWigner
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
TMath::Sign
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:175
TMath::Landau
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
TMath::Voigt
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
TMath::Student
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:2623
TMath::PiOver2
constexpr Double_t PiOver2()
Definition TMath.h:51
TMath::BetaDistI
Double_t BetaDistI(Double_t x, Double_t p, Double_t q)
Computes the cumulative distribution function of the Beta distribution, i.e.
Definition TMath.cxx:2090
TMath::FloorNint
Int_t FloorNint(Double_t x)
Returns the nearest integer of TMath::Floor(x).
Definition TMath.h:686
TMath::ACosH
Double_t ACosH(Double_t)
Returns the nonnegative area hyperbolic cosine of x.
Definition TMath.cxx:81
TMath::BesselK0
Double_t BesselK0(Double_t x)
Modified Bessel function I_0(x)
Definition TMath.cxx:1460
TMath::BesselY0
Double_t BesselY0(Double_t x)
Bessel function J1(x) for any real x.
Definition TMath.cxx:1705
TMath::BetaCf
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
TMath::ErfInverse
Double_t ErfInverse(Double_t x)
Returns the inverse error function.
Definition TMath.cxx:208
TMath::LaplaceDist
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:2364
TMath::Log
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:756
TMath::Erfc
Double_t Erfc(Double_t x)
Computes the complementary error function erfc(x).
Definition TMath.cxx:199
TMath::VavilovI
Double_t VavilovI(Double_t x, Double_t kappa, Double_t beta2)
Returns the value of the Vavilov cumulative distribution function (lower tail integral of the probabi...
Definition TMath.cxx:2815
TMath::Beta
Double_t Beta(Double_t p, Double_t q)
Calculates Beta-function Gamma(p)*Gamma(q)/Gamma(p+q).
Definition TMath.cxx:2011
TMath::Poisson
Double_t Poisson(Double_t x, Double_t par)
Computes the Poisson distribution function for (x,par).
Definition TMath.cxx:587
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
TMath::BesselJ0
Double_t BesselJ0(Double_t x)
Modified Bessel function K_1(x)
Definition TMath.cxx:1634
TMath::Gamma
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition TMath.cxx:353
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
TMath::StruveL0
Double_t StruveL0(Double_t x)
Struve functions of order 1.
Definition TMath.cxx:1923
TMath::NormQuantile
Double_t NormQuantile(Double_t p)
Computes quantiles for standard normal distribution N(0, 1) at probability p.
Definition TMath.cxx:2456
TMath::GamCf
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
TMath::Hypot
Double_t Hypot(Double_t x, Double_t y)
Returns sqrt(x*x + y*y)
Definition TMath.cxx:59
TMath::Cos
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:594
TMath::Quantiles
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=nullptr, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities.
Definition TMath.cxx:1207
TMath::Pi
constexpr Double_t Pi()
Definition TMath.h:37
TMath::StruveH0
Double_t StruveH0(Double_t x)
Bessel function Y1(x) for positive x.
Definition TMath.cxx:1777
TMath::LnGamma
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition TMath.cxx:509
TMath::KolmogorovProb
Double_t KolmogorovProb(Double_t z)
Calculates the Kolmogorov distribution function,.
Definition TMath.cxx:679
TMath::RootsCubic
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
TMath::ChisquareQuantile
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition TMath.cxx:2193
TMath::Sin
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:588
TMath::FDist
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:2277
TMath::SignalingNaN
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition TMath.h:910
TMath::BreitWignerRelativistic
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
TMath::BubbleLow
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
TMath::BesselK
Double_t BesselK(Int_t n, Double_t x)
Integer order modified Bessel function I_n(x)
Definition TMath.cxx:1561
TMath::BesselJ1
Double_t BesselJ1(Double_t x)
Bessel function J0(x) for any real x.
Definition TMath.cxx:1669
TMath::BetaIncomplete
Double_t BetaIncomplete(Double_t x, Double_t a, Double_t b)
Calculates the incomplete Beta-function.
Definition TMath.cxx:2103
TMath::StruveH1
Double_t StruveH1(Double_t x)
Struve functions of order 0.
Definition TMath.cxx:1846
TMath::Freq
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition TMath.cxx:270
TMath::LandauI
Double_t LandauI(Double_t x)
Returns the cumulative (lower tail integral) of the Landau distribution function at point x.
Definition TMath.cxx:2845
TMath::ATanH
Double_t ATanH(Double_t)
Returns the area hyperbolic tangent of x.
Definition TMath.cxx:95
TMath::BesselI0
Double_t BesselI0(Double_t x)
Integer order modified Bessel function K_n(x)
Definition TMath.cxx:1426
TMath::VavilovSet
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:2854
TMath::Log10
Double_t Log10(Double_t x)
Returns the common (base-10) logarithm of x.
Definition TMath.h:762
TMath::StudentI
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:2646
TMath::StudentQuantile
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:2674
TMath::BesselY1
Double_t BesselY1(Double_t x)
Bessel function Y0(x) for positive x.
Definition TMath.cxx:1739
TMath::Abs
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TMath::GammaDist
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:2347
TMath::HC
constexpr Double_t HC()
in
Definition TMath.h:233
TMath::ErfcInverse
Double_t ErfcInverse(Double_t x)
Returns the inverse of the complementary error function.
Definition TMath.cxx:242
TMath::Even
Bool_t Even(Long_t a)
Returns true if a is even.
Definition TMathBase.h:113
v
@ v
Definition rootcling_impl.cxx:3687
TMath::Limits::Epsilon
static T Epsilon()
Returns minimum double representation.
Definition TMath.h:941
m
TMarker m
Definition textangle.C:8
sum
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345