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