Logo ROOT  
Reference Guide
VavilovAccurate.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: B. List 29.4.2010
3 
4 
5  /**********************************************************************
6  * *
7  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
8  * *
9  * This library is free software; you can redistribute it and/or *
10  * modify it under the terms of the GNU General Public License *
11  * as published by the Free Software Foundation; either version 2 *
12  * of the License, or (at your option) any later version. *
13  * *
14  * This library is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
17  * General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this library (see file COPYING); if not, write *
21  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
22  * 330, Boston, MA 02111-1307 USA, or contact the author. *
23  * *
24  **********************************************************************/
25 
26 // Implementation file for class VavilovAccurate
27 //
28 // Created by: blist at Thu Apr 29 11:19:00 2010
29 //
30 // Last update: Thu Apr 29 11:19:00 2010
31 //
32 
33 
34 #include "Math/VavilovAccurate.h"
35 #include "Math/SpecFuncMathCore.h"
36 #include "Math/SpecFuncMathMore.h"
37 #include "Math/QuantFuncMathCore.h"
38 
39 #include <cassert>
40 #include <iostream>
41 #include <cmath>
42 #include <limits>
43 
44 
45 namespace ROOT {
46 namespace Math {
47 
48 VavilovAccurate *VavilovAccurate::fgInstance = nullptr;
49 
50 
51 VavilovAccurate::VavilovAccurate(double kappa, double beta2, double epsilonPM, double epsilon)
52 {
53  Set (kappa, beta2, epsilonPM, epsilon);
54 }
55 
56 
58 {
59  // desctructor (clean up resources)
60 }
61 
62 void VavilovAccurate::SetKappaBeta2 (double kappa, double beta2) {
63  Set (kappa, beta2);
64 }
65 
66 void VavilovAccurate::Set(double kappa, double beta2, double epsilonPM, double epsilon) {
67  // Method described in
68  // B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
69  // <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>.
70  fQuantileInit = false;
71 
72  fKappa = kappa;
73  fBeta2 = beta2;
74  fEpsilonPM = epsilonPM; // epsilon_+ = epsilon_-: determines support (T0, T1)
75  fEpsilon = epsilon;
76 
77  static const double eu = 0.577215664901532860606; // Euler's constant
78  static const double pi2 = 6.28318530717958647693, // 2pi
79  rpi = 0.318309886183790671538, // 1/pi
80  pih = 1.57079632679489661923; // pi/2
81  double h1 = -std::log(fEpsilon)-1.59631259113885503887; // -ln(fEpsilon) + ln(2/pi**2)
82  double deltaEpsilon = 0.001;
83  static const double logdeltaEpsilon = -std::log(deltaEpsilon); // 3 ln 10 = -ln(.001);
84  double logEpsilonPM = std::log(fEpsilonPM);
85  static const double eps = 1e-5; // accuracy of root finding for x0
86 
87  double xp[9] = {0,
88  9.29, 2.47, 0.89, 0.36, 0.15, 0.07, 0.03, 0.02};
89  double xq[7] = {0,
90  0.012, 0.03, 0.08, 0.26, 0.87, 3.83};
91 
92  if (kappa < 0.001) {
93  std::cerr << "VavilovAccurate::Set: kappa = " << kappa << " - out of range" << std::endl;
94  if (kappa < 0.001) kappa = 0.001;
95  }
96  if (beta2 < 0 || beta2 > 1) {
97  std::cerr << "VavilovAccurate::Set: beta2 = " << beta2 << " - out of range" << std::endl;
98  if (beta2 < 0) beta2 = -beta2;
99  if (beta2 > 1) beta2 = 1;
100  }
101 
102  // approximation of x_-
103  fH[5] = 1-beta2*(1-eu)-logEpsilonPM/kappa; // eq. 3.9
104  fH[6] = beta2;
105  fH[7] = 1-beta2;
106  double h4 = logEpsilonPM/kappa-(1+beta2*eu);
107  double logKappa = std::log(kappa);
108  double kappaInv = 1/kappa;
109  // Calculate T0 from Eq. (3.6), using x_- = fH[5]
110 // double e1h5 = (fH[5] > 40 ) ? 0 : -ROOT::Math::expint (-fH[5]);
111 // fT0 = (h4-fH[5]*logKappa-(fH[5]+beta2)*(std::log(fH[5])+e1h5)+std::exp(-fH[5]))/fH[5];
112  fT0 = (h4-fH[5]*logKappa-(fH[5]+beta2)*E1plLog(fH[5])+std::exp(-fH[5]))/fH[5];
113  int lp = 1;
114  while (lp < 9 && kappa < xp[lp]) ++lp;
115  int lq = 1;
116  while (lq < 7 && kappa >= xq[lq]) ++lq;
117  // Solve eq. 3.7 to get fH[0] = x_+
118  double delta = 0;
119  int ifail = 0;
120  do {
121  ifail = Rzero(-lp-0.5-delta,lq-7.5+delta,fH[0],eps,1000,&ROOT::Math::VavilovAccurate::G116f2);
122  delta += 0.5;
123  } while (ifail == 2);
124 
125  double q = 1/fH[0];
126  // Calculate T1 from Eq. (3.6)
127 // double e1h0 = (fH[0] > 40 ) ? 0 : -ROOT::Math::expint (-fH[0]);
128 // fT1 = h4*q-logKappa-(1+beta2*q)*(std::log(std::fabs(fH[0]))+e1h0)+std::exp(-fH[0])*q;
129  fT1 = h4*q-logKappa-(1+beta2*q)*E1plLog(fH[0])+std::exp(-fH[0])*q;
130 
131  fT = fT1-fT0; // Eq. (2.5)
132  fOmega = pi2/fT; // Eq. (2.5)
133  fH[1] = kappa*(2+beta2*eu)+h1;
134  if(kappa >= 0.07) fH[1] += logdeltaEpsilon; // reduce fEpsilon by a factor .001 for large kappa
135  fH[2] = beta2*kappa;
136  fH[3] = kappaInv*fOmega;
137  fH[4] = pih*fOmega;
138 
139  // Solve log(eq. (4.10)) to get fX0 = N
141 // if (ifail) {
142 // std::cerr << "Rzero failed for x0: ifail=" << ifail << ", kappa=" << kappa << ", beta2=" << beta2 << std::endl;
143 // std::cerr << "G116f1(" << 5. << ")=" << G116f1(5.) << ", G116f1(" << MAXTERMS << ")=" << G116f1(MAXTERMS) << std::endl;
144 // std::cerr << "fH[0]=" << fH[0] << ", fH[1]=" << fH[1] << ", fH[2]=" << fH[2] << ", fH[3]=" << fH[3] << ", fH[4]=" << fH[4] << std::endl;
145 // std::cerr << "fH[5]=" << fH[5] << ", fH[6]=" << fH[6] << ", fH[7]=" << fH[7] << std::endl;
146 // std::cerr << "fT0=" << fT0 << ", fT1=" << fT1 << std::endl;
147 // std::cerr << "G116f2(" << fH[0] << ")=" << G116f2(fH[0]) << std::endl;
148 // }
149  if (ifail == 2) {
150  fX0 = (G116f1(5) > G116f1(MAXTERMS)) ? MAXTERMS : 5;
151  }
152  if (fX0 < 5) fX0 = 5;
153  else if (fX0 > MAXTERMS) fX0 = MAXTERMS;
154  int n = int(fX0+1);
155  // logKappa=log(kappa)
156  double d = rpi*std::exp(kappa*(1+beta2*(eu-logKappa)));
157  fA_pdf[n] = rpi*fOmega;
158  fA_cdf[n] = 0;
159  q = -1;
160  double q2 = 2;
161  for (int k = 1; k < n; ++k) {
162  int l = n-k;
163  double x = fOmega*k;
164  double x1 = kappaInv*x;
165  double c1 = std::log(x)-ROOT::Math::cosint(x1);
166  double c2 = ROOT::Math::sinint(x1);
167  double c3 = std::sin(x1);
168  double c4 = std::cos(x1);
169  double xf1 = kappa*(beta2*c1-c4)-x*c2;
170  double xf2 = x*(c1 + fT0) + kappa*(c3+beta2*c2);
171  double d1 = q*d*fOmega*std::exp(xf1);
172  double s = std::sin(xf2);
173  double c = std::cos(xf2);
174  fA_pdf[l] = d1*c;
175  fB_pdf[l] = -d1*s;
176  d1 = q*d*std::exp(xf1)/k;
177  fA_cdf[l] = d1*s;
178  fB_cdf[l] = d1*c;
179  fA_cdf[n] += q2*fA_cdf[l];
180  q = -q;
181  q2 = -q2;
182  }
183 }
184 
185 void VavilovAccurate::InitQuantile() const {
186  fQuantileInit = true;
187 
188  fNQuant = 16;
189  // for kappa<0.02: use landau_quantile as first approximation
190  if (fKappa < 0.02) return;
191  else if (fKappa < 0.05) fNQuant = 32;
192 
193  // crude approximation for the median:
194 
195  double estmedian = -4.22784335098467134e-01-std::log(fKappa)-fBeta2;
196  if (estmedian>1.3) estmedian = 1.3;
197 
198  // distribute test values evenly below and above the median
199  for (int i = 1; i < fNQuant/2; ++i) {
200  double x = fT0 + i*(estmedian-fT0)/(fNQuant/2);
201  fQuant[i] = Cdf(x);
202  fLambda[i] = x;
203  }
204  for (int i = fNQuant/2; i < fNQuant-1; ++i) {
205  double x = estmedian + (i-fNQuant/2)*(fT1-estmedian)/(fNQuant/2-1);
206  fQuant[i] = Cdf(x);
207  fLambda[i] = x;
208  }
209 
210  fQuant[0] = 0;
211  fLambda[0] = fT0;
212  fQuant[fNQuant-1] = 1;
213  fLambda[fNQuant-1] = fT1;
214 
215 }
216 
217 double VavilovAccurate::Pdf (double x) const {
218 
219  static const double pi = 3.14159265358979323846; // pi
220 
221  int n = int(fX0);
222  double f;
223  if (x < fT0) {
224  f = 0;
225  } else if (x <= fT1) {
226  double y = x-fT0;
227  double u = fOmega*y-pi;
228  double cof = 2*cos(u);
229  double a1 = 0;
230  double a0 = fA_pdf[1];
231  double a2 = 0;
232  for (int k = 2; k <= n+1; ++k) {
233  a2 = a1;
234  a1 = a0;
235  a0 = fA_pdf[k]+cof*a1-a2;
236  }
237  double b1 = 0;
238  double b0 = fB_pdf[1];
239  for (int k = 2; k <= n; ++k) {
240  double b2 = b1;
241  b1 = b0;
242  b0 = fB_pdf[k]+cof*b1-b2;
243  }
244  f = 0.5*(a0-a2)+b0*sin(u);
245  } else {
246  f = 0;
247  }
248  return f;
249 }
250 
251 double VavilovAccurate::Pdf (double x, double kappa, double beta2) {
252  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
253  return Pdf (x);
254 }
255 
256 double VavilovAccurate::Cdf (double x) const {
257 
258  static const double pi = 3.14159265358979323846; // pi
259 
260  int n = int(fX0);
261  double f;
262  if (x < fT0) {
263  f = 0;
264  } else if (x <= fT1) {
265  double y = x-fT0;
266  double u = fOmega*y-pi;
267  double cof = 2*cos(u);
268  double a1 = 0;
269  double a0 = fA_cdf[1];
270  double a2 = 0;
271  for (int k = 2; k <= n+1; ++k) {
272  a2 = a1;
273  a1 = a0;
274  a0 = fA_cdf[k]+cof*a1-a2;
275  }
276  double b1 = 0;
277  double b0 = fB_cdf[1];
278  for (int k = 2; k <= n; ++k) {
279  double b2 = b1;
280  b1 = b0;
281  b0 = fB_cdf[k]+cof*b1-b2;
282  }
283  f = 0.5*(a0-a2)+b0*sin(u);
284  f += y/fT;
285  } else {
286  f = 1;
287  }
288  return f;
289 }
290 
291 double VavilovAccurate::Cdf (double x, double kappa, double beta2) {
292  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
293  return Cdf (x);
294 }
295 
296 double VavilovAccurate::Cdf_c (double x) const {
297 
298  static const double pi = 3.14159265358979323846; // pi
299 
300  int n = int(fX0);
301  double f;
302  if (x < fT0) {
303  f = 1;
304  } else if (x <= fT1) {
305  double y = fT1-x;
306  double u = fOmega*y-pi;
307  double cof = 2*cos(u);
308  double a1 = 0;
309  double a0 = fA_cdf[1];
310  double a2 = 0;
311  for (int k = 2; k <= n+1; ++k) {
312  a2 = a1;
313  a1 = a0;
314  a0 = fA_cdf[k]+cof*a1-a2;
315  }
316  double b1 = 0;
317  double b0 = fB_cdf[1];
318  for (int k = 2; k <= n; ++k) {
319  double b2 = b1;
320  b1 = b0;
321  b0 = fB_cdf[k]+cof*b1-b2;
322  }
323  f = -0.5*(a0-a2)+b0*sin(u);
324  f += y/fT;
325  } else {
326  f = 0;
327  }
328  return f;
329 }
330 
331 double VavilovAccurate::Cdf_c (double x, double kappa, double beta2) {
332  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
333  return Cdf_c (x);
334 }
335 
336 double VavilovAccurate::Quantile (double z) const {
337  if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
338 
339  if (!fQuantileInit) InitQuantile();
340 
341  double x;
342  if (fKappa < 0.02) {
344  if (x < fT0+5*fEpsilon) x = fT0+5*fEpsilon;
345  else if (x > fT1-10*fEpsilon) x = fT1-10*fEpsilon;
346  }
347  else {
348  // yes, I know what a binary search is, but linear search is faster for small n!
349  int i = 1;
350  while (z > fQuant[i]) ++i;
351  assert (i < fNQuant);
352 
353  assert (i >= 1);
354  assert (i < fNQuant);
355 
356  // initial solution
357  double f = (z-fQuant[i-1])/(fQuant[i]-fQuant[i-1]);
358  assert (f >= 0);
359  assert (f <= 1);
360  assert (fQuant[i] > fQuant[i-1]);
361 
362  x = f*fLambda[i] + (1-f)*fLambda[i-1];
363  }
364  if (fabs(x-fT0) < fEpsilon || fabs(x-fT1) < fEpsilon) return x;
365 
366  assert (x > fT0 && x < fT1);
367  double dx;
368  int n = 0;
369  do {
370  ++n;
371  double y = Cdf(x)-z;
372  double y1 = Pdf (x);
373  dx = - y/y1;
374  x = x + dx;
375  // protect against shooting beyond the support
376  if (x < fT0) x = 0.5*(fT0+x-dx);
377  else if (x > fT1) x = 0.5*(fT1+x-dx);
378  assert (x > fT0 && x < fT1);
379  } while (fabs(dx) > fEpsilon && n < 100);
380  return x;
381 }
382 
383 double VavilovAccurate::Quantile (double z, double kappa, double beta2) {
384  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
385  return Quantile (z);
386 }
387 
388 double VavilovAccurate::Quantile_c (double z) const {
389  if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
390 
391  if (!fQuantileInit) InitQuantile();
392 
393  double z1 = 1-z;
394 
395  double x;
396  if (fKappa < 0.02) {
398  if (x < fT0+5*fEpsilon) x = fT0+5*fEpsilon;
399  else if (x > fT1-10*fEpsilon) x = fT1-10*fEpsilon;
400  }
401  else {
402  // yes, I know what a binary search is, but linear search is faster for small n!
403  int i = 1;
404  while (z1 > fQuant[i]) ++i;
405  assert (i < fNQuant);
406 
407 // int i0=0, i1=fNQuant, i;
408 // for (int it = 0; it < LOG2fNQuant; ++it) {
409 // i = (i0+i1)/2;
410 // if (z > fQuant[i]) i0 = i;
411 // else i1 = i;
412 // }
413 // assert (i1-i0 == 1);
414 
415  assert (i >= 1);
416  assert (i < fNQuant);
417 
418  // initial solution
419  double f = (z1-fQuant[i-1])/(fQuant[i]-fQuant[i-1]);
420  assert (f >= 0);
421  assert (f <= 1);
422  assert (fQuant[i] > fQuant[i-1]);
423 
424  x = f*fLambda[i] + (1-f)*fLambda[i-1];
425  }
426  if (fabs(x-fT0) < fEpsilon || fabs(x-fT1) < fEpsilon) return x;
427 
428  assert (x > fT0 && x < fT1);
429  double dx;
430  int n = 0;
431  do {
432  ++n;
433  double y = Cdf_c(x)-z;
434  double y1 = -Pdf (x);
435  dx = - y/y1;
436  x = x + dx;
437  // protect against shooting beyond the support
438  if (x < fT0) x = 0.5*(fT0+x-dx);
439  else if (x > fT1) x = 0.5*(fT1+x-dx);
440  assert (x > fT0 && x < fT1);
441  } while (fabs(dx) > fEpsilon && n < 100);
442  return x;
443 }
444 
445 double VavilovAccurate::Quantile_c (double z, double kappa, double beta2) {
446  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
447  return Quantile_c (z);
448 }
449 
450 VavilovAccurate *VavilovAccurate::GetInstance() {
451  if (!fgInstance) fgInstance = new VavilovAccurate (1, 1);
452  return fgInstance;
453 }
454 
455 VavilovAccurate *VavilovAccurate::GetInstance(double kappa, double beta2) {
456  if (!fgInstance) fgInstance = new VavilovAccurate (kappa, beta2);
457  else if (kappa != fgInstance->fKappa || beta2 != fgInstance->fBeta2) fgInstance->Set (kappa, beta2);
458  return fgInstance;
459 }
460 
461 double vavilov_accurate_pdf (double x, double kappa, double beta2) {
462  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
463  return vavilov->Pdf (x);
464 }
465 
466 double vavilov_accurate_cdf_c (double x, double kappa, double beta2) {
467  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
468  return vavilov->Cdf_c (x);
469 }
470 
471 double vavilov_accurate_cdf (double x, double kappa, double beta2) {
472  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
473  return vavilov->Cdf (x);
474 }
475 
476 double vavilov_accurate_quantile (double z, double kappa, double beta2) {
477  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
478  return vavilov->Quantile (z);
479 }
480 
481 double vavilov_accurate_quantile_c (double z, double kappa, double beta2) {
482  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
483  return vavilov->Quantile_c (z);
484 }
485 
486 double VavilovAccurate::G116f1 (double x) const {
487  // fH[1] = kappa*(2+beta2*eu) -ln(fEpsilon) + ln(2/pi**2)
488  // fH[2] = beta2*kappa
489  // fH[3] = omwga/kappa
490  // fH[4] = pi/2 *fOmega
491  // log of Eq. (4.10)
492  return fH[1]+fH[2]*std::log(fH[3]*x)-fH[4]*x;
493 }
494 
495 double VavilovAccurate::G116f2 (double x) const {
496  // fH[5] = 1-beta2*(1-eu)-logEpsilonPM/kappa; // eq. 3.9
497  // fH[6] = beta2;
498  // fH[7] = 1-beta2;
499  // Eq. 3.7 of Schorr
500 // return fH[5]-x+fH[6]*(std::log(std::fabs(x))-ROOT::Math::expint (-x))-fH[7]*std::exp(-x);
501  return fH[5]-x+fH[6]*E1plLog(x)-fH[7]*std::exp(-x);
502 }
503 
504 int VavilovAccurate::Rzero (double a, double b, double& x0,
505  double eps, int mxf, double (VavilovAccurate::*f)(double)const) const {
506 
507  double xa, xb, fa, fb, r;
508 
509  if (a <= b) {
510  xa = a;
511  xb = b;
512  } else {
513  xa = b;
514  xb = a;
515  }
516  fa = (this->*f)(xa);
517  fb = (this->*f)(xb);
518 
519  if(fa*fb > 0) {
520  r = -2*(xb-xa);
521  x0 = 0;
522 // std::cerr << "VavilovAccurate::Rzero: fail=2, f(" << a << ")=" << (this->*f) (a)
523 // << ", f(" << b << ")=" << (this->*f) (b) << std::endl;
524  return 2;
525  }
526  int mc = 0;
527 
528  bool recalcF12 = true;
529  bool recalcFab = true;
530  bool fail = false;
531 
532 
533  double x1=0, x2=0, f1=0, f2=0, fx=0, ee=0;
534  do {
535  if (recalcF12) {
536  x0 = 0.5*(xa+xb);
537  r = x0-xa;
538  ee = eps*(std::abs(x0)+1);
539  if(r <= ee) break;
540  f1 = fa;
541  x1 = xa;
542  f2 = fb;
543  x2 = xb;
544  }
545  if (recalcFab) {
546  fx = (this->*f)(x0);
547  ++mc;
548  if(mc > mxf) {
549  fail = true;
550  break;
551  }
552  if(fx*fa > 0) {
553  xa = x0;
554  fa = fx;
555  } else {
556  xb = x0;
557  fb = fx;
558  }
559  }
560  recalcF12 = true;
561  recalcFab = true;
562 
563  double u1 = f1-f2;
564  double u2 = x1-x2;
565  double u3 = f2-fx;
566  double u4 = x2-x0;
567  if(u2 == 0 || u4 == 0) continue;
568  double f3 = fx;
569  double x3 = x0;
570  u1 = u1/u2;
571  u2 = u3/u4;
572  double ca = u1-u2;
573  double cb = (x1+x2)*u2-(x2+x0)*u1;
574  double cc = (x1-x0)*f1-x1*(ca*x1+cb);
575  if(ca == 0) {
576  if(cb == 0) continue;
577  x0 = -cc/cb;
578  } else {
579  u3 = cb/(2*ca);
580  u4 = u3*u3-cc/ca;
581  if(u4 < 0) continue;
582  x0 = -u3 + (x0+u3 >= 0 ? +1 : -1)*std::sqrt(u4);
583  }
584  if(x0 < xa || x0 > xb) continue;
585 
586  recalcF12 = false;
587 
588  r = std::abs(x0-x3) < std::abs(x0-x2) ? std::abs(x0-x3) : std::abs(x0-x2);
589  ee = eps*(std::abs(x0)+1);
590  if (r > ee) {
591  f1 = f2;
592  x1 = x2;
593  f2 = f3;
594  x2 = x3;
595  continue;
596  }
597 
598  recalcFab = false;
599 
600  fx = (this->*f) (x0);
601  if (fx == 0) break;
602  double xx, ff;
603  if (fx*fa < 0) {
604  xx = x0-ee;
605  if (xx <= xa) break;
606  ff = (this->*f)(xx);
607  fb = ff;
608  xb = xx;
609  } else {
610  xx = x0+ee;
611  if(xx >= xb) break;
612  ff = (this->*f)(xx);
613  fa = ff;
614  xa = xx;
615  }
616  if (fx*ff <= 0) break;
617  mc += 2;
618  if (mc > mxf) {
619  fail = true;
620  break;
621  }
622  f1 = f3;
623  x1 = x3;
624  f2 = fx;
625  x2 = x0;
626  x0 = xx;
627  fx = ff;
628  }
629  while (true);
630 
631  if (fail) {
632  r = -0.5*std::abs(xb-xa);
633  x0 = 0;
634  std::cerr << "VavilovAccurate::Rzero: fail=" << fail << ", f(" << a << ")=" << (this->*f) (a)
635  << ", f(" << b << ")=" << (this->*f) (b) << std::endl;
636  return 1;
637  }
638 
639  r = ee;
640  return 0;
641 }
642 
643 // Calculates log(|x|)+E_1(x)
644 double VavilovAccurate::E1plLog (double x) {
645  static const double eu = 0.577215664901532860606; // Euler's constant
646  double absx = std::fabs(x);
647  if (absx < 1E-4) {
648  return (x-0.25*x)*x-eu;
649  }
650  else if (x > 35) {
651  return log (x);
652  }
653  else if (x < -50) {
654  return -ROOT::Math::expint (-x);
655  }
656  return log(absx) -ROOT::Math::expint (-x);
657 }
658 
659 double VavilovAccurate::GetLambdaMin() const {
660  return fT0;
661 }
662 
663 double VavilovAccurate::GetLambdaMax() const {
664  return fT1;
665 }
666 
667 double VavilovAccurate::GetKappa() const {
668  return fKappa;
669 }
670 
671 double VavilovAccurate::GetBeta2() const {
672  return fBeta2;
673 }
674 
675 double VavilovAccurate::Mode() const {
676  double x = -4.22784335098467134e-01-std::log(fKappa)-fBeta2;
677  if (x>-0.223172) x = -0.223172;
678  double eps = 0.01;
679  double dx;
680 
681  do {
682  double p0 = Pdf (x - eps);
683  double p1 = Pdf (x);
684  double p2 = Pdf (x + eps);
685  double y1 = 0.5*(p2-p0)/eps;
686  double y2 = (p2-2*p1+p0)/(eps*eps);
687  dx = - y1/y2;
688  x += dx;
689  if (fabs(dx) < eps) eps = 0.1*fabs(dx);
690  } while (fabs(dx) > 1E-5);
691  return x;
692 }
693 
694 double VavilovAccurate::Mode(double kappa, double beta2) {
695  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
696  return Mode();
697 }
698 
699 double VavilovAccurate::GetEpsilonPM() const {
700  return fEpsilonPM;
701 }
702 
703 double VavilovAccurate::GetEpsilon() const {
704  return fEpsilon;
705 }
706 
707 double VavilovAccurate::GetNTerms() const {
708  return fX0;
709 }
710 
711 
712 
713 } // namespace Math
714 } // namespace ROOT
c
#define c(i)
Definition: RSha256.hxx:119
l
auto * l
Definition: textangle.C:4
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::VavilovAccurate::fQuant
double fQuant[kNquantMax]
Definition: VavilovAccurate.h:397
ROOT::Math::VavilovAccurate::GetLambdaMin
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
Definition: VavilovAccurate.cxx:697
e
#define e(i)
Definition: RSha256.hxx:121
c3
return c3
Definition: legend3.C:15
f
#define f(i)
Definition: RSha256.hxx:122
ROOT::Math::sinint
double sinint(double x)
Calculates the sine integral.
Definition: SpecFuncMathCore.cxx:122
ROOT::Math::VavilovAccurate::fT1
double fT1
Definition: VavilovAccurate.h:390
ROOT::Math::VavilovAccurate::Cdf_c
double Cdf_c(double x) const
Evaluate the Vavilov complementary cumulative probability density function.
Definition: VavilovAccurate.cxx:334
r
ROOT::R::TRInterface & r
Definition: Object.C:4
ROOT::Math::vavilov_accurate_quantile_c
double vavilov_accurate_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function.
Definition: VavilovAccurate.cxx:519
ROOT::Math::VavilovAccurate::InitQuantile
void InitQuantile() const
Definition: VavilovAccurate.cxx:223
exp
double exp(double)
ROOT::Math::VavilovAccurate::fH
double fH[8]
Definition: VavilovAccurate.h:390
ROOT::Math::VavilovAccurate::MAXTERMS
@ MAXTERMS
Definition: VavilovAccurate.h:389
ROOT::Math::VavilovAccurate::fEpsilonPM
double fEpsilonPM
Definition: VavilovAccurate.h:392
log
double log(double)
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:168
sin
double sin(double)
x
Double_t x[n]
Definition: legend1.C:17
cos
double cos(double)
ROOT::Math::VavilovAccurate::GetBeta2
virtual double GetBeta2() const
Return the current value of .
Definition: VavilovAccurate.cxx:709
ROOT::Math::VavilovAccurate::fT0
double fT0
Definition: VavilovAccurate.h:390
SpecFuncMathMore.h
SpecFuncMathCore.h
ROOT::Math::VavilovAccurate::Rzero
int Rzero(double a, double b, double &x0, double eps, int mxf, double(VavilovAccurate::*f)(double) const) const
Definition: VavilovAccurate.cxx:542
ROOT::Math::VavilovAccurate::fA_pdf
double fA_pdf[MAXTERMS+1]
Definition: VavilovAccurate.h:390
ROOT::Math::VavilovAccurate::fOmega
double fOmega
Definition: VavilovAccurate.h:390
ROOT::Math::vavilov_accurate_cdf_c
double vavilov_accurate_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
Definition: VavilovAccurate.cxx:504
ROOT::Math::VavilovAccurate::fQuantileInit
bool fQuantileInit
Definition: VavilovAccurate.h:394
ROOT::Math::VavilovAccurate::VavilovAccurate
VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5)
Initialize an object to calculate the Vavilov distribution.
Definition: VavilovAccurate.cxx:89
b
#define b(i)
Definition: RSha256.hxx:118
ROOT::Math::VavilovAccurate::G116f2
double G116f2(double x) const
Definition: VavilovAccurate.cxx:533
h1
TH1F * h1
Definition: legend1.C:5
ROOT::Math::VavilovAccurate::Pdf
double Pdf(double x) const
Evaluate the Vavilov probability density function.
Definition: VavilovAccurate.cxx:255
ROOT::Math::VavilovAccurate::fB_pdf
double fB_pdf[MAXTERMS+1]
Definition: VavilovAccurate.h:390
QuantFuncMathCore.h
ROOT::Math::VavilovAccurate::fLambda
double fLambda[kNquantMax]
Definition: VavilovAccurate.h:398
q
float * q
Definition: THbookFile.cxx:89
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:346
ROOT::Math::VavilovAccurate::fA_cdf
double fA_cdf[MAXTERMS+1]
Definition: VavilovAccurate.h:390
ROOT::Math::VavilovAccurate::GetInstance
static VavilovAccurate * GetInstance()
Returns a static instance of class VavilovFast.
Definition: VavilovAccurate.cxx:488
ROOT::Math::VavilovAccurate::GetKappa
virtual double GetKappa() const
Return the current value of .
Definition: VavilovAccurate.cxx:705
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
x3
static const double x3[11]
Definition: RooGaussKronrodIntegrator1D.cxx:392
ROOT::Math::VavilovAccurate::Quantile_c
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
Definition: VavilovAccurate.cxx:426
ROOT::Math::VavilovAccurate::E1plLog
static double E1plLog(double x)
Definition: VavilovAccurate.cxx:682
epsilon
REAL epsilon
Definition: triangle.c:617
ROOT::Math::expint
double expint(double x)
Calculates the exponential integral.
Definition: SpecFuncMathMore.cxx:285
a
auto * a
Definition: textangle.C:12
ROOT::Math::vavilov_accurate_quantile
double vavilov_accurate_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
Definition: VavilovAccurate.cxx:514
ROOT::Math::VavilovAccurate::fX0
double fX0
Definition: VavilovAccurate.h:390
y
Double_t y[n]
Definition: legend1.C:17
sqrt
double sqrt(double)
ROOT::Math::VavilovAccurate::fT
double fT
Definition: VavilovAccurate.h:390
TGeant4Unit::pi2
static constexpr double pi2
Definition: TGeant4SystemOfUnits.h:76
lq
int * lq
Definition: THbookFile.cxx:88
ROOT::Math::VavilovAccurate::fKappa
double fKappa
Definition: VavilovAccurate.h:391
ROOT::Math::VavilovAccurate::fB_cdf
double fB_cdf[MAXTERMS+1]
Definition: VavilovAccurate.h:390
ROOT::Math::eu
static const double eu
Definition: Vavilov.cxx:82
ROOT::Math::VavilovAccurate::fEpsilon
double fEpsilon
Definition: VavilovAccurate.h:392
f1
TF1 * f1
Definition: legend1.C:11
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:73
ROOT::Math::VavilovAccurate::~VavilovAccurate
virtual ~VavilovAccurate()
Destructor.
Definition: VavilovAccurate.cxx:95
ROOT::Math::VavilovAccurate
Class describing a Vavilov distribution.
Definition: VavilovAccurate.h:169
ROOT::Math::vavilov_accurate_pdf
double vavilov_accurate_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Definition: VavilovAccurate.cxx:499
ROOT::Math::VavilovAccurate::GetEpsilonPM
double GetEpsilonPM() const
Return the current value of .
Definition: VavilovAccurate.cxx:737
ROOT::Math::VavilovAccurate::G116f1
double G116f1(double x) const
Definition: VavilovAccurate.cxx:524
d
#define d(i)
Definition: RSha256.hxx:120
c2
return c2
Definition: legend2.C:14
ROOT::Math::VavilovAccurate::Cdf
double Cdf(double x) const
Evaluate the Vavilov cumulative probability density function.
Definition: VavilovAccurate.cxx:294
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:364
VavilovAccurate.h
ROOT::Math::VavilovAccurate::GetEpsilon
double GetEpsilon() const
Return the current value of .
Definition: VavilovAccurate.cxx:741
ROOT::Math::landau_quantile
double landau_quantile(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the lower tail of the Landau distribution (lan...
Definition: QuantFuncMathCore.cxx:189
ROOT::Math::VavilovAccurate::SetKappaBeta2
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
Definition: VavilovAccurate.cxx:100
ROOT::Math::VavilovAccurate::Quantile
double Quantile(double z) const
Evaluate the inverse of the Vavilov cumulative probability density function.
Definition: VavilovAccurate.cxx:374
ROOT::Math::VavilovAccurate::Mode
virtual double Mode() const
Return the value of where the pdf is maximal.
Definition: VavilovAccurate.cxx:713
ROOT::Math::VavilovAccurate::fNQuant
int fNQuant
Definition: VavilovAccurate.h:395
ROOT::Math::VavilovAccurate::Set
void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5)
(Re)Initialize the object
Definition: VavilovAccurate.cxx:104
ROOT::Math::VavilovAccurate::GetNTerms
double GetNTerms() const
Return the number of terms used in the series expansion.
Definition: VavilovAccurate.cxx:745
ROOT::Math::VavilovAccurate::fBeta2
double fBeta2
Definition: VavilovAccurate.h:391
ROOT
VSD Structures.
Definition: StringConv.hxx:21
TMath::E
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:102
ROOT::Math::VavilovAccurate::fgInstance
static VavilovAccurate * fgInstance
Definition: VavilovAccurate.h:402
ROOT::Math::VavilovAccurate::GetLambdaMax
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
Definition: VavilovAccurate.cxx:701
Math
ROOT::Math::vavilov_accurate_cdf
double vavilov_accurate_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
Definition: VavilovAccurate.cxx:509
int
c1
return c1
Definition: legend1.C:41
ROOT::Math::cosint
double cosint(double x)
Calculates the real part of the cosine integral Re(Ci).
Definition: SpecFuncMathCore.cxx:212