77 static const double eu = 0.577215664901532860606;
78 static const double pi2 = 6.28318530717958647693,
79 rpi = 0.318309886183790671538,
80 pih = 1.57079632679489661923;
82 double deltaEpsilon = 0.001;
83 static const double logdeltaEpsilon = -
std::log(deltaEpsilon);
85 static const double eps = 1
e-5;
88 9.29, 2.47, 0.89, 0.36, 0.15, 0.07, 0.03, 0.02};
90 0.012, 0.03, 0.08, 0.26, 0.87, 3.83};
93 std::cerr <<
"VavilovAccurate::Set: kappa = " << kappa <<
" - out of range" << std::endl;
94 if (kappa < 0.001) kappa = 0.001;
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;
103 fH[5] = 1-beta2*(1-
eu)-logEpsilonPM/kappa;
106 double h4 = logEpsilonPM/kappa-(1+beta2*
eu);
108 double kappaInv = 1/kappa;
114 while (lp < 9 && kappa < xp[lp]) ++lp;
116 while (lq < 7 && kappa >= xq[
lq]) ++
lq;
123 }
while (ifail == 2);
133 fH[1] = kappa*(2+beta2*
eu)+
h1;
134 if(kappa >= 0.07)
fH[1] += logdeltaEpsilon;
156 double d = rpi*
std::exp(kappa*(1+beta2*(
eu-logKappa)));
161 for (
int k = 1; k <
n; ++k) {
164 double x1 = kappaInv*
x;
169 double xf1 = kappa*(beta2*
c1-c4)-
x*
c2;
170 double xf2 =
x*(
c1 +
fT0) + kappa*(
c3+beta2*
c2);
190 if (
fKappa < 0.02)
return;
196 if (estmedian>1.3) estmedian = 1.3;
199 for (
int i = 1; i <
fNQuant/2; ++i) {
219 static const double pi = 3.14159265358979323846;
225 }
else if (
x <=
fT1) {
228 double cof = 2*
cos(u);
232 for (
int k = 2; k <=
n+1; ++k) {
239 for (
int k = 2; k <=
n; ++k) {
244 f = 0.5*(a0-a2)+b0*
sin(u);
258 static const double pi = 3.14159265358979323846;
264 }
else if (
x <=
fT1) {
267 double cof = 2*
cos(u);
271 for (
int k = 2; k <=
n+1; ++k) {
278 for (
int k = 2; k <=
n; ++k) {
283 f = 0.5*(a0-a2)+b0*
sin(u);
298 static const double pi = 3.14159265358979323846;
304 }
else if (
x <=
fT1) {
307 double cof = 2*
cos(u);
311 for (
int k = 2; k <=
n+1; ++k) {
318 for (
int k = 2; k <=
n; ++k) {
323 f = -0.5*(a0-a2)+b0*
sin(u);
337 if (z < 0 || z > 1)
return std::numeric_limits<double>::signaling_NaN();
350 while (z >
fQuant[i]) ++i;
389 if (z < 0 || z > 1)
return std::numeric_limits<double>::signaling_NaN();
404 while (z1 >
fQuant[i]) ++i;
463 return vavilov->
Pdf (
x);
468 return vavilov->
Cdf_c (
x);
473 return vavilov->
Cdf (
x);
507 double xa, xb, fa, fb,
r;
528 bool recalcF12 =
true;
529 bool recalcFab =
true;
533 double x1=0,
x2=0,
f1=0, f2=0, fx=0, ee=0;
567 if(u2 == 0 || u4 == 0)
continue;
573 double cb = (
x1+
x2)*u2-(
x2+x0)*u1;
574 double cc = (
x1-x0)*
f1-
x1*(ca*
x1+cb);
576 if(cb == 0)
continue;
582 x0 = -u3 + (x0+u3 >= 0 ? +1 : -1)*
std::sqrt(u4);
584 if(x0 < xa || x0 > xb)
continue;
600 fx = (this->*
f) (x0);
616 if (fx*ff <= 0)
break;
634 std::cerr <<
"VavilovAccurate::Rzero: fail=" << fail <<
", f(" <<
a <<
")=" << (this->*
f) (a)
635 <<
", f(" <<
b <<
")=" << (this->*
f) (b) << std::endl;
645 static const double eu = 0.577215664901532860606;
648 return (
x-0.25*
x)*
x-
eu;
677 if (
x>-0.223172)
x = -0.223172;
682 double p0 =
Pdf (
x - eps);
684 double p2 =
Pdf (
x + eps);
685 double y1 = 0.5*(p2-p0)/eps;
686 double y2 = (p2-2*p1+p0)/(eps*eps);
689 if (
fabs(dx) < eps) eps = 0.1*
fabs(dx);
690 }
while (
fabs(dx) > 1
E-5);
static const double x3[11]
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
Class describing a Vavilov distribution.
double GetEpsilon() const
Return the current value of .
void InitQuantile() const
double Quantile_c(double z) const override
Evaluate the inverse of the complementary Vavilov cumulative probability density function.
double Cdf(double x) const override
Evaluate the Vavilov cumulative probability density function.
static VavilovAccurate * GetInstance()
Returns a static instance of class VavilovFast.
double G116f1(double x) const
double GetBeta2() const override
Return the current value of .
int Rzero(double a, double b, double &x0, double eps, int mxf, double(VavilovAccurate::*f)(double) const) const
void SetKappaBeta2(double kappa, double beta2) override
Change and and recalculate coefficients if necessary.
static VavilovAccurate * fgInstance
double fA_cdf[MAXTERMS+1]
double fA_pdf[MAXTERMS+1]
double Quantile(double z) const override
Evaluate the inverse of the Vavilov cumulative probability density function.
double fLambda[kNquantMax]
double GetKappa() const override
Return the current value of .
static double E1plLog(double x)
double GetLambdaMax() const override
Return the maximum value of for which is nonzero in the current approximation.
~VavilovAccurate() override
Destructor.
double Pdf(double x) const override
Evaluate the Vavilov probability density function.
double GetLambdaMin() const override
Return the minimum value of for which is nonzero in the current approximation.
double GetEpsilonPM() const
Return the current value of .
double fB_cdf[MAXTERMS+1]
double G116f2(double x) const
double Cdf_c(double x) const override
Evaluate the Vavilov complementary cumulative probability density function.
double fQuant[kNquantMax]
void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5)
(Re)Initialize the object
VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5)
Initialize an object to calculate the Vavilov distribution.
double fB_pdf[MAXTERMS+1]
double GetNTerms() const
Return the number of terms used in the series expansion.
double Mode() const override
Return the value of where the pdf is maximal.
double vavilov_accurate_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
double vavilov_accurate_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cumulative probability density function.
double vavilov_accurate_cdf(double x, double kappa, double beta2)
The Vavilov cumulative probability density function.
double vavilov_accurate_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cumulative probability density function.
double landau_quantile(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the lower tail of the Landau distribution (lan...
double vavilov_accurate_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cumulative probability density function.
double expint(double x)
Calculates the exponential integral.
double sinint(double x)
Calculates the sine integral.
double cosint(double x)
Calculates the real part of the cosine integral Re(Ci).
RVec< PromoteType< T > > abs(const RVec< T > &v)
RVec< PromoteType< T > > cos(const RVec< T > &v)
RVec< PromoteType< T > > log(const RVec< T > &v)
RVec< PromoteType< T > > exp(const RVec< T > &v)
RVec< PromoteType< T > > sin(const RVec< T > &v)
Namespace for new Math classes and functions.
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double s
static constexpr double pi
static constexpr double pi2
constexpr Double_t E()
Base of natural log: .