CMS 3D CMS Logo
HankelFunction.cc
Go to the documentation of this file.
1 /*
2 
3 Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
4 amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
5 November. 2, 2005
6 
7 */
8 
9 #include <TMath.h>
10 
11 #include "GeneratorInterface/Hydjet2Interface/interface/HankelFunction.h"
12 
13 //compute Hankel function of zeroth order
14 enum { kNe = 2, kNCoeff = 9 };
15 
16 const double i0Coefficient[kNCoeff][kNe] = {
17  //coefficients to compute function I0
18  {1.0, 0.39894228},
19  {3.5156229, 0.01328592},
20  {3.0899424, 0.00225319},
21  {1.2067492, -0.00157565},
22  {0.2659732, 0.00916281},
23  {0.0360768, -0.02057706},
24  {0.0045813, 0.02635537},
25  {0., -0.01647633},
26  {0., 0.00392377}};
27 
28 const double i1Coefficient[kNCoeff][kNe] = {
29  //coefficients to compute function I1
30  {0.5, 0.39894228},
31  {0.87890594, -0.03988024},
32  {0.51498869, -0.00362018},
33  {0.15084934, 0.00163801},
34  {0.02658733, -0.01031555},
35  {0.00301532, 0.02282967},
36  {0.00032411, -0.02895312},
37  {0., 0.01787654},
38  {0., -0.00420059}};
39 
40 const double k0Coefficient[kNCoeff][kNe] = {
41  //coefficients to compute modified Hankel function of the zeroth order K0
42  {-0.57721566, 1.25331414},
43  {0.42278420, -0.07832358},
44  {0.23069756, 0.02189568},
45  {0.03488590, -0.01062446},
46  {0.00262698, 0.00587872},
47  {0.00010750, -0.00251540},
48  {0.0000074, 0.00053208},
49  {0., 0.},
50  {0., 0.}};
51 
52 const double k1Coefficient[kNCoeff][kNe] = {
53  //coefficients to compute modified Hankel function of the first order K1
54  {1.0, 1.25331414},
55  {0.15443144, 0.23498619},
56  {-0.67278579, -0.03655620},
57  {-0.18156897, 0.01504268},
58  {-0.01919402, -0.00780353},
59  {-0.00110404, 0.00325614},
60  {-0.00004686, -0.00068245},
61  {0., 0.},
62  {0., 0.}};
63 
64 double BesselI0(double x) {
65  // (C) Copr. 1986-92 Numerical Recipes Software +7.
66  //compute Bessel function of zeroth order
67 
68  const double p1 = i0Coefficient[0][0];
69  const double p2 = i0Coefficient[1][0];
70  const double p3 = i0Coefficient[2][0];
71  const double p4 = i0Coefficient[3][0];
72  const double p5 = i0Coefficient[4][0];
73  const double p6 = i0Coefficient[5][0];
74  const double p7 = i0Coefficient[6][0];
75 
76  const double q1 = i0Coefficient[0][1];
77  const double q2 = i0Coefficient[1][1];
78  const double q3 = i0Coefficient[2][1];
79  const double q4 = i0Coefficient[3][1];
80  const double q5 = i0Coefficient[4][1];
81  const double q6 = i0Coefficient[5][1];
82  const double q7 = i0Coefficient[6][1];
83  const double q8 = i0Coefficient[7][1];
84  const double q9 = i0Coefficient[8][1];
85 
86  double i0 = 0.;
87 
88  if (TMath::Abs(x) < 3.75) {
89  double y = (x / 3.75) * (x / 3.75);
90  i0 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
91  } else {
92  double ax = TMath::Abs(x);
93  double y = 3.75 / ax;
94  i0 = (TMath::Exp(ax) / TMath::Sqrt(ax)) *
95  (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
96  }
97 
98  return i0;
99 }
100 
101 double BesselI1(double x) {
102  // (C) Copr. 1986-92 Numerical Recipes Software +7.
103  //compute Bessel function of first order
104 
105  const double p1 = i1Coefficient[0][0];
106  const double p2 = i1Coefficient[1][0];
107  const double p3 = i1Coefficient[2][0];
108  const double p4 = i1Coefficient[3][0];
109  const double p5 = i1Coefficient[4][0];
110  const double p6 = i1Coefficient[5][0];
111  const double p7 = i1Coefficient[6][0];
112 
113  const double q1 = i1Coefficient[0][1];
114  const double q2 = i1Coefficient[1][1];
115  const double q3 = i1Coefficient[2][1];
116  const double q4 = i1Coefficient[3][1];
117  const double q5 = i1Coefficient[4][1];
118  const double q6 = i1Coefficient[5][1];
119  const double q7 = i1Coefficient[6][1];
120  const double q8 = i1Coefficient[7][1];
121  const double q9 = i1Coefficient[8][1];
122 
123  double i1 = 0.;
124 
125  if (TMath::Abs(x) < 3.75) {
126  double y = (x / 3.75) * (x / 3.75);
127  i1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
128  } else {
129  double ax = TMath::Abs(x);
130  double y = 3.75 / ax;
131  i1 = (TMath::Exp(ax) / TMath::Sqrt(ax)) *
132  (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
133  if (x < 0.)
134  i1 = -i1;
135  }
136 
137  return i1;
138 }
139 
140 double HankelK0(double x) {
141  const double p1 = k0Coefficient[0][0];
142  const double p2 = k0Coefficient[1][0];
143  const double p3 = k0Coefficient[2][0];
144  const double p4 = k0Coefficient[3][0];
145  const double p5 = k0Coefficient[4][0];
146  const double p6 = k0Coefficient[5][0];
147  const double p7 = k0Coefficient[6][0];
148 
149  const double q1 = k0Coefficient[0][1];
150  const double q2 = k0Coefficient[1][1];
151  const double q3 = k0Coefficient[2][1];
152  const double q4 = k0Coefficient[3][1];
153  const double q5 = k0Coefficient[4][1];
154  const double q6 = k0Coefficient[5][1];
155  const double q7 = k0Coefficient[6][1];
156 
157  double k0 = 0.;
158  if (x <= 2.0) {
159  double y = x * x / 4.0;
160  k0 = (-TMath::Log(x / 2.0) * BesselI0(x)) + (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
161  } else {
162  double y = (2.0 / x);
163  k0 = (TMath::Exp(-x) / TMath::Sqrt(x)) * (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
164  }
165 
166  return k0;
167 }
168 
169 double HankelK1(double x) {
170  // (C) Copr. 1986-92 Numerical Recipes Software +7.
171  // compute modified Hankel function of the first order
172  const double p1 = k1Coefficient[0][0];
173  const double p2 = k1Coefficient[1][0];
174  const double p3 = k1Coefficient[2][0];
175  const double p4 = k1Coefficient[3][0];
176  const double p5 = k1Coefficient[4][0];
177  const double p6 = k1Coefficient[5][0];
178  const double p7 = k1Coefficient[6][0];
179 
180  const double q1 = k1Coefficient[0][1];
181  const double q2 = k1Coefficient[1][1];
182  const double q3 = k1Coefficient[2][1];
183  const double q4 = k1Coefficient[3][1];
184  const double q5 = k1Coefficient[4][1];
185  const double q6 = k1Coefficient[5][1];
186  const double q7 = k1Coefficient[6][1];
187 
188  double k1 = 0.;
189 
190  if (x <= 2.0) {
191  double y = x * x / 4.0;
192  k1 = (TMath::Log(x / 2.0) * BesselI1(x)) +
193  (1.0 / x) * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
194  } else {
195  double y = 2.0 / x;
196  k1 = (TMath::Exp(-x) / TMath::Sqrt(x)) * (q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
197  }
198 
199  return k1;
200 }
201 
202 double HankelKn(int n, double x) {
203  // (C) Copr. 1986-92 Numerical Recipes Software +7.
204  // compute modified Hankel function of the first order
205  if (n < 2)
206  throw "Bad argument n in Kn";
207 
208  double tox = 2.0 / x;
209  double km = HankelK0(x);
210  double k = HankelK1(x);
211  double kp = 0.;
212 
213  for (int c = 1; c <= n - 1; ++c) {
214  kp = km + c * tox * k;
215  km = k;
216  k = kp;
217  }
218 
219  return k;
220 }
k1Coefficient
const double k1Coefficient[kNCoeff][kNe]
Definition: HankelFunction.cc:52
HankelK0
double HankelK0(double x)
Definition: HankelFunction.cc:140
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
reco::ParticleMasses::k0
const double k0
Definition: ParticleMasses.h:7
testProducerWithPsetDescEmpty_cfi.i1
i1
Definition: testProducerWithPsetDescEmpty_cfi.py:45
BesselI1
double BesselI1(double x)
Definition: HankelFunction.cc:101
HankelKn
double HankelKn(int n, double x)
Definition: HankelFunction.cc:202
Abs
T Abs(T a)
Definition: MathUtil.h:49
p2
double p2[4]
Definition: TauolaWrapper.h:90
dqmdumpme.k
k
Definition: dqmdumpme.py:60
q2
double q2[4]
Definition: TauolaWrapper.h:88
q1
double q1[4]
Definition: TauolaWrapper.h:87
HankelFunction.h
kNe
Definition: HankelFunction.cc:14
p4
double p4[4]
Definition: TauolaWrapper.h:92
p1
double p1[4]
Definition: TauolaWrapper.h:89
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
BesselI0
double BesselI0(double x)
Definition: HankelFunction.cc:64
genVertex_cff.x
x
Definition: genVertex_cff.py:12
detailsBasic3DVector::y
float float y
Definition: extBasic3DVector.h:14
i1Coefficient
const double i1Coefficient[kNCoeff][kNe]
Definition: HankelFunction.cc:28
kNCoeff
Definition: HankelFunction.cc:14
p3
double p3[4]
Definition: TauolaWrapper.h:91
k0Coefficient
const double k0Coefficient[kNCoeff][kNe]
Definition: HankelFunction.cc:40
HankelK1
double HankelK1(double x)
Definition: HankelFunction.cc:169
i0Coefficient
const double i0Coefficient[kNCoeff][kNe]
Definition: HankelFunction.cc:16
kp
int kp
Definition: CascadeWrapper.h:13
Generated for CMSSW Reference Manual by doxygen 1.8.16