CMS 3D CMS Logo
VVIObj.cc
Go to the documentation of this file.
1 //
2 // VVIObj.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/SiStripRecHitConverter/interface/VVIObj.h"
17 #else
18 #include "VVIObj.h"
19 #endif
20 
21 #include <cmath>
22 #include <algorithm>
23 #include <functional>
24 
25 namespace sistripvvi {
26 
27  namespace VVIObjDetails {
28  void sincosint(double x, double& sint, double& cint);
29  double cosint(double x);
30  double sinint(double x);
31  double expint(double x);
32 
33  inline double f1(double x, double const* h_) { return h_[0] + h_[1] * std::log(h_[2] * x) - h_[3] * x; }
34  inline double f2(double x, double const* h_) {
35  return h_[4] - x + h_[5] * (std::log(std::abs(x)) + expint(x)) - h_[6] * std::exp(-x);
36  }
37  template <typename F>
38  int dzero(double a, double b, double& x0, double& rv, double eps, int mxf, F func);
39  } // namespace VVIObjDetails
40 
41  // ***************************************************************************************************************************************
47  // ***************************************************************************************************************************************
48 
49  VVIObj::VVIObj(double kappa, double beta2, int mode) : mode_(mode) {
50  const double xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
51  const double xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
52  double h_[7];
53  double q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
54  int lp, lq, k, l, n;
55 
56  // Make sure that the inputs are reasonable
57 
58  if (kappa < 0.01)
59  kappa = 0.01;
60  if (kappa > 10.)
61  kappa = 10.;
62  if (beta2 < 0.)
63  beta2 = 0.;
64  if (beta2 > 1.)
65  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])
79  break;
80  }
81  ll = -lp - 1.5;
82  for (lq = 0; lq < 7; ++lq) {
83  if (kappa <= xq[lq])
84  break;
85  }
86  ul = lq - 6.5;
87  // double (*fp2)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f2);
88  VVIObjDetails::dzero(ll, ul, u, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f2, std::placeholders::_1, h_));
89  q = 1. / u;
90  t1_ = h4 * q - h5 - (beta2 * q + 1) * (log((fabs(u))) + VVIObjDetails::expint(u)) + exp(-u) * q;
91  t_ = t1_ - t0_;
92  omega_ = 6.2831853000000004 / t_;
93  h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
94  if (kappa >= .07) {
95  h_[0] += 6.90775527;
96  }
97  h_[1] = beta2 * kappa;
98  h_[2] = h6 * omega_;
99  h_[3] = omega_ * 1.5707963250000001;
100  // double (*fp1)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f1);
101  VVIObjDetails::dzero(5., 155., x0_, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f1, std::placeholders::_1, h_));
102  n = x0_ + 1.;
103  d = exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
104  a_[n - 1] = 0.;
105  if (mode_ == 0) {
106  a_[n - 1] = omega_ * .31830988654751274;
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  VVIObjDetails::sincosint(x1, c2, c1);
115  c1 = log(x) - c1;
116  c3 = sin(x1);
117  c4 = cos(x1);
118  xf1 = kappa * (beta2 * c1 - c4) - x * c2;
119  xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
120  if (mode_ == 0) {
121  d1 = q * d * omega_ * exp(xf1);
122  a_[l - 1] = d1 * cos(xf2);
123  b_[l - 1] = -d1 * sin(xf2);
124  } else {
125  d1 = q * d * exp(xf1) / k;
126  a_[l - 1] = d1 * sin(xf2);
127  b_[l - 1] = d1 * cos(xf2);
128  a_[n - 1] += q2 * a_[l - 1];
129  }
130  q = -q;
131  q2 = -q2;
132  }
133 
134  } // VVIObj
135 
136  // *************************************************************************************************************************************
140  // *************************************************************************************************************************************
141 
142  double VVIObj::fcn(double x) const {
143  // Local variables
144 
145  double f, u, y, a0, a1;
146  double a2 = 0.;
147  double b1, b0, b2, cof;
148  int k, n, n1;
149 
150  n = x0_;
151  if (x < t0_) {
152  f = 0.;
153  } else if (x <= t1_) {
154  y = x - t0_;
155  u = omega_ * y - 3.141592653589793;
156  cof = cos(u) * 2.;
157  a1 = 0.;
158  a0 = a_[0];
159  n1 = n + 1;
160  for (k = 2; k <= n1; ++k) {
161  a2 = a1;
162  a1 = a0;
163  a0 = a_[k - 1] + cof * a1 - a2;
164  }
165  b1 = 0.;
166  b0 = b_[0];
167  for (k = 2; k <= n; ++k) {
168  b2 = b1;
169  b1 = b0;
170  b0 = b_[k - 1] + cof * b1 - b2;
171  }
172  f = (a0 - a2) * .5 + b0 * sin(u);
173  if (mode_ != 0) {
174  f += y / t_;
175  }
176  } else {
177  f = 0.;
178  if (mode_ != 0) {
179  f = 1.;
180  }
181  }
182  return f;
183  } // fcn
184 
185  // *************************************************************************************************************************************
189  // *************************************************************************************************************************************
190 
191  void VVIObj::limits(double& xl, double& xu) const {
192  xl = t0_;
193  xu = t1_;
194  return;
195  } // limits
196 
197  namespace VVIObjDetails {
198  double cosint(double x) {
199  // Initialized data
200 
201  const double zero = 0.;
202  const double one = 1.;
203  const double two = 2.;
204  const double eight = 8.;
205  const double ce = .57721566490153;
206  const double c__[14] = {1.9405491464836,
207  .9413409132865,
208  -.579845034293,
209  .3091572011159,
210  -.0916101792208,
211  .0164437407515,
212  -.0019713091952,
213  1.692538851e-4,
214  -1.09393296e-5,
215  5.522386e-7,
216  -2.23995e-8,
217  7.465e-10,
218  -2.08e-11,
219  5e-13};
220  const double p[23] = {
221  .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
222  6.0069212e-7, -1.2586794e-7, 2.932563e-8, -7.45696e-9, 2.04105e-9, -5.9502e-10,
223  1.8323e-10, -5.921e-11, 1.997e-11, -7e-12, 2.54e-12, -9.5e-13,
224  3.7e-13, -1.4e-13, 6e-14, -2e-14, 1e-14};
225  const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
226  -4.2110196e-7, 6.907245e-8, -1.318321e-8, 2.83697e-9, -6.7329e-10,
227  1.734e-10, -4.787e-11, 1.403e-11, -4.33e-12, 1.4e-12,
228  -4.7e-13, 1.7e-13, -6e-14, 2e-14, -1e-14};
229 
230  // System generated locals
231  double d__1;
232 
233  // Local variables
234  double h__;
235  int i__;
236  double r__, y, b0, b1, b2, pp, qq, alfa;
237 
238  // If x==0, return same
239 
240  if (x == zero) {
241  return zero;
242  }
243  if (fabs(x) <= eight) {
244  y = x / eight;
245  // Computing 2nd power
246  d__1 = y;
247  h__ = two * (d__1 * d__1) - one;
248  alfa = -two * h__;
249  b1 = zero;
250  b2 = zero;
251  for (i__ = 13; i__ >= 0; --i__) {
252  b0 = c__[i__] - alfa * b1 - b2;
253  b2 = b1;
254  b1 = b0;
255  }
256  b1 = ce + log((fabs(x))) - b0 + h__ * b2;
257  } else {
258  r__ = one / x;
259  y = eight * r__;
260  // Computing 2nd power
261  d__1 = y;
262  h__ = two * (d__1 * d__1) - one;
263  alfa = -two * h__;
264  b1 = zero;
265  b2 = zero;
266  for (i__ = 22; i__ >= 0; --i__) {
267  b0 = p[i__] - alfa * b1 - b2;
268  b2 = b1;
269  b1 = b0;
270  }
271  pp = b0 - h__ * b2;
272  b1 = zero;
273  b2 = zero;
274  for (i__ = 19; i__ >= 0; --i__) {
275  b0 = q[i__] - alfa * b1 - b2;
276  b2 = b1;
277  b1 = b0;
278  }
279  qq = b0 - h__ * b2;
280  b1 = r__ * (qq * sin(x) - r__ * pp * cos(x));
281  }
282  return b1;
283  } // cosint
284 
285  double sinint(double x) {
286  // Initialized data
287 
288  const double zero = 0.;
289  const double one = 1.;
290  const double two = 2.;
291  const double eight = 8.;
292  const double pih = 1.5707963267949;
293  const double s[14] = {1.9522209759531,
294  -.6884042321257,
295  .4551855132256,
296  -.1804571236838,
297  .0410422133759,
298  -.0059586169556,
299  6.001427414e-4,
300  -4.44708329e-5,
301  2.5300782e-6,
302  -1.141308e-7,
303  4.1858e-9,
304  -1.273e-10,
305  3.3e-12,
306  -1e-13};
307  const double p[23] = {
308  .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
309  6.0069212e-7, -1.2586794e-7, 2.932563e-8, -7.45696e-9, 2.04105e-9, -5.9502e-10,
310  1.8323e-10, -5.921e-11, 1.997e-11, -7e-12, 2.54e-12, -9.5e-13,
311  3.7e-13, -1.4e-13, 6e-14, -2e-14, 1e-14};
312  const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
313  -4.2110196e-7, 6.907245e-8, -1.318321e-8, 2.83697e-9, -6.7329e-10,
314  1.734e-10, -4.787e-11, 1.403e-11, -4.33e-12, 1.4e-12,
315  -4.7e-13, 1.7e-13, -6e-14, 2e-14, -1e-14};
316 
317  // System generated locals
318  double d__1;
319 
320  // Local variables
321  double h__;
322  int i__;
323  double r__, y, b0, b1, b2, pp, qq, alfa;
324 
325  if (fabs(x) <= eight) {
326  y = x / eight;
327  d__1 = y;
328  h__ = two * (d__1 * d__1) - one;
329  alfa = -two * h__;
330  b1 = zero;
331  b2 = zero;
332  for (i__ = 13; i__ >= 0; --i__) {
333  b0 = s[i__] - alfa * b1 - b2;
334  b2 = b1;
335  b1 = b0;
336  }
337  b1 = y * (b0 - b2);
338  } else {
339  r__ = one / x;
340  y = eight * r__;
341  d__1 = y;
342  h__ = two * (d__1 * d__1) - one;
343  alfa = -two * h__;
344  b1 = zero;
345  b2 = zero;
346  for (i__ = 22; i__ >= 0; --i__) {
347  b0 = p[i__] - alfa * b1 - b2;
348  b2 = b1;
349  b1 = b0;
350  }
351  pp = b0 - h__ * b2;
352  b1 = zero;
353  b2 = zero;
354  for (i__ = 19; i__ >= 0; --i__) {
355  b0 = q[i__] - alfa * b1 - b2;
356  b2 = b1;
357  b1 = b0;
358  }
359  qq = b0 - h__ * b2;
360  d__1 = fabs(pih);
361  if (x < 0.)
362  d__1 = -d__1;
363  b1 = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
364  }
365 
366  return b1;
367  } // sinint
368 
369  void sincosint(double x, double& sint, double& cint) {
370  // Initialized data
371 
372  const double zero = 0.;
373  const double one = 1.;
374  const double two = 2.;
375  const double eight = 8.;
376  const double ce = .57721566490153;
377  const double pih = 1.5707963267949;
378  const double s__[14] = {1.9522209759531,
379  -.6884042321257,
380  .4551855132256,
381  -.1804571236838,
382  .0410422133759,
383  -.0059586169556,
384  6.001427414e-4,
385  -4.44708329e-5,
386  2.5300782e-6,
387  -1.141308e-7,
388  4.1858e-9,
389  -1.273e-10,
390  3.3e-12,
391  -1e-13};
392 
393  const double c__[14] = {1.9405491464836,
394  .9413409132865,
395  -.579845034293,
396  .3091572011159,
397  -.0916101792208,
398  .0164437407515,
399  -.0019713091952,
400  1.692538851e-4,
401  -1.09393296e-5,
402  5.522386e-7,
403  -2.23995e-8,
404  7.465e-10,
405  -2.08e-11,
406  5e-13};
407 
408  const double p[23] = {
409  .96074783975204, -.0371138962124, .00194143988899, -1.7165988425e-4, 2.112637753e-5, -3.27163257e-6,
410  6.0069212e-7, -1.2586794e-7, 2.932563e-8, -7.45696e-9, 2.04105e-9, -5.9502e-10,
411  1.8323e-10, -5.921e-11, 1.997e-11, -7e-12, 2.54e-12, -9.5e-13,
412  3.7e-13, -1.4e-13, 6e-14, -2e-14, 1e-14};
413  const double q[20] = {.98604065696238, -.0134717382083, 4.5329284117e-4, -3.067288652e-5, 3.13199198e-6,
414  -4.2110196e-7, 6.907245e-8, -1.318321e-8, 2.83697e-9, -6.7329e-10,
415  1.734e-10, -4.787e-11, 1.403e-11, -4.33e-12, 1.4e-12,
416  -4.7e-13, 1.7e-13, -6e-14, 2e-14, -1e-14};
417 
418  // System generated locals
419  double d__1;
420 
421  // Local variables
422  double h__;
423  int i__;
424  double r__, y, b0, b1, b2, pp, qq, alfa;
425 
426  sint = 0;
427  cint = 0;
428 
429  if (fabs(x) <= eight) {
430  y = x / eight;
431  // Computing 2nd power
432  d__1 = y;
433  h__ = two * (d__1 * d__1) - one;
434  alfa = -two * h__;
435 
436  // cos
437  if (x != 0) {
438  b1 = zero;
439  b2 = zero;
440  for (i__ = 13; i__ >= 0; --i__) {
441  b0 = c__[i__] - alfa * b1 - b2;
442  b2 = b1;
443  b1 = b0;
444  }
445  cint = ce + log((fabs(x))) - b0 + h__ * b2;
446  }
447  // sin
448  b1 = zero;
449  b2 = zero;
450  for (i__ = 13; i__ >= 0; --i__) {
451  b0 = s__[i__] - alfa * b1 - b2;
452  b2 = b1;
453  b1 = b0;
454  }
455  sint = y * (b0 - b2);
456 
457  } else {
458  r__ = one / x;
459  y = eight * r__;
460  // Computing 2nd power
461  d__1 = y;
462  h__ = two * (d__1 * d__1) - one;
463  alfa = -two * h__;
464  b1 = zero;
465  b2 = zero;
466  for (i__ = 22; i__ >= 0; --i__) {
467  b0 = p[i__] - alfa * b1 - b2;
468  b2 = b1;
469  b1 = b0;
470  }
471  pp = b0 - h__ * b2;
472  b1 = zero;
473  b2 = zero;
474  for (i__ = 19; i__ >= 0; --i__) {
475  b0 = q[i__] - alfa * b1 - b2;
476  b2 = b1;
477  b1 = b0;
478  }
479  qq = b0 - h__ * b2;
480  // cos
481  cint = r__ * (qq * sin(x) - r__ * pp * cos(x));
482  // sin
483  d__1 = pih;
484  if (x < 0.)
485  d__1 = -d__1;
486  sint = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
487  }
488  }
489 
490  double expint(double x) {
491  // Initialized data
492 
493  const double zero = 0.;
494  const double q2[7] = {
495  .10340013040487, 3.319092135933, 20.449478501379, 41.280784189142, 32.426421069514, 10.041164382905, 1.};
496  const double p3[6] = {
497  -2.3909964453136, -147.98219500504, -254.3763397689, -119.55761038372, -19.630408535939, -.9999999999036};
498  const double q3[6] = {177.60070940351, 530.68509610812, 462.23027156148, 156.81843364539, 21.630408494238, 1.};
499  const double p4[8] = {-8.6693733995107,
500  -549.14226552109,
501  -4210.0161535707,
502  -249301.39345865,
503  -119623.66934925,
504  -22174462.775885,
505  3892804.213112,
506  -391546073.8091};
507  const double q4[8] = {34.171875,
508  -1607.0892658722,
509  35730.029805851,
510  -483547.43616216,
511  4285596.2461175,
512  -24903337.574054,
513  89192576.757561,
514  -165254299.72521};
515  const double a1[8] = {-2.1808638152072,
516  -21.901023385488,
517  9.3081638566217,
518  25.076281129356,
519  -33.184253199722,
520  60.121799083008,
521  -43.253113287813,
522  1.0044310922808};
523  const double b1[8] = {0.,
524  3.9370770185272,
525  300.89264837292,
526  -6.2504116167188,
527  1003.6743951673,
528  14.325673812194,
529  2736.2411988933,
530  .52746885196291};
531  const double a2[8] = {-3.4833465360285,
532  -18.65454548834,
533  -8.2856199414064,
534  -32.34673303054,
535  17.960168876925,
536  1.7565631546961,
537  -1.9502232128966,
538  .99999429607471};
539  const double b2[8] = {0.,
540  69.500065588743,
541  57.283719383732,
542  25.777638423844,
543  760.76114800773,
544  28.951672792514,
545  -3.4394226689987,
546  1.0008386740264};
547  const double a3[6] = {
548  -27.780928934438, -10.10479081576, -9.1483008216736, -5.0223317461851, -3.0000077799358, 1.0000000000704};
549  const double one = 1.;
550  const double b3[6] = {0., 122.39993926823, 2.7276100778779, -7.1897518395045, -2.9990118065262, 1.999999942826};
551  const double two = 2.;
552  const double three = 3.;
553  const double x0 = .37250741078137;
554  const double xl[6] = {-24., -12., -6., 0., 1., 4.};
555  const double p1[5] = {4.293125234321, 39.894153870321, 292.52518866921, 425.69682638592, -434.98143832952};
556  const double q1[5] = {1., 18.899288395003, 150.95038744251, 568.05252718987, 753.58564359843};
557  const double p2[7] = {.43096783946939,
558  6.9052252278444,
559  23.019255939133,
560  24.378408879132,
561  9.0416155694633,
562  .99997957705159,
563  4.656271079751e-7};
564 
565  /* Local variables */
566  double v, y, ap, bp, aq, dp, bq, dq;
567 
568  if (x <= xl[0]) {
569  ap = a3[0] - x;
570  for (int i__ = 2; i__ <= 5; ++i__) {
571  /* L1: */
572  ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
573  }
574  y = exp(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
575  } else if (x <= xl[1]) {
576  ap = a2[0] - x;
577  for (int i__ = 2; i__ <= 7; ++i__) {
578  ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
579  }
580  y = exp(-x) / x * (a2[7] + b2[7] / ap);
581  } else if (x <= xl[2]) {
582  ap = a1[0] - x;
583  for (int i__ = 2; i__ <= 7; ++i__) {
584  ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
585  }
586  y = exp(-x) / x * (a1[7] + b1[7] / ap);
587  } else if (x < xl[3]) {
588  v = -two * (x / three + one);
589  bp = zero;
590  dp = p4[0];
591  for (int i__ = 2; i__ <= 8; ++i__) {
592  ap = bp;
593  bp = dp;
594  dp = p4[i__ - 1] - ap + v * bp;
595  }
596  bq = zero;
597  dq = q4[0];
598  for (int i__ = 2; i__ <= 8; ++i__) {
599  aq = bq;
600  bq = dq;
601  dq = q4[i__ - 1] - aq + v * bq;
602  }
603  y = -log(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
604  } else if (x == xl[3]) {
605  return zero;
606  } else if (x < xl[4]) {
607  ap = p1[0];
608  aq = q1[0];
609  for (int i__ = 2; i__ <= 5; ++i__) {
610  ap = p1[i__ - 1] + x * ap;
611  aq = q1[i__ - 1] + x * aq;
612  }
613  y = -log(x) + ap / aq;
614  } else if (x <= xl[5]) {
615  y = one / x;
616  ap = p2[0];
617  aq = q2[0];
618  for (int i__ = 2; i__ <= 7; ++i__) {
619  ap = p2[i__ - 1] + y * ap;
620  aq = q2[i__ - 1] + y * aq;
621  }
622  y = exp(-x) * ap / aq;
623  } else {
624  y = one / x;
625  ap = p3[0];
626  aq = q3[0];
627  for (int i__ = 2; i__ <= 6; ++i__) {
628  ap = p3[i__ - 1] + y * ap;
629  aq = q3[i__ - 1] + y * aq;
630  }
631  y = exp(-x) * y * (one + y * ap / aq);
632  }
633  return y;
634  } // expint
635 
636  template <typename F>
637  int dzero(double a, double b, double& x0, double& rv, double eps, int mxf, F func) {
638  /* System generated locals */
639  double d__1, d__2, d__3, d__4;
640 
641  // Local variables
642  double f1, f2, f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc, fa, fb, ee, ff;
643  int mc;
644  double xa, xb, fx, xx, su4;
645 
646  xa = std::min(a, b);
647  xb = std::max(a, b);
648  fa = func(xa);
649  fb = func(xb);
650  if (fa * fb > 0.) {
651  rv = (xb - xa) * -2;
652  x0 = 0.;
653  return 1;
654  }
655  mc = 0;
656  L1:
657  x0 = (xa + xb) * .5;
658  rv = x0 - xa;
659  ee = eps * (fabs(x0) + 1);
660  if (rv <= ee) {
661  rv = ee;
662  ff = func(x0);
663  return 0;
664  }
665  f1 = fa;
666  x1 = xa;
667  f2 = fb;
668  x2 = xb;
669  L2:
670  fx = func(x0);
671  ++mc;
672  if (mc > mxf) {
673  rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
674  x0 = 0.;
675  return 0;
676  }
677  if (fx * fa > 0.) {
678  xa = x0;
679  fa = fx;
680  } else {
681  xb = x0;
682  fb = fx;
683  }
684  L3:
685  u1 = f1 - f2;
686  u2 = x1 - x2;
687  u3 = f2 - fx;
688  u4 = x2 - x0;
689  if (u2 == 0. || u4 == 0.) {
690  goto L1;
691  }
692  f3 = fx;
693  x3 = x0;
694  u1 /= u2;
695  u2 = u3 / u4;
696  ca = u1 - u2;
697  cb = (x1 + x2) * u2 - (x2 + x0) * u1;
698  cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
699  if (ca == 0.) {
700  if (cb == 0.) {
701  goto L1;
702  }
703  x0 = -cc / cb;
704  } else {
705  u3 = cb / (ca * 2);
706  u4 = u3 * u3 - cc / ca;
707  if (u4 < 0.) {
708  goto L1;
709  }
710  su4 = fabs(u4);
711  if (x0 + u3 < 0.f) {
712  su4 = -su4;
713  }
714  x0 = -u3 + su4;
715  }
716  if (x0 < xa || x0 > xb) {
717  goto L1;
718  }
719  // Computing MIN
720  d__3 = (d__1 = x0 - x3, fabs(d__1)), d__4 = (d__2 = x0 - x2, fabs(d__2));
721  rv = std::min(d__3, d__4);
722  ee = eps * (fabs(x0) + 1);
723  if (rv > ee) {
724  f1 = f2;
725  x1 = x2;
726  f2 = f3;
727  x2 = x3;
728  goto L2;
729  }
730  fx = func(x0);
731  if (fx == 0.) {
732  rv = ee;
733  ff = func(x0);
734  return 0;
735  }
736  if (fx * fa < 0.) {
737  xx = x0 - ee;
738  if (xx <= xa) {
739  rv = ee;
740  ff = func(x0);
741  return 0;
742  }
743  ff = func(xx);
744  fb = ff;
745  xb = xx;
746  } else {
747  xx = x0 + ee;
748  if (xx >= xb) {
749  rv = ee;
750  ff = func(x0);
751  return 0;
752  }
753  ff = func(xx);
754  fa = ff;
755  xa = xx;
756  }
757  if (fx * ff > 0.) {
758  mc += 2;
759  if (mc > mxf) {
760  rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
761  x0 = 0.;
762  return 0;
763  }
764  f1 = f3;
765  x1 = x3;
766  f2 = fx;
767  x2 = x0;
768  x0 = xx;
769  fx = ff;
770  goto L3;
771  }
772  /* L4: */
773  rv = ee;
774  ff = func(x0);
775  return 0;
776  } // dzero
777 
778  } // namespace VVIObjDetails
779 } // namespace sistripvvi
EcalMonitorTask_cff.func
func
Definition: EcalMonitorTask_cff.py:10
CaloTowersParam_cfi.mc
mc
Definition: CaloTowersParam_cfi.py:8
b3
weight_default_t b3[10]
Definition: b3.h:9
b1
weight_default_t b1[25]
Definition: b1.h:9
sistripvvi::VVIObjDetails::expint
double expint(double x)
Private version of the sine integral.
Definition: VVIObj.cc:490
sistripvvi::VVIObj::b_
double b_[155]
Definition: VVIObj.h:43
sistripvvi::VVIObjDetails::sinint
double sinint(double x)
Private version of the cosine integral.
Definition: VVIObj.cc:285
SiStripOfflineCRack_cfg.p2
p2
Definition: SiStripOfflineCRack_cfg.py:140
multPhiCorr_741_25nsDY_cfi.fx
fx
Definition: multPhiCorr_741_25nsDY_cfi.py:9
sistripvvi::VVIObjDetails::sincosint
void sincosint(double x, double &sint, double &cint)
Definition: VVIObj.cc:369
Calorimetry_cff.dp
dp
Definition: Calorimetry_cff.py:158
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:36
testProducerWithPsetDescEmpty_cfi.u3
u3
Definition: testProducerWithPsetDescEmpty_cfi.py:50
testProducerWithPsetDescEmpty_cfi.u1
u1
Definition: testProducerWithPsetDescEmpty_cfi.py:49
sistripvvi::VVIObjDetails::dzero
int dzero(double a, double b, double &x0, double &rv, double eps, int mxf, F func)
Definition: VVIObj.cc:637
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
LaserDQM_cfg.p1
p1
Definition: LaserDQM_cfg.py:42
ALCARECOPromptCalibProdSiPixelAli0T_cff.mode
mode
Definition: ALCARECOPromptCalibProdSiPixelAli0T_cff.py:96
sistripvvi::VVIObj::limits
void limits(double &xl, double &xu) const
density (mode=0) or distribution (mode=1) function
Definition: VVIObj.cc:191
sistripvvi::VVIObjDetails::cosint
double cosint(double x)
Private version of the cosine and sine integral.
Definition: VVIObj.cc:198
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:33
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
testProducerWithPsetDescEmpty_cfi.a2
a2
Definition: testProducerWithPsetDescEmpty_cfi.py:35
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
sistripvvi::VVIObj::omega_
double omega_
Definition: VVIObj.h:40
b2
weight_default_t b2[10]
Definition: b2.h:9
ztail.d
d
Definition: ztail.py:151
sistripvvi::VVIObjDetails::f2
double f2(double x, double const *h_)
Definition: VVIObj.cc:34
a0
static constexpr float a0
Definition: L1EGammaCrystalsEmulatorProducer.cc:79
chargedHadronTrackResolutionFilter_cfi.p3
p3
Definition: chargedHadronTrackResolutionFilter_cfi.py:12
b
double b
Definition: hdecay.h:118
a
double a
Definition: hdecay.h:119
HLT_2022v15_cff.eps
eps
Definition: HLT_2022v15_cff.py:8703
hgcalLayerClusters_cfi.kappa
kappa
Definition: hgcalLayerClusters_cfi.py:35
x
float x
Definition: beamSpotDipStandalone.cc:55
sistripvvi::VVIObj::fcn
double fcn(double x) const
Definition: VVIObj.cc:142
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:80
F
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:163
sistripvvi::VVIObj::VVIObj
VVIObj(double kappa=0.01, double beta2=1., int mode=0)
Constructor.
Definition: VVIObj.cc:49
d1
static constexpr float d1
Definition: L1EGammaCrystalsEmulatorProducer.cc:82
sistripvvi::VVIObj::a_
double a_[155]
Definition: VVIObj.h:42
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
Generated for CMSSW Reference Manual by doxygen 1.8.14