CMS 3D CMS Logo
List of all members | Public Member Functions | Private Attributes
VVIObj Class Reference

#include <VVIObj.h>

Public Member Functions

double fcn (double x) const
 
void limits (double &xl, double &xu) const
 density (mode=0) or distribution (mode=1) function More...
 
 VVIObj (double kappa=0.01, double beta2=1., int mode=0)
 Constructor. More...
 

Private Attributes

double a_ [155]
 
double b_ [155]
 
const int mode_
 returns the limits on the non-zero (mode=0) or normalized region (mode=1) More...
 
double omega_
 
double t0_
 
double t1_
 
double t_
 
double 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 VVIObj.h.

Constructor & Destructor Documentation

◆ VVIObj()

VVIObj::VVIObj ( double  kappa = 0.01,
double  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 42 of file VVIObj.cc.

42  : mode_(mode) {
43  const double xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
44  const double xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
45  double h_[7];
46  double q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
47  int lp, lq, k, l, n;
48 
49  // Make sure that the inputs are reasonable
50 
51  if (kappa < 0.01)
52  kappa = 0.01;
53  if (kappa > 10.)
54  kappa = 10.;
55  if (beta2 < 0.)
56  beta2 = 0.;
57  if (beta2 > 1.)
58  beta2 = 1.;
59 
60  h_[4] = 1. - beta2 * 0.42278433999999998 + 7.6 / kappa;
61  h_[5] = beta2;
62  h_[6] = 1. - beta2;
63  h4 = -7.6 / kappa - (beta2 * .57721566 + 1);
64  h5 = log(kappa);
65  h6 = 1. / kappa;
66  t0_ = (h4 - h_[4] * h5 - (h_[4] + beta2) * (log(h_[4]) + VVIObjDetails::expint(h_[4])) + exp(-h_[4])) / h_[4];
67 
68  // Set up limits for the root search
69 
70  for (lp = 0; lp < 9; ++lp) {
71  if (kappa >= xp[lp])
72  break;
73  }
74  ll = -lp - 1.5;
75  for (lq = 0; lq < 7; ++lq) {
76  if (kappa <= xq[lq])
77  break;
78  }
79  ul = lq - 6.5;
80  auto f2 = [h_](double x) {
81  return h_[4] - x + h_[5] * (std::log(std::abs(x)) + VVIObjDetails::expint(x)) - h_[6] * std::exp(-x);
82  };
83  VVIObjDetails::dzero(ll, ul, u, rv, 1.e-5, 1000, f2);
84  q = 1. / u;
85  t1_ = h4 * q - h5 - (beta2 * q + 1) * (log((fabs(u))) + VVIObjDetails::expint(u)) + exp(-u) * q;
86  t_ = t1_ - t0_;
87  omega_ = 6.2831853000000004 / t_;
88  h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
89  if (kappa >= .07) {
90  h_[0] += 6.90775527;
91  }
92  h_[1] = beta2 * kappa;
93  h_[2] = h6 * omega_;
94  h_[3] = omega_ * 1.5707963250000001;
95  auto f1 = [h_](double x) { return h_[0] + h_[1] * std::log(h_[2] * x) - h_[3] * x; };
96  VVIObjDetails::dzero(5., 155., x0_, rv, 1.e-5, 1000, f1);
97  n = x0_ + 1.;
98  d = exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
99  a_[n - 1] = 0.;
100  if (mode_ == 0) {
101  a_[n - 1] = omega_ * .31830988654751274;
102  }
103  q = -1.;
104  q2 = 2.;
105  for (k = 1; k < n; ++k) {
106  l = n - k;
107  x = omega_ * k;
108  x1 = h6 * x;
109  VVIObjDetails::sincosint(x1, c2, c1);
110  c1 = log(x) - c1;
111  c3 = sin(x1);
112  c4 = cos(x1);
113  xf1 = kappa * (beta2 * c1 - c4) - x * c2;
114  xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
115  if (mode_ == 0) {
116  d1 = q * d * omega_ * exp(xf1);
117  a_[l - 1] = d1 * cos(xf2);
118  b_[l - 1] = -d1 * sin(xf2);
119  } else {
120  d1 = q * d * exp(xf1) / k;
121  a_[l - 1] = d1 * sin(xf2);
122  b_[l - 1] = d1 * cos(xf2);
123  a_[n - 1] += q2 * a_[l - 1];
124  }
125  q = -q;
126  q2 = -q2;
127  }
128 
129 } // VVIObj

References a_, funct::abs(), b_, alignmentValidation::c1, funct::cos(), ztail::d, d1, VVIObjDetails::dzero(), MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, VVIObjDetails::expint(), sistripvvi::VVIObjDetails::f1(), sistripvvi::VVIObjDetails::f2(), dqmdumpme::k, HLTEgPhaseIITestSequence_cff::kappa, cmsLHEtoEOSManager::l, dqm-mbProfile::log, mode_, dqmiodumpmetadata::n, omega_, submitPVResolutionJobs::q, q2, funct::sin(), VVIObjDetails::sincosint(), t0_, t1_, t_, x, x0_, and testProducerWithPsetDescEmpty_cfi::x1.

Member Function Documentation

◆ fcn()

double VVIObj::fcn ( double  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 137 of file VVIObj.cc.

137  {
138  // Local variables
139 
140  double f, u, y, a0, a1;
141  double a2 = 0.;
142  double b1, b0, b2, cof;
143  int k, n, n1;
144 
145  n = x0_;
146  if (x < t0_) {
147  f = 0.;
148  } else if (x <= t1_) {
149  y = x - t0_;
150  u = omega_ * y - 3.141592653589793;
151  cof = cos(u) * 2.;
152  a1 = 0.;
153  a0 = a_[0];
154  n1 = n + 1;
155  for (k = 2; k <= n1; ++k) {
156  a2 = a1;
157  a1 = a0;
158  a0 = a_[k - 1] + cof * a1 - a2;
159  }
160  b1 = 0.;
161  b0 = b_[0];
162  for (k = 2; k <= n; ++k) {
163  b2 = b1;
164  b1 = b0;
165  b0 = b_[k - 1] + cof * b1 - b2;
166  }
167  f = (a0 - a2) * .5 + b0 * sin(u);
168  if (mode_ != 0) {
169  f += y / t_;
170  }
171  } else {
172  f = 0.;
173  if (mode_ != 0) {
174  f = 1.;
175  }
176  }
177  return f;
178 } // fcn

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

Referenced by SiPixelTemplateSplit::PixelTempSplit().

◆ limits()

void VVIObj::limits ( double &  xl,
double &  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 186 of file VVIObj.cc.

186  {
187  xl = t0_;
188  xu = t1_;
189  return;
190 } // limits

References t0_, and t1_.

Member Data Documentation

◆ a_

double VVIObj::a_[155]
private

Definition at line 41 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

◆ b_

double VVIObj::b_[155]
private

Definition at line 42 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

◆ mode_

const int VVIObj::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 35 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

◆ omega_

double VVIObj::omega_
private

Definition at line 39 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

◆ t0_

double VVIObj::t0_
private

Definition at line 36 of file VVIObj.h.

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

◆ t1_

double VVIObj::t1_
private

Definition at line 37 of file VVIObj.h.

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

◆ t_

double VVIObj::t_
private

Definition at line 38 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

◆ x0_

double VVIObj::x0_
private

Definition at line 40 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

DDAxes::y
VVIObjDetails::expint
double expint(double x)
Private version of the sine integral.
Definition: VVIObj.cc:485
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
VVIObj::a_
double a_[155]
Definition: VVIObj.h:41
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
ALCARECOPromptCalibProdSiPixelAli0T_cff.mode
mode
Definition: ALCARECOPromptCalibProdSiPixelAli0T_cff.py:96
DDAxes::x
HLTEgPhaseIITestSequence_cff.kappa
kappa
Definition: HLTEgPhaseIITestSequence_cff.py:1343
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
VVIObj::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObj.h:35
VVIObj::t1_
double t1_
Definition: VVIObj.h:37
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
VVIObj::t_
double t_
Definition: VVIObj.h:38
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
dqmdumpme.k
k
Definition: dqmdumpme.py:60
VVIObjDetails::dzero
int dzero(double a, double b, double &x0, double &rv, double eps, int mxf, F func)
Private version of the exponential integral.
Definition: VVIObj.cc:632
q2
double q2[4]
Definition: TauolaWrapper.h:88
VVIObjDetails::sincosint
void sincosint(double x, double &sint, double &cint)
Definition: VVIObj.cc:364
VVIObj::x0_
double x0_
Definition: VVIObj.h:40
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
submitPVResolutionJobs.q
q
Definition: submitPVResolutionJobs.py:84
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
sistripvvi::VVIObjDetails::f2
double f2(double x, double const *h_)
Definition: VVIObj.cc:34
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:204
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
VVIObj::t0_
double t0_
Definition: VVIObj.h:36
ztail.d
d
Definition: ztail.py:151
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
sistripvvi::VVIObjDetails::f1
double f1(double x, double const *h_)
Private version of the exponential integral.
Definition: VVIObj.cc:33
VVIObj::b_
double b_[155]
Definition: VVIObj.h:42
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
d1
static constexpr float d1
Definition: L1EGammaCrystalsEmulatorProducer.cc:85
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
VVIObj::omega_
double omega_
Definition: VVIObj.h:39
Generated for CMSSW Reference Manual by doxygen 1.8.16