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