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