CMS 3D CMS Logo

Functions
sistripvvi::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)
 
double expint (double x)
 Private version of the sine integral. More...
 
double f1 (double x, double const *h_)
 Private version of the exponential integral. More...
 
double f2 (double x, double const *h_)
 
void sincosint (double x, double &sint, double &cint)
 
double sinint (double x)
 Private version of the cosine integral. More...
 

Function Documentation

double sistripvvi::VVIObjDetails::cosint ( double  x)

Private version of the cosine and sine integral.

Definition at line 197 of file VVIObj.cc.

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

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

Definition at line 574 of file VVIObj.cc.

References f, f1(), 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 f2(), and sistripvvi::VVIObj::VVIObj().

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

Private version of the sine integral.

Definition at line 452 of file VVIObj.cc.

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

Referenced by f2(), and sistripvvi::VVIObj::VVIObj().

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

Private version of the exponential integral.

Definition at line 35 of file VVIObj.cc.

References cmsBatch::log.

Referenced by VVIObjDetails::dzero(), dzero(), VVIObj::VVIObj(), and sistripvvi::VVIObj::VVIObj().

35 { return h_[0]+h_[1]*std::log(h_[2]*x)-h_[3]*x;}
double sistripvvi::VVIObjDetails::f2 ( double  x,
double const *  h_ 
)
inline

Definition at line 36 of file VVIObj.cc.

References a, funct::abs(), b, dzero(), JetChargeProducer_cfi::exp, expint(), F(), jets_cff::func, and cmsBatch::log.

Referenced by VVIObjDetails::dzero(), dzero(), VVIObj::VVIObj(), and sistripvvi::VVIObj::VVIObj().

36 { return h_[4]-x+h_[5]*(std::log(std::abs(x))+expint(x))-h_[6]*std::exp(-x);}
double expint(double x)
Private version of the sine integral.
Definition: VVIObj.cc:452
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void sistripvvi::VVIObjDetails::sincosint ( double  x,
double &  sint,
double &  cint 
)

Definition at line 349 of file VVIObj.cc.

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

Referenced by sistripvvi::VVIObj::VVIObj().

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

Private version of the cosine integral.

Definition at line 275 of file VVIObj.cc.

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

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