CMS 3D CMS Logo

Functions
VVIObjDetails Namespace Reference

Functions

double cosint (double x)
 Private version of the cosine and sine integral. More...
 
template<typename F >
int dzero (double a, double b, double &x0, double &rv, double eps, int mxf, F func)
 Private version of the exponential integral. More...
 
double expint (double x)
 Private version of the sine integral. More...
 
void sincosint (double x, double &sint, double &cint)
 
double sinint (double x)
 Private version of the cosine integral. More...
 

Function Documentation

double VVIObjDetails::cosint ( double  x)

Private version of the cosine and sine integral.

Definition at line 193 of file VVIObj.cc.

References funct::cos(), MillePedeFileConverter_cfg::e, cmsBatch::log, AlCaHLTBitMon_ParallelJobs::p, createTree::pp, lumiQueryAPI::q, funct::sin(), x, and y.

193  {
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,.9413409132865,
202  -.579845034293,.3091572011159,-.0916101792208,.0164437407515,
203  -.0019713091952,1.692538851e-4,-1.09393296e-5,5.522386e-7,
204  -2.23995e-8,7.465e-10,-2.08e-11,5e-13 };
205  const double p[23] = { .96074783975204,-.0371138962124,
206  .00194143988899,-1.7165988425e-4,2.112637753e-5,-3.27163257e-6,
207  6.0069212e-7,-1.2586794e-7,2.932563e-8,-7.45696e-9,2.04105e-9,
208  -5.9502e-10,1.8323e-10,-5.921e-11,1.997e-11,-7e-12,2.54e-12,
209  -9.5e-13,3.7e-13,-1.4e-13,6e-14,-2e-14,1e-14 };
210  const double q[20] = { .98604065696238,-.0134717382083,
211  4.5329284117e-4,-3.067288652e-5,3.13199198e-6,-4.2110196e-7,
212  6.907245e-8,-1.318321e-8,2.83697e-9,-6.7329e-10,1.734e-10,
213  -4.787e-11,1.403e-11,-4.33e-12,1.4e-12,-4.7e-13,1.7e-13,-6e-14,
214  2e-14,-1e-14 };
215 
216  // System generated locals
217  double d__1;
218 
219  // Local variables
220  double h__;
221  int i__;
222  double r__, y, b0, b1, b2, pp, qq, alfa;
223 
224  // If x==0, return same
225 
226  if (x == zero) {
227  return zero;
228  }
229  if (fabs(x) <= eight) {
230  y = x / eight;
231  // Computing 2nd power
232  d__1 = y;
233  h__ = two * (d__1 * d__1) - one;
234  alfa = -two * h__;
235  b1 = zero;
236  b2 = zero;
237  for (i__ = 13; i__ >= 0; --i__) {
238  b0 = c__[i__] - alfa * b1 - b2;
239  b2 = b1;
240  b1 = b0;
241  }
242  b1 = ce + log((fabs(x))) - b0 + h__ * b2;
243  } else {
244  r__ = one / x;
245  y = eight * r__;
246  // Computing 2nd power
247  d__1 = y;
248  h__ = two * (d__1 * d__1) - one;
249  alfa = -two * h__;
250  b1 = zero;
251  b2 = zero;
252  for (i__ = 22; i__ >= 0; --i__) {
253  b0 = p[i__] - alfa * b1 - b2;
254  b2 = b1;
255  b1 = b0;
256  }
257  pp = b0 - h__ * b2;
258  b1 = zero;
259  b2 = zero;
260  for (i__ = 19; i__ >= 0; --i__) {
261  b0 = q[i__] - alfa * b1 - b2;
262  b2 = b1;
263  b1 = b0;
264  }
265  qq = b0 - h__ * b2;
266  b1 = r__ * (qq * sin(x) - r__ * pp * cos(x));
267  }
268  return b1;
269  } // cosint
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
template<typename F >
int VVIObjDetails::dzero ( double  a,
double  b,
double &  x0,
double &  rv,
double  eps,
int  mxf,
F  func 
)

Private version of the exponential integral.

Definition at line 570 of file VVIObj.cc.

References f, sistripvvi::VVIObjDetails::f1(), sistripvvi::VVIObjDetails::f2(), connectstrParser::f3, benchmark_cfg::fa, benchmark_cfg::fb, alignCSCRings::ff, jets_cff::func, SiStripPI::max, CaloTowersParam_cfi::mc, min(), MetAnalyzer::u1, MetAnalyzer::u2, globals_cff::x1, globals_cff::x2, and geometryCSVtoXML::xx.

Referenced by VVIObj::VVIObj().

571  {
572  /* System generated locals */
573  double d__1, d__2, d__3, d__4;
574 
575  // Local variables
576  double f1, f2, f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc, fa, fb, ee, ff;
577  int mc;
578  double xa, xb, fx, xx, su4;
579 
580  xa = std::min(a,b);
581  xb = std::max(a,b);
582  fa = func(xa);
583  fb = func(xb);
584  if (fa * fb > 0.) {
585  rv = (xb - xa) * -2;
586  x0 = 0.;
587  return 1;
588  }
589  mc = 0;
590  L1:
591  x0 = (xa + xb) * .5;
592  rv = x0 - xa;
593  ee = eps * (fabs(x0) + 1);
594  if (rv <= ee) {
595  rv = ee;
596  ff = func(x0);
597  return 0;
598  }
599  f1 = fa;
600  x1 = xa;
601  f2 = fb;
602  x2 = xb;
603  L2:
604  fx = func(x0);
605  ++mc;
606  if (mc > mxf) {
607  rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
608  x0 = 0.;
609  return 0;
610  }
611  if (fx * fa > 0.) {
612  xa = x0;
613  fa = fx;
614  } else {
615  xb = x0;
616  fb = fx;
617  }
618  L3:
619  u1 = f1 - f2;
620  u2 = x1 - x2;
621  u3 = f2 - fx;
622  u4 = x2 - x0;
623  if (u2 == 0. || u4 == 0.) {goto L1;}
624  f3 = fx;
625  x3 = x0;
626  u1 /= u2;
627  u2 = u3 / u4;
628  ca = u1 - u2;
629  cb = (x1 + x2) * u2 - (x2 + x0) * u1;
630  cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
631  if (ca == 0.) {
632  if (cb == 0.) {goto L1;}
633  x0 = -cc / cb;
634  } else {
635  u3 = cb / (ca * 2);
636  u4 = u3 * u3 - cc / ca;
637  if (u4 < 0.) {goto L1;}
638  su4 = fabs(u4);
639  if (x0 + u3 < 0.f) {su4 = -su4;}
640  x0 = -u3 + su4;
641  }
642  if (x0 < xa || x0 > xb) {goto L1;}
643  // Computing MIN
644  d__3 = (d__1 = x0 - x3, fabs(d__1)), d__4 = (d__2 = x0 - x2, fabs(d__2));
645  rv = std::min(d__3,d__4);
646  ee = eps * (fabs(x0) + 1);
647  if (rv > ee) {
648  f1 = f2;
649  x1 = x2;
650  f2 = f3;
651  x2 = x3;
652  goto L2;
653  }
654  fx = func(x0);
655  if (fx == 0.) {
656  rv = ee;
657  ff = func(x0);
658  return 0;
659  }
660  if (fx * fa < 0.) {
661  xx = x0 - ee;
662  if (xx <= xa) {
663  rv = ee;
664  ff = func(x0);
665  return 0;
666  }
667  ff = func(xx);
668  fb = ff;
669  xb = xx;
670  } else {
671  xx = x0 + ee;
672  if (xx >= xb) {
673  rv = ee;
674  ff = func(x0);
675  return 0;
676  }
677  ff = func(xx);
678  fa = ff;
679  xa = xx;
680  }
681  if (fx * ff > 0.) {
682  mc += 2;
683  if (mc > mxf) {
684  rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
685  x0 = 0.;
686  return 0;
687  }
688  f1 = f3;
689  x1 = x3;
690  f2 = fx;
691  x2 = x0;
692  x0 = xx;
693  fx = ff;
694  goto L3;
695  }
696  /* L4: */
697  rv = ee;
698  ff = func(x0);
699  return 0;
700  } // dzero
double f1(double x, double const *h_)
Private version of the exponential integral.
Definition: VVIObj.cc:35
double f[11][100]
T min(T a, T b)
Definition: MathUtil.h:58
double f2(double x, double const *h_)
Definition: VVIObj.cc:36
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
double VVIObjDetails::expint ( double  x)

Private version of the sine integral.

Definition at line 448 of file VVIObj.cc.

References JetChargeProducer_cfi::exp, cmsBatch::log, p1, p2, p3, p4, q1, q2, findQualityFiles::v, x, and y.

Referenced by VVIObj::VVIObj().

448  {
449 
450  // Initialized data
451 
452  const double zero = 0.;
453  const double q2[7] = { .10340013040487,3.319092135933,
454  20.449478501379,41.280784189142,32.426421069514,10.041164382905,
455  1. };
456  const double p3[6] = { -2.3909964453136,-147.98219500504,
457  -254.3763397689,-119.55761038372,-19.630408535939,-.9999999999036
458  };
459  const double q3[6] = { 177.60070940351,530.68509610812,
460  462.23027156148,156.81843364539,21.630408494238,1. };
461  const double p4[8] = { -8.6693733995107,-549.14226552109,
462  -4210.0161535707,-249301.39345865,-119623.66934925,
463  -22174462.775885,3892804.213112,-391546073.8091 };
464  const double q4[8] = { 34.171875,-1607.0892658722,35730.029805851,
465  -483547.43616216,4285596.2461175,-24903337.574054,89192576.757561,
466  -165254299.72521 };
467  const double a1[8] = { -2.1808638152072,-21.901023385488,
468  9.3081638566217,25.076281129356,-33.184253199722,60.121799083008,
469  -43.253113287813,1.0044310922808 };
470  const double b1[8] = { 0.,3.9370770185272,300.89264837292,
471  -6.2504116167188,1003.6743951673,14.325673812194,2736.2411988933,
472  .52746885196291 };
473  const double a2[8] = { -3.4833465360285,-18.65454548834,
474  -8.2856199414064,-32.34673303054,17.960168876925,1.7565631546961,
475  -1.9502232128966,.99999429607471 };
476  const double b2[8] = { 0.,69.500065588743,57.283719383732,
477  25.777638423844,760.76114800773,28.951672792514,-3.4394226689987,
478  1.0008386740264 };
479  const double a3[6] = { -27.780928934438,-10.10479081576,
480  -9.1483008216736,-5.0223317461851,-3.0000077799358,
481  1.0000000000704 };
482  const double one = 1.;
483  const double b3[6] = { 0.,122.39993926823,2.7276100778779,
484  -7.1897518395045,-2.9990118065262,1.999999942826 };
485  const double two = 2.;
486  const double three = 3.;
487  const double x0 = .37250741078137;
488  const double xl[6] = { -24.,-12.,-6.,0.,1.,4. };
489  const double p1[5] = { 4.293125234321,39.894153870321,
490  292.52518866921,425.69682638592,-434.98143832952 };
491  const double q1[5] = { 1.,18.899288395003,150.95038744251,
492  568.05252718987,753.58564359843 };
493  const double p2[7] = { .43096783946939,6.9052252278444,
494  23.019255939133,24.378408879132,9.0416155694633,.99997957705159,
495  4.656271079751e-7 };
496 
497  /* Local variables */
498  double v, y, ap, bp, aq, dp, bq, dq;
499 
500  if (x <= xl[0]) {
501  ap = a3[0] - x;
502  for ( int i__ = 2; i__ <= 5; ++i__) {
503  /* L1: */
504  ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
505  }
506  y = exp(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
507  } else if (x <= xl[1]) {
508  ap = a2[0] - x;
509  for ( int i__ = 2; i__ <= 7; ++i__) {
510  ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
511  }
512  y = exp(-x) / x * (a2[7] + b2[7] / ap);
513  } else if (x <= xl[2]) {
514  ap = a1[0] - x;
515  for ( int i__ = 2; i__ <= 7; ++i__) {
516  ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
517  }
518  y = exp(-x) / x * (a1[7] + b1[7] / ap);
519  } else if (x < xl[3]) {
520  v = -two * (x / three + one);
521  bp = zero;
522  dp = p4[0];
523  for ( int i__ = 2; i__ <= 8; ++i__) {
524  ap = bp;
525  bp = dp;
526  dp = p4[i__ - 1] - ap + v * bp;
527  }
528  bq = zero;
529  dq = q4[0];
530  for ( int i__ = 2; i__ <= 8; ++i__) {
531  aq = bq;
532  bq = dq;
533  dq = q4[i__ - 1] - aq + v * bq;
534  }
535  y = -log(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
536  } else if (x == xl[3]) {
537  return zero;
538  } else if (x < xl[4]) {
539  ap = p1[0];
540  aq = q1[0];
541  for ( int i__ = 2; i__ <= 5; ++i__) {
542  ap = p1[i__ - 1] + x * ap;
543  aq = q1[i__ - 1] + x * aq;
544  }
545  y = -log(x) + ap / aq;
546  } else if (x <= xl[5]) {
547  y = one / x;
548  ap = p2[0];
549  aq = q2[0];
550  for ( int i__ = 2; i__ <= 7; ++i__) {
551  ap = p2[i__ - 1] + y * ap;
552  aq = q2[i__ - 1] + y * aq;
553  }
554  y = exp(-x) * ap / aq;
555  } else {
556  y = one / x;
557  ap = p3[0];
558  aq = q3[0];
559  for ( int i__ = 2; i__ <= 6; ++i__) {
560  ap = p3[i__ - 1] + y * ap;
561  aq = q3[i__ - 1] + y * aq;
562  }
563  y = exp(-x) * y * (one + y * ap / aq);
564  }
565  return y;
566 } // expint
double q2[4]
Definition: TauolaWrapper.h:88
double p4[4]
Definition: TauolaWrapper.h:92
double p2[4]
Definition: TauolaWrapper.h:90
double q1[4]
Definition: TauolaWrapper.h:87
double p1[4]
Definition: TauolaWrapper.h:89
double p3[4]
Definition: TauolaWrapper.h:91
void VVIObjDetails::sincosint ( double  x,
double &  sint,
double &  cint 
)

Definition at line 345 of file VVIObj.cc.

References funct::cos(), MillePedeFileConverter_cfg::e, cmsBatch::log, AlCaHLTBitMon_ParallelJobs::p, createTree::pp, lumiQueryAPI::q, funct::sin(), x, and y.

Referenced by VVIObj::VVIObj().

345  {
346  // Initialized data
347 
348  const double zero = 0.;
349  const double one = 1.;
350  const double two = 2.;
351  const double eight = 8.;
352  const double ce = .57721566490153;
353  const double pih = 1.5707963267949;
354  const double s__[14] = { 1.9522209759531,-.6884042321257,
355  .4551855132256,-.1804571236838,.0410422133759,-.0059586169556,
356  6.001427414e-4,-4.44708329e-5,2.5300782e-6,-1.141308e-7,4.1858e-9,
357  -1.273e-10,3.3e-12,-1e-13 };
358 
359  const double c__[14] = { 1.9405491464836,.9413409132865,
360  -.579845034293,.3091572011159,-.0916101792208,.0164437407515,
361  -.0019713091952,1.692538851e-4,-1.09393296e-5,5.522386e-7,
362  -2.23995e-8,7.465e-10,-2.08e-11,5e-13 };
363 
364  const double p[23] = { .96074783975204,-.0371138962124,
365  .00194143988899,-1.7165988425e-4,2.112637753e-5,-3.27163257e-6,
366  6.0069212e-7,-1.2586794e-7,2.932563e-8,-7.45696e-9,2.04105e-9,
367  -5.9502e-10,1.8323e-10,-5.921e-11,1.997e-11,-7e-12,2.54e-12,
368  -9.5e-13,3.7e-13,-1.4e-13,6e-14,-2e-14,1e-14 };
369  const double q[20] = { .98604065696238,-.0134717382083,
370  4.5329284117e-4,-3.067288652e-5,3.13199198e-6,-4.2110196e-7,
371  6.907245e-8,-1.318321e-8,2.83697e-9,-6.7329e-10,1.734e-10,
372  -4.787e-11,1.403e-11,-4.33e-12,1.4e-12,-4.7e-13,1.7e-13,-6e-14,
373  2e-14,-1e-14 };
374 
375  // System generated locals
376  double d__1;
377 
378  // Local variables
379  double h__;
380  int i__;
381  double r__, y, b0, b1, b2, pp, qq, alfa;
382 
383  sint=0;
384  cint=0;
385 
386 
387  if (fabs(x) <= eight) {
388  y = x / eight;
389  // Computing 2nd power
390  d__1 = y;
391  h__ = two * (d__1 * d__1) - one;
392  alfa = -two * h__;
393 
394  // cos
395  if (x!=0) {
396  b1 = zero;
397  b2 = zero;
398  for (i__ = 13; i__ >= 0; --i__) {
399  b0 = c__[i__] - alfa * b1 - b2;
400  b2 = b1;
401  b1 = b0;
402  }
403  cint = ce + log((fabs(x))) - b0 + h__ * b2;
404  }
405  // sin
406  b1 = zero;
407  b2 = zero;
408  for (i__ = 13; i__ >= 0; --i__) {
409  b0 = s__[i__] - alfa * b1 - b2;
410  b2 = b1;
411  b1 = b0;
412  }
413  sint = y * (b0 - b2);
414 
415  } else {
416  r__ = one / x;
417  y = eight * r__;
418  // Computing 2nd power
419  d__1 = y;
420  h__ = two * (d__1 * d__1) - one;
421  alfa = -two * h__;
422  b1 = zero;
423  b2 = zero;
424  for (i__ = 22; i__ >= 0; --i__) {
425  b0 = p[i__] - alfa * b1 - b2;
426  b2 = b1;
427  b1 = b0;
428  }
429  pp = b0 - h__ * b2;
430  b1 = zero;
431  b2 = zero;
432  for (i__ = 19; i__ >= 0; --i__) {
433  b0 = q[i__] - alfa * b1 - b2;
434  b2 = b1;
435  b1 = b0;
436  }
437  qq = b0 - h__ * b2;
438  // cos
439  cint = r__ * (qq * sin(x) - r__ * pp * cos(x));
440  // sin
441  d__1 = pih;
442  if(x < 0.) d__1 = -d__1;
443  sint = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
444  }
445  }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
double VVIObjDetails::sinint ( double  x)

Private version of the cosine integral.

Definition at line 271 of file VVIObj.cc.

References funct::cos(), MillePedeFileConverter_cfg::e, AlCaHLTBitMon_ParallelJobs::p, createTree::pp, lumiQueryAPI::q, alignCSCRings::s, funct::sin(), x, and y.

271  {
272  // Initialized data
273 
274  const double zero = 0.;
275  const double one = 1.;
276  const double two = 2.;
277  const double eight = 8.;
278  const double pih = 1.5707963267949;
279  const double s[14] = { 1.9522209759531,-.6884042321257,
280  .4551855132256,-.1804571236838,.0410422133759,-.0059586169556,
281  6.001427414e-4,-4.44708329e-5,2.5300782e-6,-1.141308e-7,4.1858e-9,
282  -1.273e-10,3.3e-12,-1e-13 };
283  const double p[23] = { .96074783975204,-.0371138962124,
284  .00194143988899,-1.7165988425e-4,2.112637753e-5,-3.27163257e-6,
285  6.0069212e-7,-1.2586794e-7,2.932563e-8,-7.45696e-9,2.04105e-9,
286  -5.9502e-10,1.8323e-10,-5.921e-11,1.997e-11,-7e-12,2.54e-12,
287  -9.5e-13,3.7e-13,-1.4e-13,6e-14,-2e-14,1e-14 };
288  const double q[20] = { .98604065696238,-.0134717382083,
289  4.5329284117e-4,-3.067288652e-5,3.13199198e-6,-4.2110196e-7,
290  6.907245e-8,-1.318321e-8,2.83697e-9,-6.7329e-10,1.734e-10,
291  -4.787e-11,1.403e-11,-4.33e-12,1.4e-12,-4.7e-13,1.7e-13,-6e-14,
292  2e-14,-1e-14 };
293 
294  // System generated locals
295  double d__1;
296 
297  // Local variables
298  double h__;
299  int i__;
300  double r__, y, b0, b1, b2, pp, qq, alfa;
301 
302  if (fabs(x) <= eight) {
303  y = x / eight;
304  d__1 = y;
305  h__ = two * (d__1 * d__1) - one;
306  alfa = -two * h__;
307  b1 = zero;
308  b2 = zero;
309  for (i__ = 13; i__ >= 0; --i__) {
310  b0 = s[i__] - alfa * b1 - b2;
311  b2 = b1;
312  b1 = b0;
313  }
314  b1 = y * (b0 - b2);
315  } else {
316  r__ = one / x;
317  y = eight * r__;
318  d__1 = y;
319  h__ = two * (d__1 * d__1) - one;
320  alfa = -two * h__;
321  b1 = zero;
322  b2 = zero;
323  for (i__ = 22; i__ >= 0; --i__) {
324  b0 = p[i__] - alfa * b1 - b2;
325  b2 = b1;
326  b1 = b0;
327  }
328  pp = b0 - h__ * b2;
329  b1 = zero;
330  b2 = zero;
331  for (i__ = 19; i__ >= 0; --i__) {
332  b0 = q[i__] - alfa * b1 - b2;
333  b2 = b1;
334  b1 = b0;
335  }
336  qq = b0 - h__ * b2;
337  d__1 = fabs(pih);
338  if(x < 0.) d__1 = -d__1;
339  b1 = d__1 - r__ * (r__ * pp * sin(x) + qq * cos(x));
340  }
341 
342  return b1;
343  } // sinint
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Cos< T >::type cos(const T &t)
Definition: Cos.h:22