| File: | libraries/AMPTOOLS_AMPS/DblRegge_FastEta.cc |
| Location: | line 325, column 9 |
| Description: | Value stored to 'X' 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; |
| 325 | X = X0; |
Value stored to 'X' is never read | |
| 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 |