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