CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Attributes
VVIObjF Class Reference

#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 More...
 
 VVIObjF (float kappa=0.01, float beta2=1., int mode=0)
 Constructor. More...
 

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) More...
 
float omega_
 
float t0_
 
float t1_
 
float t_
 
float x0_
 

Detailed Description

Port of CERNLIB routines vvidis/vviden (G116) to calculate higher quality Vavilov density and distribution functions

Definition at line 24 of file VVIObjF.h.

Constructor & Destructor Documentation

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

Parameters
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_, funct::abs(), b_, EnergyCorrector::c, alignmentValidation::c1, counter::c2, VVIObjFDetails::dzero(), alignCSCRings::e, VVIObjFDetails::expint(), f, python.connectstrParser::f1, python.connectstrParser::f2, myMath::fast_expf(), myMath::fast_logf(), roll_playback::k, kappa, prof2calltree::l, mode_, n, omega_, lumiQueryAPI::q, q2, alignCSCRings::s, VVIObjFDetails::sincosint(), t0_, t1_, t_, x, and x0_.

45  : mode_(mode) {
46 
47  const float xp[9] = { 9.29,2.47,.89,.36,.15,.07,.03,.02,0.0 };
48  const float xq[7] = { .012,.03,.08,.26,.87,3.83,11.0 };
49  float h_[7];
50  float q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
51  int lp, lq, k, l, n;
52 
53  // Make sure that the inputs are reasonable
54 
55  if(kappa < 0.01f) kappa = 0.01f;
56  if(kappa > 10.f) kappa = 10.f;
57  if(beta2 < 0.f) beta2 = 0.f;
58  if(beta2 > 1.f) beta2 = 1.f;
59 
60  float invKappa = 1.f/kappa;
61  h_[4] = 1.f - beta2*0.42278433999999998f + (7.6f*invKappa);
62  h_[5] = beta2;
63  h_[6] = 1.f - beta2;
64  h4 = - (7.6f*invKappa) - (beta2 * .57721566f + 1.f);
65  h5 = vdt::fast_logf(kappa);
66  h6 = invKappa;
67  t0_ = (h4 - h_[4]*h5 - (h_[4] + beta2)*(vdt::fast_logf(h_[4]) + VVIObjFDetails::expint(h_[4])) + vdt::fast_expf(-h_[4]))/h_[4];
68 
69  // Set up limits for the root search
70 
71  for (lp = 0; lp < 9; ++lp) {
72  if (kappa >= xp[lp]) break;
73  }
74  ll = -float(lp) - 1.5f;
75  for (lq = 0; lq < 7; ++lq) {
76  if (kappa <= xq[lq]) break;
77  }
78  ul = lq - 6.5f;
79  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);};
80  VVIObjFDetails::dzero(ll, ul, u, rv, 1.e-3f, 100, f2);
81  q = 1./u;
82  t1_ = h4 * q - h5 - (beta2 * q + 1.f) * (vdt::fast_logf((fabs(u))) + VVIObjFDetails::expint(u)) + vdt::fast_expf(-u) * q;
83  t_ = t1_ - t0_;
84  omega_ = 6.2831853000000004f/t_;
85  h_[0] = kappa * (beta2 * .57721566f + 2.f) + 9.9166128600000008f;
86  if (kappa >= .07) {h_[0] += 6.90775527f;}
87  h_[1] = beta2 * kappa;
88  h_[2] = h6 * omega_;
89  h_[3] = omega_ * 1.5707963250000001f;
90  auto f1 = [h_](float x){ return h_[0]+h_[1]*vdt::fast_logf(h_[2]*x)-h_[3]*x;};
91  VVIObjFDetails::dzero(5.f, 155.f, x0_, rv, 1.e-3f, 100, f1);
92  n = x0_ + 1.;
93  d = vdt::fast_expf(kappa * (beta2 * (.57721566f - h5) + 1.f)) * .31830988654751274f;
94  a_[n - 1] = 0.f;
95  if (mode_ == 0) {
96  a_[n - 1] = omega_ * .31830988654751274f;
97  }
98  q = -1.;
99  q2 = 2.;
100  for (k = 1; k < n; ++k) {
101  l = n - k;
102  x = omega_ * k;
103  x1 = h6 * x;
104  VVIObjFDetails::sincosint(x1,c2,c1);
105  c1 = vdt::fast_logf(x) - c1;
106  vdt::fast_sincosf(x1,c3,c4);
107  xf1 = kappa * (beta2 * c1 - c4) - x * c2;
108  xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
109  float s,c; vdt::fast_sincosf(xf2,s,c);
110  if (mode_ == 0) {
111  d1 = q * d * omega_ * vdt::fast_expf(xf1);
112  a_[l - 1] = d1 * c;
113  b_[l - 1] = -d1 * s;
114  } else {
115  d1 = q * d * vdt::fast_expf(xf1)/k;
116  a_[l - 1] = d1 * s;
117  b_[l - 1] = d1 * c;
118  a_[n - 1] += q2 * a_[l - 1];
119  }
120  q = -q;
121  q2 = -q2;
122  }
123 
124 } // VVIObjF
VVIObjF::omega_
float omega_
Definition: VVIObjF.h:41
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:39
alignBH_cfg.mode
list mode
Definition: alignBH_cfg.py:7
q2
double q2[4]
Definition: TauolaWrapper.h:88
counter.c2
tuple c2
Definition: counter.py:145
VVIObjF::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObjF.h:37
roll_playback.k
tuple k
Definition: roll_playback.py:143
VVIObjFDetails::sincosint
void sincosint(float x, float &sint, float &cint)
Definition: VVIObjF.cc:193
VVIObjF::x0_
float x0_
Definition: VVIObjF.h:42
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
VVIObjFDetails::dzero
int dzero(float a, float b, float &x0, float &rv, float eps, int mxf, F func)
Private version of the exponential integral.
Definition: VVIObjF.cc:321
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:38
VVIObjFDetails::expint
float expint(float x)
Private version of the cosine and sine integral.
Definition: VVIObjF.cc:198
VVIObjF::b_
float b_[155]
Definition: VVIObjF.h:44
myMath::fast_expf
float fast_expf(float x)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:27
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:72
x
Definition: DDAxes.h:10
myMath::fast_logf
float fast_logf(float x)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:28
VVIObjF::t_
float t_
Definition: VVIObjF.h:40
VVIObjF::a_
float a_[155]
Definition: VVIObjF.h:43

Member Function Documentation

float VVIObjF::fcn ( float  x) const

Vavilov function method Returns density fcn (mode=0) or distribution fcn (mode=1)

Parameters
x- (input) Argument of function [typically defined as (Q-mpv)/sigma]

Definition at line 133 of file VVIObjF.cc.

References a_, b_, f, roll_playback::k, mode_, n, omega_, t0_, t1_, t_, x0_, and detailsBasic3DVector::y.

Referenced by SiPixelTemplateReco::PixelTempReco2D().

133  {
134 
135  // Local variables
136 
137  float f, u, y, a0, a1;
138  float a2 = 0.;
139  float b1, b0, b2, cof;
140  int k, n, n1;
141 
142  n = x0_;
143  if (x < t0_) {
144  f = 0.f;
145  } else if (x <= t1_) {
146  y = x - t0_;
147  u = omega_ * y - 3.141592653589793f;
148  float su,cu; vdt::fast_sincosf(u,su,cu);
149  cof = cu * 2.f;
150  a1 = 0.;
151  a0 = a_[0];
152  n1=n+1;
153  for (k = 2; k <= n1; ++k) {
154  a2 = a1;
155  a1 = a0;
156  a0 = a_[k - 1] + cof * a1 - a2;
157  }
158  b1 = 0.;
159  b0 = b_[0];
160  for (k = 2; k <= n; ++k) {
161  b2 = b1;
162  b1 = b0;
163  b0 = b_[k - 1] + cof * b1 - b2;
164  }
165  f = (a0 - a2) * .5f + b0 * su;
166  if (mode_ != 0) {f += y / t_;}
167  } else {
168  f = 0.f;
169  if (mode_ != 0) {f = 1.f;}
170  }
171  return f;
172 } // fcn
VVIObjF::omega_
float omega_
Definition: VVIObjF.h:41
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:39
VVIObjF::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObjF.h:37
roll_playback.k
tuple k
Definition: roll_playback.py:143
VVIObjF::x0_
float x0_
Definition: VVIObjF.h:42
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:38
VVIObjF::b_
float b_[155]
Definition: VVIObjF.h:44
x
Definition: DDAxes.h:10
VVIObjF::t_
float t_
Definition: VVIObjF.h:40
VVIObjF::a_
float a_[155]
Definition: VVIObjF.h:43
void VVIObjF::limits ( float &  xl,
float &  xu 
) const

density (mode=0) or distribution (mode=1) function

Vavilov limits method

Parameters
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.

References t0_, and t1_.

183  {
184 
185  xl = t0_;
186  xu = t1_;
187  return;
188 } // limits
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:39
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:38

Member Data Documentation

float VVIObjF::a_[155]
private

Definition at line 43 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

float VVIObjF::b_[155]
private

Definition at line 44 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

const int VVIObjF::mode_
private

returns the limits on the non-zero (mode=0) or normalized region (mode=1)

set to 0 to calculate the density function and to 1 to calculate the distribution function

Definition at line 37 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

float VVIObjF::omega_
private

Definition at line 41 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

float VVIObjF::t0_
private

Definition at line 38 of file VVIObjF.h.

Referenced by fcn(), limits(), and VVIObjF().

float VVIObjF::t1_
private

Definition at line 39 of file VVIObjF.h.

Referenced by fcn(), limits(), and VVIObjF().

float VVIObjF::t_
private

Definition at line 40 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

float VVIObjF::x0_
private

Definition at line 42 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

Generated for CMSSW Reference Manual by doxygen 1.8.5