CMS 3D CMS Logo
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, float beta2, int mode)
 Constructor. More...
 
 VVIObjF (float kappa)
 Constructor with kappa only. 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 25 of file VVIObjF.h.

Constructor & Destructor Documentation

◆ VVIObjF() [1/2]

VVIObjF::VVIObjF ( float  kappa,
float  beta2,
int  mode 
)

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 43 of file VVIObjF.cc.

References a_, funct::abs(), b_, c, alignmentValidation::c1, ztail::d, d1, VVIObjFDetails::dzero(), MillePedeFileConverter_cfg::e, VVIObjFDetails::expint(), f, DeadROC_duringRun::f1, DeadROC_duringRun::f2, myMath::fast_expf(), myMath::fast_logf(), dqmMemoryStats::float, dqmdumpme::k, hgcalLayerClusters_cfi::kappa, MainPageGenerator::l, mode_, dqmiodumpmetadata::n, omega_, submitPVResolutionJobs::q, alignCSCRings::s, VVIObjFDetails::sincosint(), t0_, t1_, t_, x, x0_, and testProducerWithPsetDescEmpty_cfi::x1.

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

◆ VVIObjF() [2/2]

VVIObjF::VVIObjF ( float  kappa)

Constructor with kappa only.

Constructor Set Vavilov parameter kappa and calculate the distribution fcn

Parameters
kappa- (input) Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10. (Gaussian-like)]

Definition at line 144 of file VVIObjF.cc.

References a_, funct::abs(), b_, c, alignmentValidation::c1, ztail::d, d1, VVIObjFDetails::dzero(), MillePedeFileConverter_cfg::e, VVIObjFDetails::expint(), f, DeadROC_duringRun::f1, DeadROC_duringRun::f2, myMath::fast_expf(), myMath::fast_logf(), dqmMemoryStats::float, dqmdumpme::k, hgcalLayerClusters_cfi::kappa, MainPageGenerator::l, dqmiodumpmetadata::n, omega_, submitPVResolutionJobs::q, alignCSCRings::s, VVIObjFDetails::sincosint(), t0_, t1_, t_, x, x0_, and testProducerWithPsetDescEmpty_cfi::x1.

144  : mode_(1) {
145  const float xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
146  const float xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
147  float h_[5];
148  float q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
149  int lp, lq, k, l, n;
150 
151  // Make sure that the inputs are reasonable
152 
153  if (kappa < 0.01f)
154  kappa = 0.01f;
155  if (kappa > 10.f)
156  kappa = 10.f;
157 
158  float invKappa = 1.f / kappa;
159  h_[4] = 0.57721566f + (7.6f * invKappa);
160  h4 = -(7.6f * invKappa) - 1.57721566f;
161  h5 = vdt::fast_logf(kappa);
162  h6 = invKappa;
163  t0_ = (h4 - h_[4] * h5 - (h_[4] + 1.f) * (vdt::fast_logf(h_[4]) + VVIObjFDetails::expint(h_[4])) +
164  vdt::fast_expf(-h_[4])) /
165  h_[4];
166 
167  // Set up limits for the root search
168 
169  for (lp = 0; lp < 9; ++lp) {
170  if (kappa >= xp[lp])
171  break;
172  }
173  ll = -float(lp) - 1.5f;
174  for (lq = 0; lq < 7; ++lq) {
175  if (kappa <= xq[lq])
176  break;
177  }
178  ul = lq - 6.5f;
179  auto f2 = [h_](float x) { return h_[4] - x + (vdt::fast_logf(std::abs(x)) + VVIObjFDetails::expint(x)); };
180  VVIObjFDetails::dzero(ll, ul, u, rv, 1.e-3f, 100, f2);
181  q = 1. / u;
182  t1_ = h4 * q - h5 - (q + 1.f) * (vdt::fast_logf((fabs(u))) + VVIObjFDetails::expint(u)) + vdt::fast_expf(-u) * q;
183  t_ = t1_ - t0_;
184  omega_ = 6.2831853000000004f / t_;
185  h_[0] = kappa * 2.57721566f + 9.9166128600000008f;
186  if (kappa >= .07) {
187  h_[0] += 6.90775527f;
188  }
189  h_[1] = kappa;
190  h_[2] = h6 * omega_;
191  h_[3] = omega_ * 1.5707963250000001f;
192  auto f1 = [h_](float x) { return h_[0] + h_[1] * vdt::fast_logf(h_[2] * x) - h_[3] * x; };
193  VVIObjFDetails::dzero(5.f, 155.f, x0_, rv, 1.e-3f, 100, f1);
194  n = x0_ + 1.;
195  d = vdt::fast_expf(kappa * ((0.57721566f - h5) + 1.f)) * .31830988654751274f;
196  a_[n - 1] = 0.f;
197  q = -1.;
198  q2 = 2.;
199  for (k = 1; k < n; ++k) {
200  l = n - k;
201  x = omega_ * k;
202  x1 = h6 * x;
203  VVIObjFDetails::sincosint(x1, c2, c1);
204  c1 = vdt::fast_logf(x) - c1;
205  vdt::fast_sincosf(x1, c3, c4);
206  xf1 = kappa * (c1 - c4) - x * c2;
207  xf2 = x * c1 + kappa * (c3 + c2) + t0_ * x;
208  float s, c;
209  vdt::fast_sincosf(xf2, s, c);
210  d1 = q * d * vdt::fast_expf(xf1) / k;
211  a_[l - 1] = d1 * s;
212  b_[l - 1] = d1 * c;
213  a_[n - 1] += q2 * a_[l - 1];
214  q = -q;
215  q2 = -q2;
216  }
217 
218 } // VVIObjF with kappa only
VVIObjF::omega_
float omega_
Definition: VVIObjF.h:40
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:38
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
DeadROC_duringRun.f1
f1
Definition: DeadROC_duringRun.py:219
VVIObjF::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObjF.h:36
VVIObjFDetails::sincosint
void sincosint(float x, float &sint, float &cint)
Definition: VVIObjF.cc:285
VVIObjF::x0_
float x0_
Definition: VVIObjF.h:41
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:433
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:37
ztail.d
d
Definition: ztail.py:151
VVIObjFDetails::expint
float expint(float x)
Private version of the cosine and sine integral.
Definition: VVIObjF.cc:287
DeadROC_duringRun.f2
f2
Definition: DeadROC_duringRun.py:220
VVIObjF::b_
float b_[155]
Definition: VVIObjF.h:43
hgcalLayerClusters_cfi.kappa
kappa
Definition: hgcalLayerClusters_cfi.py:35
myMath::fast_expf
float fast_expf(float x)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:27
d1
static constexpr float d1
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
myMath::fast_logf
float fast_logf(float x)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:28
VVIObjF::t_
float t_
Definition: VVIObjF.h:39
VVIObjF::a_
float a_[155]
Definition: VVIObjF.h:42

Member Function Documentation

◆ fcn()

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 226 of file VVIObjF.cc.

References a0, testProducerWithPsetDescEmpty_cfi::a2, a_, b0, b1, b2, b_, f, dqmdumpme::k, mode_, dqmiodumpmetadata::n, omega_, t0_, t1_, t_, x, x0_, and y.

Referenced by SiPixelTemplateReco::PixelTempReco1D(), and SiPixelTemplateReco2D::PixelTempReco2D().

226  {
227  // Local variables
228 
229  float f, u, y, a0, a1;
230  float a2 = 0.;
231  float b1, b0, b2, cof;
232  int k, n, n1;
233 
234  n = x0_;
235  if (x < t0_) {
236  f = 0.f;
237  } else if (x <= t1_) {
238  y = x - t0_;
239  u = omega_ * y - 3.141592653589793f;
240  float su, cu;
241  vdt::fast_sincosf(u, su, cu);
242  cof = cu * 2.f;
243  a1 = 0.;
244  a0 = a_[0];
245  n1 = n + 1;
246  for (k = 2; k <= n1; ++k) {
247  a2 = a1;
248  a1 = a0;
249  a0 = a_[k - 1] + cof * a1 - a2;
250  }
251  b1 = 0.;
252  b0 = b_[0];
253  for (k = 2; k <= n; ++k) {
254  b2 = b1;
255  b1 = b0;
256  b0 = b_[k - 1] + cof * b1 - b2;
257  }
258  f = (a0 - a2) * .5f + b0 * su;
259  if (mode_ != 0) {
260  f += y / t_;
261  }
262  } else {
263  f = 0.f;
264  if (mode_ != 0) {
265  f = 1.f;
266  }
267  }
268  return f;
269 } // fcn
VVIObjF::omega_
float omega_
Definition: VVIObjF.h:40
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:38
VVIObjF::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObjF.h:36
VVIObjF::x0_
float x0_
Definition: VVIObjF.h:41
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:37
VVIObjF::b_
float b_[155]
Definition: VVIObjF.h:43
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:79
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:80
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:80
VVIObjF::t_
float t_
Definition: VVIObjF.h:39
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:80
VVIObjF::a_
float a_[155]
Definition: VVIObjF.h:42

◆ limits()

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 277 of file VVIObjF.cc.

References t0_, and t1_.

277  {
278  xl = t0_;
279  xu = t1_;
280  return;
281 } // limits
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:38
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:37

Member Data Documentation

◆ a_

float VVIObjF::a_[155]
private

Definition at line 42 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

◆ b_

float VVIObjF::b_[155]
private

Definition at line 43 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

◆ mode_

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 36 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

◆ omega_

float VVIObjF::omega_
private

Definition at line 40 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

◆ t0_

float VVIObjF::t0_
private

Definition at line 37 of file VVIObjF.h.

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

◆ t1_

float VVIObjF::t1_
private

Definition at line 38 of file VVIObjF.h.

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

◆ t_

float VVIObjF::t_
private

Definition at line 39 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

◆ x0_

float VVIObjF::x0_
private

Definition at line 41 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

Generated for CMSSW Reference Manual by doxygen 1.8.14