| File: | libraries/AMPTOOLS_AMPS/dblReggeMod.cc |
| Location: | line 303, column 2 |
| 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 "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; |
| 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; |
Value stored to 'A' is never read | |
| 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 | } |