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