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