CMS 3D CMS Logo
List of all members | Public Member Functions | Private Attributes
sistripvvi::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 25 of file VVIObj.h.

Constructor & Destructor Documentation

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 52 of file VVIObj.cc.

References a_, b_, alignmentValidation::c1, funct::cos(), edmIntegrityCheck::d, sistripvvi::VVIObjDetails::dzero(), MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, sistripvvi::VVIObjDetails::expint(), sistripvvi::VVIObjDetails::f1(), sistripvvi::VVIObjDetails::f2(), gen::k, kappa, checklumidiff::l, cmsBatch::log, mode_, gen::n, omega_, lumiQueryAPI::q, q2, funct::sin(), sistripvvi::VVIObjDetails::sincosint(), t0_, t1_, t_, x0_, and globals_cff::x1.

52  : mode_(mode) {
53 
54  const double xp[9] = { 9.29,2.47,.89,.36,.15,.07,.03,.02,0.0 };
55  const double xq[7] = { .012,.03,.08,.26,.87,3.83,11.0 };
56  double h_[7];
57  double q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
58  int lp, lq, k, l, n;
59 
60  // Make sure that the inputs are reasonable
61 
62  if(kappa < 0.01) kappa = 0.01;
63  if(kappa > 10.) kappa = 10.;
64  if(beta2 < 0.) beta2 = 0.;
65  if(beta2 > 1.) beta2 = 1.;
66 
67  h_[4] = 1. - beta2*0.42278433999999998 + 7.6/kappa;
68  h_[5] = beta2;
69  h_[6] = 1. - beta2;
70  h4 = -7.6/kappa - (beta2 * .57721566 + 1);
71  h5 = log(kappa);
72  h6 = 1./kappa;
73  t0_ = (h4 - h_[4]*h5 - (h_[4] + beta2)*(log(h_[4]) + VVIObjDetails::expint(h_[4])) + exp(-h_[4]))/h_[4];
74 
75  // Set up limits for the root search
76 
77  for (lp = 0; lp < 9; ++lp) {
78  if (kappa >= xp[lp]) break;
79  }
80  ll = -lp - 1.5;
81  for (lq = 0; lq < 7; ++lq) {
82  if (kappa <= xq[lq]) break;
83  }
84  ul = lq - 6.5;
85  // double (*fp2)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f2);
86  VVIObjDetails::dzero(ll, ul, u, rv, 1.e-5, 1000, boost::bind(&VVIObjDetails::f2, _1,h_));
87  q = 1./u;
88  t1_ = h4 * q - h5 - (beta2 * q + 1) * (log((fabs(u))) + VVIObjDetails::expint(u)) + exp(-u) * q;
89  t_ = t1_ - t0_;
90  omega_ = 6.2831853000000004/t_;
91  h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
92  if (kappa >= .07) {h_[0] += 6.90775527;}
93  h_[1] = beta2 * kappa;
94  h_[2] = h6 * omega_;
95  h_[3] = omega_ * 1.5707963250000001;
96  // double (*fp1)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f1);
97  VVIObjDetails::dzero(5., 155., x0_, rv, 1.e-5, 1000, boost::bind(&VVIObjDetails::f1, _1,h_));
98  n = x0_ + 1.;
99  d = exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
100  a_[n - 1] = 0.;
101  if (mode_ == 0) {
102  a_[n - 1] = omega_ * .31830988654751274;
103  }
104  q = -1.;
105  q2 = 2.;
106  for (k = 1; k < n; ++k) {
107  l = n - k;
108  x = omega_ * k;
109  x1 = h6 * x;
110  VVIObjDetails::sincosint(x1,c2,c1);
111  c1 = log(x) - c1;
112  c3 = sin(x1);
113  c4 = cos(x1);
114  xf1 = kappa * (beta2 * c1 - c4) - x * c2;
115  xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
116  if (mode_ == 0) {
117  d1 = q * d * omega_ * exp(xf1);
118  a_[l - 1] = d1 * cos(xf2);
119  b_[l - 1] = -d1 * sin(xf2);
120  } else {
121  d1 = q * d * exp(xf1)/k;
122  a_[l - 1] = d1 * sin(xf2);
123  b_[l - 1] = d1 * cos(xf2);
124  a_[n - 1] += q2 * a_[l - 1];
125  }
126  q = -q;
127  q2 = -q2;
128  }
129 
130 } // VVIObj
sistripvvi::VVIObjDetails::expint
double expint(double x)
Private version of the sine integral.
Definition: VVIObj.cc:452
sistripvvi::VVIObj::b_
double b_[155]
Definition: VVIObj.h:45
sistripvvi::VVIObjDetails::sincosint
void sincosint(double x, double &sint, double &cint)
Definition: VVIObj.cc:349
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
sistripvvi::VVIObj::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObj.h:38
sistripvvi::VVIObjDetails::dzero
int dzero(double a, double b, double &x0, double &rv, double eps, int mxf, F func)
Definition: VVIObj.cc:574
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
q2
double q2[4]
Definition: TauolaWrapper.h:88
ALCARECOPromptCalibProdSiPixelAli0T_cff.mode
mode
Definition: ALCARECOPromptCalibProdSiPixelAli0T_cff.py:96
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
sistripvvi::VVIObjDetails::f1
double f1(double x, double const *h_)
Private version of the exponential integral.
Definition: VVIObj.cc:35
sistripvvi::VVIObj::omega_
double omega_
Definition: VVIObj.h:42
globals_cff.x1
x1
Definition: globals_cff.py:24
sistripvvi::VVIObjDetails::f2
double f2(double x, double const *h_)
Definition: VVIObj.cc:36
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:34
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
sistripvvi::VVIObj::a_
double a_[155]
Definition: VVIObj.h:44

Member Function Documentation

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 139 of file VVIObj.cc.

References a_, b_, funct::cos(), f, gen::k, mode_, gen::n, omega_, funct::sin(), t0_, t1_, t_, and x0_.

Referenced by SiStripTemplateReco::StripTempReco1D(), and SiStripTemplateSplit::StripTempSplit().

139  {
140 
141  // Local variables
142 
143  double f, u, y, a0, a1;
144  double a2 = 0.;
145  double b1, b0, b2, cof;
146  int k, n, n1;
147 
148  n = x0_;
149  if (x < t0_) {
150  f = 0.;
151  } else if (x <= t1_) {
152  y = x - t0_;
153  u = omega_ * y - 3.141592653589793;
154  cof = cos(u) * 2.;
155  a1 = 0.;
156  a0 = a_[0];
157  n1=n+1;
158  for (k = 2; k <= n1; ++k) {
159  a2 = a1;
160  a1 = a0;
161  a0 = a_[k - 1] + cof * a1 - a2;
162  }
163  b1 = 0.;
164  b0 = b_[0];
165  for (k = 2; k <= n; ++k) {
166  b2 = b1;
167  b1 = b0;
168  b0 = b_[k - 1] + cof * b1 - b2;
169  }
170  f = (a0 - a2) * .5 + b0 * sin(u);
171  if (mode_ != 0) {f += y / t_;}
172  } else {
173  f = 0.;
174  if (mode_ != 0) {f = 1.;}
175  }
176  return f;
177 } // fcn
sistripvvi::VVIObj::b_
double b_[155]
Definition: VVIObj.h:45
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
sistripvvi::VVIObj::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObj.h:38
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
sistripvvi::VVIObj::omega_
double omega_
Definition: VVIObj.h:42
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
sistripvvi::VVIObj::a_
double a_[155]
Definition: VVIObj.h:44
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 188 of file VVIObj.cc.

References t0_, and t1_.

188  {
189 
190  xl = t0_;
191  xu = t1_;
192  return;
193 } // limits

Member Data Documentation

double sistripvvi::VVIObj::a_[155]
private

Definition at line 44 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

double sistripvvi::VVIObj::b_[155]
private

Definition at line 45 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

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

Referenced by fcn(), and VVIObj().

double sistripvvi::VVIObj::omega_
private

Definition at line 42 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

double sistripvvi::VVIObj::t0_
private

Definition at line 39 of file VVIObj.h.

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

double sistripvvi::VVIObj::t1_
private

Definition at line 40 of file VVIObj.h.

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

double sistripvvi::VVIObj::t_
private

Definition at line 41 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

double sistripvvi::VVIObj::x0_
private

Definition at line 43 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

Generated for CMSSW Reference Manual by doxygen 1.8.11