CMS 3D CMS Logo
VVIObjF.cc
Go to the documentation of this file.
1 //
2 // VVIObjF.cc Version 2.0
3 //
4 // Port of CERNLIB G116 Functions vviden/vvidis
5 //
6 // Created by Morris Swartz on 1/14/2010.
7 // 2010 __TheJohnsHopkinsUniversity__.
8 //
9 // V1.1 - make dzero call both fcns with a switch
10 // V1.2 - remove inappriate initializers and add methods to return non-zero/normalized region
11 // V2.0 - restructuring and speed improvements by V. Innocente
12 //
13 
14 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
15 // put CMSSW location of SimpleHelix.h here
16 #include "RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h"
17 #else
18 #include "VVIObjF.h"
19 #endif
20 
21 #include <cmath>
22 #include <algorithm>
23 #include "vdt/vdtMath.h"
24 
25 namespace VVIObjFDetails {
26  void sincosint(float x, float& sint, float& cint);
27  float expint(float x);
28 
29  template <typename F>
30  int dzero(float a, float b, float& x0, float& rv, float eps, int mxf, F func);
31 } // namespace VVIObjFDetails
32 
33 // ***************************************************************************************************************************************
39 // ***************************************************************************************************************************************
40 
41 VVIObjF::VVIObjF(float kappa, float beta2, int mode) : mode_(mode) {
42  const float xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
43  const float xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
44  float h_[7];
45  float q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
46  int lp, lq, k, l, n;
47 
48  // Make sure that the inputs are reasonable
49 
50  if (kappa < 0.01f)
51  kappa = 0.01f;
52  if (kappa > 10.f)
53  kappa = 10.f;
54  if (beta2 < 0.f)
55  beta2 = 0.f;
56  if (beta2 > 1.f)
57  beta2 = 1.f;
58 
59  float invKappa = 1.f / kappa;
60  h_[4] = 1.f - beta2 * 0.42278433999999998f + (7.6f * invKappa);
61  h_[5] = beta2;
62  h_[6] = 1.f - beta2;
63  h4 = -(7.6f * invKappa) - (beta2 * .57721566f + 1.f);
64  h5 = vdt::fast_logf(kappa);
65  h6 = invKappa;
66  t0_ = (h4 - h_[4] * h5 - (h_[4] + beta2) * (vdt::fast_logf(h_[4]) + VVIObjFDetails::expint(h_[4])) +
67  vdt::fast_expf(-h_[4])) /
68  h_[4];
69 
70  // Set up limits for the root search
71 
72  for (lp = 0; lp < 9; ++lp) {
73  if (kappa >= xp[lp])
74  break;
75  }
76  ll = -float(lp) - 1.5f;
77  for (lq = 0; lq < 7; ++lq) {
78  if (kappa <= xq[lq])
79  break;
80  }
81  ul = lq - 6.5f;
82  auto f2 = [h_](float x) {
83  return h_[4] - x + h_[5] * (vdt::fast_logf(std::abs(x)) + VVIObjFDetails::expint(x)) - h_[6] * vdt::fast_expf(-x);
84  };
85  VVIObjFDetails::dzero(ll, ul, u, rv, 1.e-3f, 100, f2);
86  q = 1. / u;
87  t1_ = h4 * q - h5 - (beta2 * q + 1.f) * (vdt::fast_logf((fabs(u))) + VVIObjFDetails::expint(u)) +
88  vdt::fast_expf(-u) * q;
89  t_ = t1_ - t0_;
90  omega_ = 6.2831853000000004f / t_;
91  h_[0] = kappa * (beta2 * .57721566f + 2.f) + 9.9166128600000008f;
92  if (kappa >= .07) {
93  h_[0] += 6.90775527f;
94  }
95  h_[1] = beta2 * kappa;
96  h_[2] = h6 * omega_;
97  h_[3] = omega_ * 1.5707963250000001f;
98  auto f1 = [h_](float x) { return h_[0] + h_[1] * vdt::fast_logf(h_[2] * x) - h_[3] * x; };
99  VVIObjFDetails::dzero(5.f, 155.f, x0_, rv, 1.e-3f, 100, f1);
100  n = x0_ + 1.;
101  d = vdt::fast_expf(kappa * (beta2 * (.57721566f - h5) + 1.f)) * .31830988654751274f;
102  a_[n - 1] = 0.f;
103  if (mode_ == 0) {
104  a_[n - 1] = omega_ * .31830988654751274f;
105  }
106  q = -1.;
107  q2 = 2.;
108  for (k = 1; k < n; ++k) {
109  l = n - k;
110  x = omega_ * k;
111  x1 = h6 * x;
112  VVIObjFDetails::sincosint(x1, c2, c1);
113  c1 = vdt::fast_logf(x) - c1;
114  vdt::fast_sincosf(x1, c3, c4);
115  xf1 = kappa * (beta2 * c1 - c4) - x * c2;
116  xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
117  float s, c;
118  vdt::fast_sincosf(xf2, s, c);
119  if (mode_ == 0) {
120  d1 = q * d * omega_ * vdt::fast_expf(xf1);
121  a_[l - 1] = d1 * c;
122  b_[l - 1] = -d1 * s;
123  } else {
124  d1 = q * d * vdt::fast_expf(xf1) / k;
125  a_[l - 1] = d1 * s;
126  b_[l - 1] = d1 * c;
127  a_[n - 1] += q2 * a_[l - 1];
128  }
129  q = -q;
130  q2 = -q2;
131  }
132 
133 } // VVIObjF
134 
135 // *************************************************************************************************************************************
139 // *************************************************************************************************************************************
140 
141 float VVIObjF::fcn(float x) const {
142  // Local variables
143 
144  float f, u, y, a0, a1;
145  float a2 = 0.;
146  float b1, b0, b2, cof;
147  int k, n, n1;
148 
149  n = x0_;
150  if (x < t0_) {
151  f = 0.f;
152  } else if (x <= t1_) {
153  y = x - t0_;
154  u = omega_ * y - 3.141592653589793f;
155  float su, cu;
156  vdt::fast_sincosf(u, su, cu);
157  cof = cu * 2.f;
158  a1 = 0.;
159  a0 = a_[0];
160  n1 = n + 1;
161  for (k = 2; k <= n1; ++k) {
162  a2 = a1;
163  a1 = a0;
164  a0 = a_[k - 1] + cof * a1 - a2;
165  }
166  b1 = 0.;
167  b0 = b_[0];
168  for (k = 2; k <= n; ++k) {
169  b2 = b1;
170  b1 = b0;
171  b0 = b_[k - 1] + cof * b1 - b2;
172  }
173  f = (a0 - a2) * .5f + b0 * su;
174  if (mode_ != 0) {
175  f += y / t_;
176  }
177  } else {
178  f = 0.f;
179  if (mode_ != 0) {
180  f = 1.f;
181  }
182  }
183  return f;
184 } // fcn
185 
186 // *************************************************************************************************************************************
190 // *************************************************************************************************************************************
191 
192 void VVIObjF::limits(float& xl, float& xu) const {
193  xl = t0_;
194  xu = t1_;
195  return;
196 } // limits
197 
198 #include "sicif.h"
199 namespace VVIObjFDetails {
200  void sincosint(float x, float& sint, float& cint) { sicif(x, sint, cint); }
201 
202  float expint(float x) {
203  // Initialized data
204 
205  const float zero = 0.;
206  const float q2[7] = {
207  .10340013040487, 3.319092135933, 20.449478501379, 41.280784189142, 32.426421069514, 10.041164382905, 1.};
208  const float p3[6] = {
209  -2.3909964453136, -147.98219500504, -254.3763397689, -119.55761038372, -19.630408535939, -.9999999999036};
210  const float q3[6] = {177.60070940351, 530.68509610812, 462.23027156148, 156.81843364539, 21.630408494238, 1.};
211  const float p4[8] = {-8.6693733995107,
212  -549.14226552109,
213  -4210.0161535707,
214  -249301.39345865,
215  -119623.66934925,
216  -22174462.775885,
217  3892804.213112,
218  -391546073.8091};
219  const float q4[8] = {34.171875,
220  -1607.0892658722,
221  35730.029805851,
222  -483547.43616216,
223  4285596.2461175,
224  -24903337.574054,
225  89192576.757561,
226  -165254299.72521};
227  const float a1[8] = {-2.1808638152072,
228  -21.901023385488,
229  9.3081638566217,
230  25.076281129356,
231  -33.184253199722,
232  60.121799083008,
233  -43.253113287813,
234  1.0044310922808};
235  const float b1[8] = {0.,
236  3.9370770185272,
237  300.89264837292,
238  -6.2504116167188,
239  1003.6743951673,
240  14.325673812194,
241  2736.2411988933,
242  .52746885196291};
243  const float a2[8] = {-3.4833465360285,
244  -18.65454548834,
245  -8.2856199414064,
246  -32.34673303054,
247  17.960168876925,
248  1.7565631546961,
249  -1.9502232128966,
250  .99999429607471};
251  const float b2[8] = {0.,
252  69.500065588743,
253  57.283719383732,
254  25.777638423844,
255  760.76114800773,
256  28.951672792514,
257  -3.4394226689987,
258  1.0008386740264};
259  const float a3[6] = {
260  -27.780928934438, -10.10479081576, -9.1483008216736, -5.0223317461851, -3.0000077799358, 1.0000000000704};
261  const float one = 1.;
262  const float b3[6] = {0., 122.39993926823, 2.7276100778779, -7.1897518395045, -2.9990118065262, 1.999999942826};
263  const float two = 2.;
264  const float three = 3.;
265  const float x0 = .37250741078137;
266  const float xl[6] = {-24., -12., -6., 0., 1., 4.};
267  const float p1[5] = {4.293125234321, 39.894153870321, 292.52518866921, 425.69682638592, -434.98143832952};
268  const float q1[5] = {1., 18.899288395003, 150.95038744251, 568.05252718987, 753.58564359843};
269  const float p2[7] = {.43096783946939,
270  6.9052252278444,
271  23.019255939133,
272  24.378408879132,
273  9.0416155694633,
274  .99997957705159,
275  4.656271079751e-7};
276 
277  // Local variables
278  float v, y, ap, bp, aq, dp, bq, dq;
279 
280  if (x <= xl[0]) {
281  ap = a3[0] - x;
282  for (int i__ = 2; i__ <= 5; ++i__) {
283  ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
284  }
285  y = vdt::fast_expf(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
286  } else if (x <= xl[1]) {
287  ap = a2[0] - x;
288  for (int i__ = 2; i__ <= 7; ++i__) {
289  ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
290  }
291  y = vdt::fast_expf(-x) / x * (a2[7] + b2[7] / ap);
292  } else if (x <= xl[2]) {
293  ap = a1[0] - x;
294  for (int i__ = 2; i__ <= 7; ++i__) {
295  ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
296  }
297  y = vdt::fast_expf(-x) / x * (a1[7] + b1[7] / ap);
298  } else if (x < xl[3]) {
299  v = -two * (x / three + one);
300  bp = zero;
301  dp = p4[0];
302  for (int i__ = 2; i__ <= 8; ++i__) {
303  ap = bp;
304  bp = dp;
305  dp = p4[i__ - 1] - ap + v * bp;
306  }
307  bq = zero;
308  dq = q4[0];
309  for (int i__ = 2; i__ <= 8; ++i__) {
310  aq = bq;
311  bq = dq;
312  dq = q4[i__ - 1] - aq + v * bq;
313  }
314  y = -vdt::fast_logf(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
315  } else if (x == xl[3]) {
316  return zero;
317  } else if (x < xl[4]) {
318  ap = p1[0];
319  aq = q1[0];
320  for (int i__ = 2; i__ <= 5; ++i__) {
321  ap = p1[i__ - 1] + x * ap;
322  aq = q1[i__ - 1] + x * aq;
323  }
324  y = -vdt::fast_logf(x) + ap / aq;
325  } else if (x <= xl[5]) {
326  y = one / x;
327  ap = p2[0];
328  aq = q2[0];
329  for (int i__ = 2; i__ <= 7; ++i__) {
330  ap = p2[i__ - 1] + y * ap;
331  aq = q2[i__ - 1] + y * aq;
332  }
333  y = vdt::fast_expf(-x) * ap / aq;
334  } else {
335  y = one / x;
336  ap = p3[0];
337  aq = q3[0];
338  for (int i__ = 2; i__ <= 6; ++i__) {
339  ap = p3[i__ - 1] + y * ap;
340  aq = q3[i__ - 1] + y * aq;
341  }
342  y = vdt::fast_expf(-x) * y * (one + y * ap / aq);
343  }
344  return y;
345  } // expint
346 
347  template <typename F>
348  int dzero(float a, float b, float& x0, float& rv, float eps, int mxf, F func) {
349  /* System generated locals */
350  float d__1, d__2, d__3, d__4;
351 
352  // Local variables
353  float f1, f2, f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc, fa, fb, ee, ff;
354  int mc;
355  float xa, xb, fx, xx, su4;
356 
357  xa = std::min(a, b);
358  xb = std::max(a, b);
359  fa = func(xa);
360  fb = func(xb);
361  if (fa * fb > 0.f) {
362  rv = (xb - xa) * -2.f;
363  x0 = 0.f;
364  return 1;
365  }
366  mc = 0;
367  L1:
368  x0 = (xa + xb) * 0.5f;
369  rv = x0 - xa;
370  ee = eps * (std::abs(x0) + 1.f);
371  if (rv <= ee) {
372  rv = ee;
373  ff = func(x0);
374  return 0;
375  }
376  f1 = fa;
377  x1 = xa;
378  f2 = fb;
379  x2 = xb;
380  L2:
381  fx = func(x0);
382  ++mc;
383  if (mc > mxf) {
384  rv = (d__1 = xb - xa, fabs(d__1)) * -0.5f;
385  x0 = 0.;
386  return 0;
387  }
388  if (fx * fa > 0.f) {
389  xa = x0;
390  fa = fx;
391  } else {
392  xb = x0;
393  fb = fx;
394  }
395  L3:
396  u1 = f1 - f2;
397  u2 = x1 - x2;
398  u3 = f2 - fx;
399  u4 = x2 - x0;
400  if (u2 == 0.f || u4 == 0.f) {
401  goto L1;
402  }
403  f3 = fx;
404  x3 = x0;
405  u1 /= u2;
406  u2 = u3 / u4;
407  ca = u1 - u2;
408  cb = (x1 + x2) * u2 - (x2 + x0) * u1;
409  cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
410  if (ca == 0.f) {
411  if (cb == 0.f) {
412  goto L1;
413  }
414  x0 = -cc / cb;
415  } else {
416  u3 = cb / (ca * 2.f);
417  u4 = u3 * u3 - cc / ca;
418  if (u4 < 0.f) {
419  goto L1;
420  }
421  su4 = std::abs(u4);
422  if (x0 + u3 < 0.f) {
423  su4 = -su4;
424  }
425  x0 = -u3 + su4;
426  }
427  if (x0 < xa || x0 > xb) {
428  goto L1;
429  }
430  // Computing MIN
431  d__3 = (d__1 = x0 - x3, std::abs(d__1));
432  d__4 = (d__2 = x0 - x2, std::abs(d__2));
433  rv = std::min(d__3, d__4);
434  ee = eps * (std::abs(x0) + 1);
435  if (rv > ee) {
436  f1 = f2;
437  x1 = x2;
438  f2 = f3;
439  x2 = x3;
440  goto L2;
441  }
442  fx = func(x0);
443  if (fx == 0.f) {
444  rv = ee;
445  ff = func(x0);
446  return 0;
447  }
448  if (fx * fa < 0.f) {
449  xx = x0 - ee;
450  if (xx <= xa) {
451  rv = ee;
452  ff = func(x0);
453  return 0;
454  }
455  ff = func(xx);
456  fb = ff;
457  xb = xx;
458  } else {
459  xx = x0 + ee;
460  if (xx >= xb) {
461  rv = ee;
462  ff = func(x0);
463  return 0;
464  }
465  ff = func(xx);
466  fa = ff;
467  xa = xx;
468  }
469  if (fx * ff > 0.f) {
470  mc += 2;
471  if (mc > mxf) {
472  rv = (d__1 = xb - xa, std::abs(d__1)) * -0.5f;
473  x0 = 0.f;
474  return 0;
475  }
476  f1 = f3;
477  x1 = x3;
478  f2 = fx;
479  x2 = x0;
480  x0 = xx;
481  fx = ff;
482  goto L3;
483  }
484  /* L4: */
485  rv = ee;
486  ff = func(x0);
487  return 0;
488  } // dzero
489 
490 } // namespace VVIObjFDetails
PixelRegions::L3
Definition: PixelRegionContainers.h:32
PixelRegions::L1
Definition: PixelRegionContainers.h:32
DDAxes::y
PixelRegions::L2
Definition: PixelRegionContainers.h:32
SiPixelPI::one
Definition: SiPixelPayloadInspectorHelper.h:39
VVIObjF::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObjF.h:35
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
VVIObjFDetails::sincosint
void sincosint(float x, float &sint, float &cint)
Definition: VVIObjF.cc:200
VVIObjFDetails
Definition: VVIObjF.cc:25
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
CaloTowersParam_cfi.mc
mc
Definition: CaloTowersParam_cfi.py:8
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:36
min
T min(T a, T b)
Definition: MathUtil.h:58
multPhiCorr_741_25nsDY_cfi.fx
fx
Definition: multPhiCorr_741_25nsDY_cfi.py:9
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
ALCARECOPromptCalibProdSiPixelAli0T_cff.mode
mode
Definition: ALCARECOPromptCalibProdSiPixelAli0T_cff.py:96
VVIObjF::omega_
float omega_
Definition: VVIObjF.h:39
VVIObjF::x0_
float x0_
Definition: VVIObjF.h:40
DDAxes::x
SiPixelPI::zero
Definition: SiPixelPayloadInspectorHelper.h:39
findQualityFiles.v
v
Definition: findQualityFiles.py:179
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
VVIObjFDetails::expint
float expint(float x)
Private version of the cosine and sine integral.
Definition: VVIObjF.cc:202
VVIObjF.h
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
F
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163
VVIObjF::a_
float a_[155]
Definition: VVIObjF.h:41
b1
static constexpr float b1
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
alignCSCRings.s
s
Definition: alignCSCRings.py:92
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:37
myMath::fast_expf
float fast_expf(float x)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:27
Calorimetry_cff.dp
dp
Definition: Calorimetry_cff.py:157
alignCSCRings.ff
ff
Definition: alignCSCRings.py:148
myMath::fast_logf
float fast_logf(float x)
Definition: EcalUncalibRecHitRatioMethodAlgo.h:28
p2
double p2[4]
Definition: TauolaWrapper.h:90
dqmdumpme.k
k
Definition: dqmdumpme.py:60
b
double b
Definition: hdecay.h:118
q2
double q2[4]
Definition: TauolaWrapper.h:88
q1
double q1[4]
Definition: TauolaWrapper.h:87
a
double a
Definition: hdecay.h:119
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
DeadROC_duringRun.f2
f2
Definition: DeadROC_duringRun.py:220
testProducerWithPsetDescEmpty_cfi.u3
u3
Definition: testProducerWithPsetDescEmpty_cfi.py:50
testProducerWithPsetDescEmpty_cfi.u1
u1
Definition: testProducerWithPsetDescEmpty_cfi.py:49
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
p4
double p4[4]
Definition: TauolaWrapper.h:92
p1
double p1[4]
Definition: TauolaWrapper.h:89
submitPVResolutionJobs.q
q
Definition: submitPVResolutionJobs.py:84
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
VVIObjF::fcn
float fcn(float x) const
Definition: VVIObjF.cc:141
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
TrackCollections2monitor_cff.func
func
Definition: TrackCollections2monitor_cff.py:359
cc
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:204
sicif
int sicif(float xx, float &si, float &ci)
Definition: sicif.h:219
SiPixelPI::two
Definition: SiPixelPayloadInspectorHelper.h:39
testProducerWithPsetDescEmpty_cfi.b3
b3
Definition: testProducerWithPsetDescEmpty_cfi.py:36
VVIObjF::b_
float b_[155]
Definition: VVIObjF.h:42
genVertex_cff.x
x
Definition: genVertex_cff.py:12
p3
double p3[4]
Definition: TauolaWrapper.h:91
benchmark_cfg.fa
fa
Definition: benchmark_cfg.py:13
VVIObjF::t_
float t_
Definition: VVIObjF.h:38
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:348
ztail.d
d
Definition: ztail.py:151
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
kappa
static const G4double kappa
Definition: UrbanMscModel93.cc:35
MetAnalyzer.u2
u2
Definition: MetAnalyzer.py:61
sicif.h
DeadROC_duringRun.f1
f1
Definition: DeadROC_duringRun.py:219
VVIObjF::limits
void limits(float &xl, float &xu) const
density (mode=0) or distribution (mode=1) function
Definition: VVIObjF.cc:192
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:81
benchmark_cfg.fb
fb
Definition: benchmark_cfg.py:14
d1
static constexpr float d1
Definition: L1EGammaCrystalsEmulatorProducer.cc:84
geometryCSVtoXML.xx
xx
Definition: geometryCSVtoXML.py:19
VVIObjF::VVIObjF
VVIObjF(float kappa=0.01, float beta2=1., int mode=0)
Constructor.
Definition: VVIObjF.cc:41
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16