CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // Copyright 2010 __TheJohnsHopkinsUniversity__. All rights reserved.
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 
22 #include <cmath>
23 #include <algorithm>
24 #include "vdt/vdtMath.h"
25 
26 namespace VVIObjFDetails {
27  void sincosint(float x, float & sint, float & cint);
28  float expint(float x);
29 
30  template<typename F>
31  int dzero(float a, float b, float& x0,
32  float& rv, float eps, int mxf, F func);
33 }
34 
35 
36 
37 // ***************************************************************************************************************************************
43 // ***************************************************************************************************************************************
44 
45 VVIObjF::VVIObjF(float kappa, float beta2, int mode) : 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
125 
126 // *************************************************************************************************************************************
130 // *************************************************************************************************************************************
131 
132 
133 float VVIObjF::fcn(float x) const {
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
173 
174 
175 
176 // *************************************************************************************************************************************
180 // *************************************************************************************************************************************
181 
182 
183 void VVIObjF::limits(float& xl, float& xu) const {
184 
185  xl = t0_;
186  xu = t1_;
187  return;
188 } // limits
189 
190 
191 #include "sicif.h"
192 namespace VVIObjFDetails {
193  void sincosint(float x, float & sint, float & cint) {
194  sicif(x,sint,cint);
195  }
196 
197 
198 float expint(float x) {
199 
200  // Initialized data
201 
202  const float zero = 0.;
203  const float q2[7] = { .10340013040487,3.319092135933,
204  20.449478501379,41.280784189142,32.426421069514,10.041164382905,
205  1. };
206  const float p3[6] = { -2.3909964453136,-147.98219500504,
207  -254.3763397689,-119.55761038372,-19.630408535939,-.9999999999036
208  };
209  const float q3[6] = { 177.60070940351,530.68509610812,
210  462.23027156148,156.81843364539,21.630408494238,1. };
211  const float p4[8] = { -8.6693733995107,-549.14226552109,
212  -4210.0161535707,-249301.39345865,-119623.66934925,
213  -22174462.775885,3892804.213112,-391546073.8091 };
214  const float q4[8] = { 34.171875,-1607.0892658722,35730.029805851,
215  -483547.43616216,4285596.2461175,-24903337.574054,89192576.757561,
216  -165254299.72521 };
217  const float a1[8] = { -2.1808638152072,-21.901023385488,
218  9.3081638566217,25.076281129356,-33.184253199722,60.121799083008,
219  -43.253113287813,1.0044310922808 };
220  const float b1[8] = { 0.,3.9370770185272,300.89264837292,
221  -6.2504116167188,1003.6743951673,14.325673812194,2736.2411988933,
222  .52746885196291 };
223  const float a2[8] = { -3.4833465360285,-18.65454548834,
224  -8.2856199414064,-32.34673303054,17.960168876925,1.7565631546961,
225  -1.9502232128966,.99999429607471 };
226  const float b2[8] = { 0.,69.500065588743,57.283719383732,
227  25.777638423844,760.76114800773,28.951672792514,-3.4394226689987,
228  1.0008386740264 };
229  const float a3[6] = { -27.780928934438,-10.10479081576,
230  -9.1483008216736,-5.0223317461851,-3.0000077799358,
231  1.0000000000704 };
232  const float one = 1.;
233  const float b3[6] = { 0.,122.39993926823,2.7276100778779,
234  -7.1897518395045,-2.9990118065262,1.999999942826 };
235  const float two = 2.;
236  const float three = 3.;
237  const float x0 = .37250741078137;
238  const float xl[6] = { -24.,-12.,-6.,0.,1.,4. };
239  const float p1[5] = { 4.293125234321,39.894153870321,
240  292.52518866921,425.69682638592,-434.98143832952 };
241  const float q1[5] = { 1.,18.899288395003,150.95038744251,
242  568.05252718987,753.58564359843 };
243  const float p2[7] = { .43096783946939,6.9052252278444,
244  23.019255939133,24.378408879132,9.0416155694633,.99997957705159,
245  4.656271079751e-7 };
246 
247  // Local variables
248  float v, y, ap, bp, aq, dp, bq, dq;
249 
250  if (x <= xl[0]) {
251  ap = a3[0] - x;
252  for ( int i__ = 2; i__ <= 5; ++i__) {
253  ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
254  }
255  y = vdt::fast_expf(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
256  } else if (x <= xl[1]) {
257  ap = a2[0] - x;
258  for ( int i__ = 2; i__ <= 7; ++i__) {
259  ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
260  }
261  y = vdt::fast_expf(-x) / x * (a2[7] + b2[7] / ap);
262  } else if (x <= xl[2]) {
263  ap = a1[0] - x;
264  for ( int i__ = 2; i__ <= 7; ++i__) {
265  ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
266  }
267  y = vdt::fast_expf(-x) / x * (a1[7] + b1[7] / ap);
268  } else if (x < xl[3]) {
269  v = -two * (x / three + one);
270  bp = zero;
271  dp = p4[0];
272  for ( int i__ = 2; i__ <= 8; ++i__) {
273  ap = bp;
274  bp = dp;
275  dp = p4[i__ - 1] - ap + v * bp;
276  }
277  bq = zero;
278  dq = q4[0];
279  for ( int i__ = 2; i__ <= 8; ++i__) {
280  aq = bq;
281  bq = dq;
282  dq = q4[i__ - 1] - aq + v * bq;
283  }
284  y = -vdt::fast_logf(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
285  } else if (x == xl[3]) {
286  return zero;
287  } else if (x < xl[4]) {
288  ap = p1[0];
289  aq = q1[0];
290  for ( int i__ = 2; i__ <= 5; ++i__) {
291  ap = p1[i__ - 1] + x * ap;
292  aq = q1[i__ - 1] + x * aq;
293  }
294  y = -vdt::fast_logf(x) + ap / aq;
295  } else if (x <= xl[5]) {
296  y = one / x;
297  ap = p2[0];
298  aq = q2[0];
299  for ( int i__ = 2; i__ <= 7; ++i__) {
300  ap = p2[i__ - 1] + y * ap;
301  aq = q2[i__ - 1] + y * aq;
302  }
303  y = vdt::fast_expf(-x) * ap / aq;
304  } else {
305  y = one / x;
306  ap = p3[0];
307  aq = q3[0];
308  for ( int i__ = 2; i__ <= 6; ++i__) {
309  ap = p3[i__ - 1] + y * ap;
310  aq = q3[i__ - 1] + y * aq;
311  }
312  y = vdt::fast_expf(-x) * y * (one + y * ap / aq);
313  }
314  return y;
315 } // expint
316 
317 
318 
319 
320  template<typename F>
321  int dzero(float a, float b, float& x0,
322  float& rv, float eps, int mxf, F func) {
323  /* System generated locals */
324  float d__1, d__2, d__3, d__4;
325 
326  // Local variables
327  float f1, f2, f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc, fa, fb, ee, ff;
328  int mc;
329  float xa, xb, fx, xx, su4;
330 
331  xa = std::min(a,b);
332  xb = std::max(a,b);
333  fa = func(xa);
334  fb = func(xb);
335  if (fa * fb > 0.f) {
336  rv = (xb - xa) * -2.f;
337  x0 = 0.f;
338  return 1;
339  }
340  mc = 0;
341  L1:
342  x0 = (xa + xb) * 0.5f;
343  rv = x0 - xa;
344  ee = eps * (std::abs(x0) + 1.f);
345  if (rv <= ee) {
346  rv = ee;
347  ff = func(x0);
348  return 0;
349  }
350  f1 = fa;
351  x1 = xa;
352  f2 = fb;
353  x2 = xb;
354  L2:
355  fx = func(x0);
356  ++mc;
357  if (mc > mxf) {
358  rv = (d__1 = xb - xa, fabs(d__1)) * -0.5f;
359  x0 = 0.;
360  return 0;
361  }
362  if (fx * fa > 0.f) {
363  xa = x0;
364  fa = fx;
365  } else {
366  xb = x0;
367  fb = fx;
368  }
369  L3:
370  u1 = f1 - f2;
371  u2 = x1 - x2;
372  u3 = f2 - fx;
373  u4 = x2 - x0;
374  if (u2 == 0.f || u4 == 0.f) {goto L1;}
375  f3 = fx;
376  x3 = x0;
377  u1 /= u2;
378  u2 = u3 / u4;
379  ca = u1 - u2;
380  cb = (x1 + x2) * u2 - (x2 + x0) * u1;
381  cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
382  if (ca == 0.f) {
383  if (cb == 0.f) {goto L1;}
384  x0 = -cc / cb;
385  } else {
386  u3 = cb / (ca * 2.f);
387  u4 = u3 * u3 - cc / ca;
388  if (u4 < 0.f) {goto L1;}
389  su4 = std::abs(u4);
390  if (x0 + u3 < 0.f) {su4 = -su4;}
391  x0 = -u3 + su4;
392  }
393  if (x0 < xa || x0 > xb) {goto L1;}
394  // Computing MIN
395  d__3 = (d__1 = x0 - x3, std::abs(d__1));
396  d__4 = (d__2 = x0 - x2, std::abs(d__2));
397  rv = std::min(d__3,d__4);
398  ee = eps * (std::abs(x0) + 1);
399  if (rv > ee) {
400  f1 = f2;
401  x1 = x2;
402  f2 = f3;
403  x2 = x3;
404  goto L2;
405  }
406  fx = func(x0);
407  if (fx == 0.f) {
408  rv = ee;
409  ff = func(x0);
410  return 0;
411  }
412  if (fx * fa < 0.f) {
413  xx = x0 - ee;
414  if (xx <= xa) {
415  rv = ee;
416  ff = func(x0);
417  return 0;
418  }
419  ff = func(xx);
420  fb = ff;
421  xb = xx;
422  } else {
423  xx = x0 + ee;
424  if (xx >= xb) {
425  rv = ee;
426  ff = func(x0);
427  return 0;
428  }
429  ff = func(xx);
430  fa = ff;
431  xa = xx;
432  }
433  if (fx * ff > 0.f) {
434  mc += 2;
435  if (mc > mxf) {
436  rv = (d__1 = xb - xa, std::abs(d__1)) * -0.5f;
437  x0 = 0.f;
438  return 0;
439  }
440  f1 = f3;
441  x1 = x3;
442  f2 = fx;
443  x2 = x0;
444  x0 = xx;
445  fx = ff;
446  goto L3;
447  }
448  /* L4: */
449  rv = ee;
450  ff = func(x0);
451  return 0;
452  } // dzero
453 
454 }
VVIObjF::omega_
float omega_
Definition: VVIObjF.h:41
benchmark_cfg.fb
string fb
Definition: benchmark_cfg.py:14
VVIObjF::t1_
float t1_
Definition: VVIObjF.h:39
abs
#define abs(x)
Definition: mlp_lapack.h:159
min
#define min(a, b)
Definition: mlp_lapack.h:161
alignBH_cfg.mode
list mode
Definition: alignBH_cfg.py:7
q2
double q2[4]
Definition: TauolaWrapper.h:88
VVIObjF::mode_
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
Definition: VVIObjF.h:37
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
VVIObjFDetails::sincosint
void sincosint(float x, float &sint, float &cint)
Definition: VVIObjF.cc:193
p4
double p4[4]
Definition: TauolaWrapper.h:92
VVIObjF::x0_
float x0_
Definition: VVIObjF.h:42
f
double f[11][100]
Definition: MuScleFitUtils.cc:79
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::VVIObjF
VVIObjF(float kappa=0.01, float beta2=1., int mode=0)
Constructor.
Definition: VVIObjF.cc:45
p2
double p2[4]
Definition: TauolaWrapper.h:90
VVIObjF::t0_
float t0_
Definition: VVIObjF.h:38
VVIObjF::fcn
float fcn(float x) const
Definition: VVIObjF.cc:133
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
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
q1
double q1[4]
Definition: TauolaWrapper.h:87
reco::dp
auto dp
Definition: deltaR.h:24
b
double b
Definition: hdecay.h:120
benchmark_cfg.fa
string fa
Definition: benchmark_cfg.py:13
p1
double p1[4]
Definition: TauolaWrapper.h:89
a
double a
Definition: hdecay.h:121
x
Definition: DDAxes.h:10
F
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:281
VVIObjF::limits
void limits(float &xl, float &xu) const
density (mode=0) or distribution (mode=1) function
Definition: VVIObjF.cc:183
VVIObjF::t_
float t_
Definition: VVIObjF.h:40
VVIObjF::a_
float a_[155]
Definition: VVIObjF.h:43
p3
double p3[4]
Definition: TauolaWrapper.h:91
sicif
int sicif(float xx, float &si, float &ci)
Definition: sicif.h:227
Generated for CMSSW Reference Manual by doxygen 1.8.5