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