File: | libraries/AMPTOOLS_AMPS/dblReggeMod.cc |
Location: | line 206, column 3 |
Description: | Value stored to 'y' is never read |
1 | #include <cassert> |
2 | #include <complex> |
3 | #include <iostream> |
4 | #include <string> |
5 | #include <sstream> |
6 | #include <cstdlib> |
7 | #include <math.h> |
8 | #include "TLorentzVector.h" |
9 | #include "TLorentzRotation.h" |
10 | |
11 | #include "IUAmpTools/Kinematics.h" |
12 | #include "dblReggeMod.h" |
13 | |
14 | dblReggeMod::dblReggeMod( const vector< string >& args ) : |
15 | UserAmplitude< dblReggeMod >( args ) |
16 | { |
17 | assert( args.size() == 5 )((args.size() == 5) ? static_cast<void> (0) : __assert_fail ("args.size() == 5", "libraries/AMPTOOLS_AMPS/dblReggeMod.cc" , 17, __PRETTY_FUNCTION__)); |
18 | // j = atoi( args[0].c_str() ); |
19 | // fast = atoi( args[1].c_str() ); |
20 | a_eta = AmpParameter(args[0]); |
21 | a_pi = AmpParameter(args[1]); |
22 | b_eta = AmpParameter( args[2] ); |
23 | b_pi = AmpParameter( args[3]); |
24 | S0 = AmpParameter( args[4] ); |
25 | |
26 | registerParameter( a_pi ); |
27 | registerParameter( a_eta ); |
28 | registerParameter( b_pi ); |
29 | registerParameter( b_eta ); |
30 | registerParameter( S0 ); |
31 | } |
32 | |
33 | complex< GDouble > |
34 | dblReggeMod::calcAmplitude( GDouble** pKin ) const { |
35 | TLorentzVector beam ( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] ); |
36 | TLorentzVector recoil ( pKin[1][1], pKin[1][2], pKin[1][3], pKin[1][0] ); |
37 | TLorentzVector p1 ( pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0] ); |
38 | TLorentzVector p2 ( pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0] ); |
39 | |
40 | TLorentzVector resonance = p1 + p2; |
41 | |
42 | double mass2[4] = {beam.M2(), p1.M2(), p2.M2(), recoil.M2()}; |
43 | |
44 | double s12 = resonance.M2(); |
45 | //double s13 = (p1 + recoil).M2(); |
46 | double s23 = (p2 + recoil).M2(); |
47 | double t1 = (beam - p1).M2(); |
48 | double t2 = (beam - p2).M2(); |
49 | double s = (recoil + p1 + p2).M2(); |
50 | double u3 = t1 + t2 +s12 - (beam.M2() + p1.M2() + p2.M2()); |
51 | |
52 | |
53 | double param[4] = {a_eta,a_pi,b_eta,b_pi}; |
54 | double inv[5] = {s, s12, s23, t1 ,u3}; |
55 | // double varGJ[5] = {s,s12,u3, theta, phi}; |
56 | |
57 | int hel[3] = {1,-1,-1}; |
58 | std::complex<double> amp = ampEtaPi0(param, hel, inv, mass2); |
59 | |
60 | return amp; |
61 | } |
62 | |
63 | |
64 | std::complex<double> dblReggeMod::ampEtaPi0(double par[4], int hel[3], double inv[5], double mass2[4]) const{ |
65 | |
66 | std::complex<double> zero (0,0); |
67 | |
68 | double s,s12,s23,t1,u3; |
69 | double m12, m22, m32, ma2, mb2; |
70 | s = inv[0]; s12 = inv[1]; s23 = inv[2]; t1 = inv[3]; u3 = inv[4]; |
71 | ma2 = mass2[0]; m12 = mass2[1]; m22 = mass2[2]; m32 = mass2[3]; |
72 | double t2 = -t1+u3-s12+ma2+m12+m22; |
73 | double s13 = s-s12-s23+m12+m22+m32; |
74 | |
75 | // scalar part |
76 | double app = 0.9; // slope of Regge trajectories alpha' |
77 | double alp0eta = app*t1 + 0.5; |
78 | double alp0pi0 = app*t2 + 0.5; |
79 | double alp1 = app*u3 + 0.5; |
80 | |
81 | |
82 | int tau[2] = {-1, -1}; // only vector exchange |
83 | double si[2] = {s12,s23}; |
84 | double alp[2] = {alp0eta, alp1}; |
85 | |
86 | |
87 | std::complex<double> ADR1 = DoubleRegge(tau, s, si, alp); // fast eta |
88 | |
89 | |
90 | si[1] = s13; alp[0] = alp0pi0; |
91 | |
92 | std::complex<double> ADR2 = DoubleRegge(tau, s, si, alp); // fast pi0 |
93 | |
94 | // helicity part |
95 | double fac1 = sqrt(-t1/mass2[2]); |
96 | double fac2 = sqrt(-t2/mass2[2]); // use the pion mass in both fac1 and fac2 |
97 | double fac3 = pow(-u3/4./mass2[0],abs((hel[1]-hel[2])/4.)); // hel[1,2] are twice the nucleon helicities! |
98 | double parity = pow(-1,(hel[1]-hel[2])/2.); |
99 | if(hel[1] == -1){fac3 = fac3*parity;} |
100 | |
101 | double Bot1 = a_eta*exp(b_eta*t1); |
102 | double Bot2 = a_pi*exp(b_pi*t2); |
103 | |
104 | return fac3*(Bot1*fac1*ADR1 + Bot2*fac2*ADR2 ); |
105 | } |
106 | |
107 | std::complex<double> dblReggeMod::V12(double alp1, double alp2, double eta) const{ |
108 | |
109 | if(alp1==alp2 ){return 0.0;} |
110 | std::complex<double> res = CHGM(-alp1, 1.-alp1+alp2, -1/eta); |
111 | res *= cgamma(alp1-alp2,0)/cgamma(-alp2,0); |
112 | |
113 | return res; |
114 | } |
115 | |
116 | std::complex<double> dblReggeMod::DoubleRegge(int tau[2], double s, double si[2], double alp[2]) const{ |
117 | std::complex<double> ui (0,1); |
118 | // signature factors: |
119 | std::complex<double> x0 = 1/2.*((double)tau[0] + exp(-ui*M_PI3.14159265358979323846*alp[0])); |
120 | std::complex<double> x1 = 1/2.*((double)tau[1] + exp(-ui*M_PI3.14159265358979323846*alp[1])); |
121 | std::complex<double> x01 = 1/2.*((double)tau[0]*tau[1] + exp(-ui*M_PI3.14159265358979323846*(alp[0]-alp[1]))); |
122 | std::complex<double> x10 = 1/2.*((double)tau[1]*tau[0] + exp(-ui*M_PI3.14159265358979323846*(alp[1]-alp[0]))); |
123 | // double Regge vertices: |
124 | double eta = S0*s/(si[0]*si[1]); |
125 | std::complex<double> V0 = V12(alp[0], alp[1], eta); |
126 | std::complex<double> V1 = V12(alp[1], alp[0], eta); |
127 | |
128 | // combine pieces: |
129 | std::complex<double> t1 = pow(s/S0,alp[1])*pow(si[0]/S0,alp[0]-alp[1])*x1*x01*V1; |
130 | std::complex<double> t0 = pow(s/S0,alp[0])*pow(si[1]/S0,alp[1]-alp[0])*x0*x10*V0; |
131 | |
132 | |
133 | return (t0+t1)*cgamma(-alp[0],0)*cgamma(-alp[1],0);; |
134 | } |
135 | |
136 | |
137 | std::complex<double> dblReggeMod::cgamma(std::complex<double> z,int OPT) const{ |
138 | std::complex<double> ui (0,1); |
139 | std::complex<double> g, infini= 1e308+ 0.0*ui; // z0,z1 |
140 | double x0,q1,q2,x,y,th,th1,th2,g0,gr,gi,gr1,gi1; |
141 | double na=0.0,t,x1 = 1,y1=0.0,sr,si; |
142 | int j,k; |
143 | |
144 | |
145 | static double a[] = { |
146 | 8.333333333333333e-02, |
147 | -2.777777777777778e-03, |
148 | 7.936507936507937e-04, |
149 | -5.952380952380952e-04, |
150 | 8.417508417508418e-04, |
151 | -1.917526917526918e-03, |
152 | 6.410256410256410e-03, |
153 | -2.955065359477124e-02, |
154 | 1.796443723688307e-01, |
155 | -1.39243221690590}; |
156 | |
157 | x = real(z); |
158 | y = imag(z); |
159 | |
160 | if (x > 171) return infini; |
161 | if ((y == 0.0) && (x == (int)x) && (x <= 0.0)) |
162 | return infini; |
163 | else if (x < 0.0) { |
164 | x1 = x; |
165 | y1 = y; |
166 | x = -x; |
167 | y = -y; |
168 | } |
169 | x0 = x; |
170 | if (x <= 7.0) { |
171 | na = (int)(7.0-x); |
172 | x0 = x+na; |
173 | } |
174 | q1 = sqrt(x0*x0+y*y); |
175 | th = atan(y/x0); |
176 | gr = (x0-0.5)*log(q1)-th*y-x0+0.5*log(2.0*M_PI3.14159265358979323846); |
177 | gi = th*(x0-0.5)+y*log(q1)-y; |
178 | for (k=0;k<10;k++){ |
179 | t = pow(q1,-1.0-2.0*k); |
180 | gr += (a[k]*t*cos((2.0*k+1.0)*th)); |
181 | gi -= (a[k]*t*sin((2.0*k+1.0)*th)); |
182 | } |
183 | if (x <= 7.0) { |
184 | gr1 = 0.0; |
185 | gi1 = 0.0; |
186 | for (j=0;j<na;j++) { |
187 | gr1 += (0.5*log((x+j)*(x+j)+y*y)); |
188 | gi1 += atan(y/(x+j)); |
189 | } |
190 | gr -= gr1; |
191 | gi -= gi1; |
192 | } |
193 | |
194 | |
195 | if (x1 <= 0.0) { |
196 | q1 = sqrt(x*x+y*y); |
197 | th1 = atan(y/x); |
198 | sr = -sin(M_PI3.14159265358979323846*x)*cosh(M_PI3.14159265358979323846*y); |
199 | si = -cos(M_PI3.14159265358979323846*x)*sinh(M_PI3.14159265358979323846*y); |
200 | q2 = sqrt(sr*sr+si*si); |
201 | th2 = atan(si/sr); |
202 | if (sr < 0.0) th2 += M_PI3.14159265358979323846; |
203 | gr = log(M_PI3.14159265358979323846/(q1*q2))-gr; |
204 | gi = -th1-th2-gi; |
205 | x = x1; |
206 | y = y1; |
Value stored to 'y' is never read | |
207 | } |
208 | |
209 | if (OPT == 0) { |
210 | g0 = exp(gr); |
211 | gr = g0*cos(gi); |
212 | gi = g0*sin(gi); |
213 | } |
214 | g = gr + ui*gi; |
215 | |
216 | return g; |
217 | } |
218 | |
219 | |
220 | |
221 | double dblReggeMod::CHGM(double A, double B, double X) const{ |
222 | double A0=A, A1=A, X0=X, HG = 0.0; |
223 | double TBA, TB, TA, Y0=0.0, Y1=0.0, RG, LA = (int) A, NL, R, M, INF = pow(10,300); |
224 | double sum1, sum2, R1, R2, HG1, HG2; |
225 | if (B == 0.0 || B == -abs( (int) B)){ |
226 | HG = INF; |
227 | } else if(A == 0.0 || X == 0.0) { |
228 | HG = 1.0; |
229 | } else if(A == -1.0){ |
230 | HG = 1.0 - X/B; |
231 | } else if(A == B){ |
232 | HG = exp(X); |
233 | } else if (A-B == 1.0){ |
234 | HG = (1.0+X/B)*exp(X); |
235 | } else if (A == 1.0 && B == 2.0){ |
236 | HG = (exp(X)-1.0)/X; |
237 | } else if(A == (int)A && A < 0.0){ |
238 | M = (int) -A; |
239 | R = 1.0; |
240 | HG = 1.0; |
241 | for (int k = 1; k<= M ; k++) { |
242 | R = R*(A+k-1.0)/k/(B+k-1.0)*X; |
243 | HG+=R; |
244 | } |
245 | } |
246 | if(HG != 0){return HG;} |
247 | |
248 | if(X<0.0){ |
249 | A = B-A; |
250 | A0 = A; |
251 | X = fabs(X); |
252 | } |
253 | if(A<2.0) {NL = 0;} |
254 | else{ |
255 | NL = 1; |
256 | LA = (int) A; |
257 | A = A-LA-1.0; |
258 | } |
259 | for (int n = 0; n<= NL; n++) { |
260 | if(A0 >= 2.0 ) { A+=1.0; } |
261 | if(X <= 30.0 + fabs(B) || A < 0.0){ |
262 | HG = 1.0; |
263 | RG = 1.0; |
264 | for (int j = 1; j<= 500; j++) { |
265 | RG = RG*(A+j-1)/(j*(B+j-1))*X; |
266 | HG += RG; |
267 | if(fabs(RG/HG) < pow(10.,-15.)) { |
268 | if(n==0) {Y0 = HG;} |
269 | if(n==1) {Y1 = HG;} |
270 | } |
271 | continue; |
272 | } |
273 | } else { |
274 | TA = tgamma(A); |
275 | TB = tgamma(B); |
276 | TBA = tgamma(B-A); |
277 | sum1 = 1.0; |
278 | sum2 = 1.0; |
279 | R1 = 1.0; |
280 | R2 = 1.0; |
281 | for (int i = 1; i<=8; i++) { |
282 | R1 = - R1*(A+i-1)*(A-B+i)/(X*i); |
283 | R2 = - R2*(B-A+i-1)*(A-i)/(X*i); |
284 | sum1+=R1; |
285 | sum2+=R2; |
286 | } |
287 | HG1 = TB/TBA*pow(X,-A)*cos(M_PI3.14159265358979323846*A)*sum1; |
288 | HG2 = TB/TA*exp(X)*pow(X,A-B)*sum2; |
289 | HG = HG1+HG2; |
290 | } |
291 | if(n==0) {Y0 = HG;} |
292 | if(n==1) {Y1 = HG;} |
293 | } |
294 | if(A0 >= 2.0){ |
295 | for (int i=1; i<=LA-1; i++) { |
296 | HG = ((2.*A-B+X)*Y1+(B-A)*Y0)/A; |
297 | Y0 = Y1; |
298 | Y1 = HG; |
299 | A += 1.; |
300 | } |
301 | } |
302 | if(X0<0.0) {HG = HG*exp(X0);} |
303 | A = A1; |
304 | X = X0; |
305 | |
306 | |
307 | return HG; |
308 | } |
309 | |
310 | void |
311 | dblReggeMod::updatePar( const AmpParameter& par ){ |
312 | |
313 | // could do expensive calculations here on parameter updates |
314 | // |
315 | } |