|
|
|
#include <VVIObjF.h>
Public Member Functions | |
| float | fcn (float x) const |
| void | limits (float &xl, float &xu) const |
| density (mode=0) or distribution (mode=1) function | |
| VVIObjF (float kappa=0.01, float beta2=1., int mode=0) | |
| Constructor. | |
Private Attributes | |
| float | a_ [155] |
| float | b_ [155] |
| const int | mode_ |
| returns the limits on the non-zero (mode=0) or normalized region (mode=1) | |
| float | omega_ |
| float | t0_ |
| float | t1_ |
| float | t_ |
| float | x0_ |
Port of CERNLIB routines vvidis/vviden (G116) to calculate higher quality Vavilov density and distribution functions
| VVIObjF::VVIObjF | ( | float | kappa = 0.01, |
| float | beta2 = 1., |
||
| int | mode = 0 |
||
| ) |
Constructor.
Constructor Set Vavilov parameters kappa and beta2 and define whether to calculate density fcn or distribution fcn
| kappa | - (input) Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10. (Gaussian-like)] |
| beta2 | - (input) Vavilov beta2 parameter (square of particle speed in v/c units) |
| mode | - (input) set to 0 to calculate the density function and to 1 to calculate the distribution function |
Definition at line 45 of file VVIObjF.cc.
References a_, abs, b_, trackerHits::c, alignmentValidation::c1, VVIObjDetails::dzero(), alignCSCRings::e, VVIObjDetails::expint(), f, python::connectstrParser::f1, python::connectstrParser::f2, gen::k, prof2calltree::l, mode_, n, omega_, lumiQueryAPI::q, q2, alignCSCRings::s, VVIObjDetails::sincosint(), t0_, t1_, t_, x, and x0_.
: mode_(mode) { const float xp[9] = { 9.29,2.47,.89,.36,.15,.07,.03,.02,0.0 }; const float xq[7] = { .012,.03,.08,.26,.87,3.83,11.0 }; float h_[7]; float q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv; int lp, lq, k, l, n; // Make sure that the inputs are reasonable if(kappa < 0.01f) kappa = 0.01f; if(kappa > 10.f) kappa = 10.f; if(beta2 < 0.f) beta2 = 0.f; if(beta2 > 1.f) beta2 = 1.f; float invKappa = 1.f/kappa; h_[4] = 1.f - beta2*0.42278433999999998f + (7.6f*invKappa); h_[5] = beta2; h_[6] = 1.f - beta2; h4 = - (7.6f*invKappa) - (beta2 * .57721566f + 1.f); h5 = vdt::fast_logf(kappa); h6 = invKappa; t0_ = (h4 - h_[4]*h5 - (h_[4] + beta2)*(vdt::fast_logf(h_[4]) + VVIObjFDetails::expint(h_[4])) + vdt::fast_expf(-h_[4]))/h_[4]; // Set up limits for the root search for (lp = 0; lp < 9; ++lp) { if (kappa >= xp[lp]) break; } ll = -float(lp) - 1.5f; for (lq = 0; lq < 7; ++lq) { if (kappa <= xq[lq]) break; } ul = lq - 6.5f; auto f2 = [h_](float x) { return h_[4]-x+h_[5]*(vdt::fast_logf(std::abs(x))+VVIObjFDetails::expint(x))-h_[6]*vdt::fast_expf(-x);}; VVIObjFDetails::dzero(ll, ul, u, rv, 1.e-3f, 100, f2); q = 1./u; t1_ = h4 * q - h5 - (beta2 * q + 1.f) * (vdt::fast_logf((fabs(u))) + VVIObjFDetails::expint(u)) + vdt::fast_expf(-u) * q; t_ = t1_ - t0_; omega_ = 6.2831853000000004f/t_; h_[0] = kappa * (beta2 * .57721566f + 2.f) + 9.9166128600000008f; if (kappa >= .07) {h_[0] += 6.90775527f;} h_[1] = beta2 * kappa; h_[2] = h6 * omega_; h_[3] = omega_ * 1.5707963250000001f; auto f1 = [h_](float x){ return h_[0]+h_[1]*vdt::fast_logf(h_[2]*x)-h_[3]*x;}; VVIObjFDetails::dzero(5.f, 155.f, x0_, rv, 1.e-3f, 100, f1); n = x0_ + 1.; d = vdt::fast_expf(kappa * (beta2 * (.57721566f - h5) + 1.f)) * .31830988654751274f; a_[n - 1] = 0.f; if (mode_ == 0) { a_[n - 1] = omega_ * .31830988654751274f; } q = -1.; q2 = 2.; for (k = 1; k < n; ++k) { l = n - k; x = omega_ * k; x1 = h6 * x; VVIObjFDetails::sincosint(x1,c2,c1); c1 = vdt::fast_logf(x) - c1; vdt::fast_sincosf(x1,c3,c4); xf1 = kappa * (beta2 * c1 - c4) - x * c2; xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x; float s,c; vdt::fast_sincosf(xf2,s,c); if (mode_ == 0) { d1 = q * d * omega_ * vdt::fast_expf(xf1); a_[l - 1] = d1 * c; b_[l - 1] = -d1 * s; } else { d1 = q * d * vdt::fast_expf(xf1)/k; a_[l - 1] = d1 * s; b_[l - 1] = d1 * c; a_[n - 1] += q2 * a_[l - 1]; } q = -q; q2 = -q2; } } // VVIObjF
| float VVIObjF::fcn | ( | float | x | ) | const |
Vavilov function method Returns density fcn (mode=0) or distribution fcn (mode=1)
| x | - (input) Argument of function [typically defined as (Q-mpv)/sigma] |
Definition at line 133 of file VVIObjF.cc.
References a_, b_, f, gen::k, mode_, n, omega_, t0_, t1_, t_, x0_, and detailsBasic3DVector::y.
Referenced by SiPixelTemplateReco::PixelTempReco2D().
{
// Local variables
float f, u, y, a0, a1;
float a2 = 0.;
float b1, b0, b2, cof;
int k, n, n1;
n = x0_;
if (x < t0_) {
f = 0.f;
} else if (x <= t1_) {
y = x - t0_;
u = omega_ * y - 3.141592653589793f;
float su,cu; vdt::fast_sincosf(u,su,cu);
cof = cu * 2.f;
a1 = 0.;
a0 = a_[0];
n1=n+1;
for (k = 2; k <= n1; ++k) {
a2 = a1;
a1 = a0;
a0 = a_[k - 1] + cof * a1 - a2;
}
b1 = 0.;
b0 = b_[0];
for (k = 2; k <= n; ++k) {
b2 = b1;
b1 = b0;
b0 = b_[k - 1] + cof * b1 - b2;
}
f = (a0 - a2) * .5f + b0 * su;
if (mode_ != 0) {f += y / t_;}
} else {
f = 0.f;
if (mode_ != 0) {f = 1.f;}
}
return f;
} // fcn
| void VVIObjF::limits | ( | float & | xl, |
| float & | xu | ||
| ) | const |
density (mode=0) or distribution (mode=1) function
Vavilov limits method
| xl | - (output) Smallest value of the argument for the density and the beginning of the normalized region for the distribution |
| xu | - (output) Largest value of the argument for the density and the end of the normalized region for the distribution |
Definition at line 183 of file VVIObjF.cc.
float VVIObjF::a_[155] [private] |
float VVIObjF::b_[155] [private] |
const int VVIObjF::mode_ [private] |
float VVIObjF::omega_ [private] |
float VVIObjF::t0_ [private] |
float VVIObjF::t1_ [private] |
float VVIObjF::t_ [private] |
float VVIObjF::x0_ [private] |
1.7.3