Bug Summary

File:libraries/AMPTOOLS_AMPS/Pi0ReggeModel.cc
Location:line 368, column 5
Description:Value stored to 'x' is never read

Annotated Source Code

1/*
2 ============================================================================
3 Name : Pi0ReggeModel.cc
4 Author : Vincent Mathieu (Adapted for C++ and GlueX sim-recon by Justin Stevens)
5 Version : v1.0 June 2015
6 Copyright : MyLab
7 Publication : arXiv:1505.02321
8 Description : Photoproduction of a pseudoscalar : gamma + N --> 0- + N'
9 Applied to gamma + p --> pi0 + p
10 To change the model, change 'CGLN_Ai'
11 ============================================================================
12*/
13
14#include "Pi0ReggeModel.h"
15
16// *********************************************************************************
17double Pi0PhotCS_S(double E, double theta, double &BeamSigma){
18/*
19 * compute the differential cross section from s-channel helicities
20 * in micro barns/Gev^2
21 */
22
23 double MP = 0.938272046;
24 double MPI = 0.1349766;
25 double mass[5] = {0.0, 0.0, MP, MPI, MP};
26 double pa[4], pb[4], pc[4], pd[4];
27 int hel[4] = {2,1,0,1}; // hel = {1,+,0,+} (x2)
28 complex<double> res[4] ;
29 double sig;
30 kin2to2(E, theta, mass, pa, pb, pc, pd);
31
32 res[0] = Pi0PhotAmpS(pa, pb, pc, hel); // hel = {1,+,0,+} (x2)
33
34 hel[1] = -1; // hel = {1,-,0,+} (x2)
35 res[1] = Pi0PhotAmpS(pa, pb, pc, hel);
36
37 hel[3] = -1; // hel = {1,-,0,-} (x2)
38 res[2] = Pi0PhotAmpS(pa, pb, pc, hel);
39
40 hel[1] = +1; // hel = {1,+,0,-} (x2)
41 res[3] = Pi0PhotAmpS(pa, pb, pc, hel);
42
43 sig = pow(abs(res[0]),2) + pow(abs(res[1]),2) + pow(abs(res[2]),2) + pow(abs(res[3]),2);
44 sig = sig/32 * 389.3;
45 //sig = sig/32/M_PI/(E*E-MP*MP)/(E*E-MP*MP) * 389.3;
46
47 // Sigma beam asymmetry from equation B3b from paper
48 BeamSigma = real(res[0]*conj(res[2]) - res[1]*conj(res[3]));
49 BeamSigma = BeamSigma/16 * 389.3;
50 BeamSigma = BeamSigma/sig;
51
52 return sig;
53}
54
55// *********************************************************************************
56complex<double> Pi0PhotAmpS(double pa[],double pb[],double pc[], int hel[]){
57 /* Regge amplitudes for gamma + p --> pi0 + p
58 * See arXiv:1505.02321 for the formulas
59 * Amplitudes in the S-CHANNEL !
60 * Helicities are defined in the c.o.m. of (gamma,p)
61 * Inputs:
62 * pa, pb, pc : four vectors in the c.o.m. frame.
63 * hel={mua, mub, muc, mud} : 2 x helicities of a,b,c,d
64 * Outputs:
65 * ampl : the amplitude ; complex number
66 */
67 complex<double> amp = 0.0;
68 complex<double> Ai[5] = {0}, Fi[5] = {0};
69 //pbc set but not used: double pab[4], pbc[4], pca[4], pd[4];
70 double pab[4], pca[4], pd[4];
71 struct Kin var;
72 double mass[5]={0};
73
74 // Only valid for Q^2 = 0 --> photon helicity should be +1 or -1
75 if ( abs(hel[0]) != 2 ) return 0.0;
76
77 int i;
78 for(i=0;i<4;i++){ // compute the Mandelstam invariant
79 pab[i] = pa[i] + pb[i]; // from four-vectors
80 // pbc set but not used: pbc[i] = pb[i] - pc[i]; //
81 pca[i] = pc[i] - pa[i]; //
82 pd[i] = pa[i] + pb[i] - pc[i]; //
83 }
84
85 double EPS = 0.0;
86 complex<double> I(0,1);
87 var.s = snorm( pab ) + I*EPS;
88 var.t = snorm( pca ) - I*EPS;
89 // photon should be real
90 mass[1] = 0.0;
91 mass[2] = sqrt( snorm( pb ) );
92 mass[3] = sqrt( snorm( pc ) );
93 mass[4] = sqrt( snorm( pd ) );
94
95 kinematics(var.s, var.t, mass, &var); // fill the structure var with all kin. quantities
96 CGLN_Ai(var.s,var.t,Ai); // call the model : CHANGE HERE FOR ANOTHER MODEL
97
98 CGLNA2F(var, Ai, Fi); // Convert CGLN Ai to CGLN Fi
99
100 // Test nucleon helicities hel[1] and hel[3]
101 if ( hel[1] == 1 && hel[3] == 1 ) {
102 // hel = {1,+,0,+}
103 amp = sqrt(2) * var.sinsh * (Fi[2] + Fi[1] )
104 + 1/sqrt(2) * var.sins * var.cossh * (Fi[3] + Fi[4]);
105 }
106 else if ( hel[1] == -1 && hel[3] == -1 ){
107 // hel = {1,-,0,-}
108 amp = -1/sqrt(2) * var.sins * var.cossh * (Fi[3] + Fi[4]);
109 }
110 else if ( hel[1] == -1 && hel[3] == 1 ){
111 // hel = {1,-,0,+}
112 amp = sqrt(2) * var.cossh * (Fi[2] - Fi[1] )
113 + 1/sqrt(2) * var.sins * var.sinsh * (Fi[3] - Fi[4]);
114 }
115 else if ( hel[1] == 1 && hel[3] == -1 ){
116 // hel = {1,+,0,-}
117 amp = 1/sqrt(2) * var.sins * var.sinsh * (Fi[3] - Fi[4]);
118 }
119 else return 0.0 ;
120
121 // if negative photon helicity, nucleon flip amplitude changes sign
122 if (hel[0] == -2 && hel[1] != hel[3] ) amp = -1.*amp;
123
124 // there is a factor 8*PI*W between Hi and helicity amplitudes
125 return amp * 8. * M_PI3.14159265358979323846 * sqrt(var.s);
126}
127
128// *********************************************************************************
129void CGLN_Ai(complex<double> s, complex<double> t, complex<double> CGLNA[]){
130 /*
131 * Model for gamma p --> pi0 p from arXiv:1505.02321
132 * Vincent Mathieu June 2015
133 */
134
135 // Model parameters
136 double alpV[3] = { 0.442457, 1.099910, 0.0}; // vector trajectory
137 double alpC[3] = { 0.461811, 0.166543, 0.0}; // vector cut trajectory
138 double alpA[3] = {-0.193332, 1.021300, 0.0}; // axial-vector trajectory
139 double g1 = 3.8873, g4 = -10.1403, g1c = -1.76187, g4c = -3.58089, g2 = -8.887; // residues
140
141 complex<double> Rv, Rc, Ra;
142 complex<double> avec, acut, aaxi;
143 // trajectories:
144 avec = alpV[0] + t*alpV[1] + t*t*alpV[2];
145 acut = alpC[0] + t*alpC[1] + t*t*alpC[2];
146 aaxi = alpA[0] + t*alpA[1] + t*t*alpA[2];
147
148 // Regge factors:
149 complex<double> I(0,1);
150 Rv = cgamma( 1.0 - avec, 0)/2. * ( 1.-exp(-1.*I*M_PI3.14159265358979323846*avec) ) * pow(s,avec-1.);
151 Rc = cgamma( 1.0 - acut, 0)/2. * ( 1.-exp(-1.*I*M_PI3.14159265358979323846*acut) ) * pow(s,acut-1.);
152 Ra = cgamma( 1.0 - aaxi, 0)/2. * ( 1.-exp(-1.*I*M_PI3.14159265358979323846*aaxi) ) * pow(s,aaxi-1.);
153 Rc = Rc / log(s);
154 // IF ONLY VECTOR POLE
155 //Rc = 0; Ra = 0;
156
157 // Scalar amplitudes:
158 CGLNA[1] = t* ( g1*Rv + g1c*Rc);
159 CGLNA[2] = g2 * Ra - ( g1*Rv + g1c*Rc);
160 CGLNA[3] = 0;
161 CGLNA[4] = g4*Rv + g4c*Rc;
162
163 return ;
164
165}
166
167// *********************************************************************************
168void CGLNA2F(struct Kin var, complex<double> Ai[], complex<double> Fi[]){
169 /*
170 * Compute CGLN Fi(s,t) from CGLN Ai(s,t)
171 */
172 complex<double> w;
173 complex<double> E2p, E2m, E4p, E4m;
174
175 w = sqrt(var.s);
176 E2p = var.E2s + var.m2 ;
177 E2m = var.E2s - var.m2 ;
178 E4p = var.E4s + var.m4 ;
179 E4m = var.E4s - var.m4 ;
180
181 Fi[1] = (w - var.m2) * Ai[1] + (var.m3*var.m3 - var.t)/2. * (Ai[3] - Ai[4]);
182 Fi[1] = Fi[1] + ( w - var.m2 ) * ( w - var.m4 ) * Ai[4];
183 Fi[1] = Fi[1] * sqrt( E2p * E4p ) / ( 8. * M_PI3.14159265358979323846 * w);
184
185 Fi[2] = -1.*(w + var.m2) * Ai[1] + (var.m3*var.m3 - var.t)/2. * (Ai[3] - Ai[4]);
186 Fi[2] = Fi[2] + ( w + var.m2 ) * ( w + var.m4 ) * Ai[4];
187 Fi[2] = Fi[2] * sqrt( E2m * E4m ) / ( 8. * M_PI3.14159265358979323846 * w);
188
189 Fi[3] = (w + var.m2) * ( (w - var.m2)*Ai[2] + Ai[3] - Ai[4] );
190 Fi[3] = Fi[3] * sqrt( E2m * E4p) * var.qs / ( 8. * M_PI3.14159265358979323846 * w);
191
192 Fi[4] = (w - var.m2) * ( -1.*(w + var.m2)*Ai[2] + Ai[3] - Ai[4] );
193 Fi[4] = Fi[4] * sqrt( E2p * E4m) * var.qs / ( 8. * M_PI3.14159265358979323846 * w);
194
195 return;
196}
197
198// *********************************************************************************
199void kin2to2(double Ecm, double theta, double mass[], double pa[],double pb[],double pc[], double pd[]){
200/*
201 * Kinematics of the 2-to-2 reaction, a + b --> c + d
202 * Inputs:
203 * Ecm center of mass energy
204 * cos cosine of the scattering angle in the center of mass
205 * mass = {ma,mb,mc,md} vector with the masses of external particles
206 * Outputs:
207 * pa, pb, pc, pd are the momenta of the particles
208 */
209 double Ea, Eb, Ec, Ed; // Energies of the particles
210 double pi, pf; // initial and final breakup momenta
211 double ma, mb, mc, md; // masses of the particles
212
213 ma = mass[1]; mb = mass[2]; mc = mass[3]; md = mass[4];
214
215 // Check that inputs are valid
216 if( ( ma<0 ) || ( mb<0 ) || ( mc<0 ) || ( md<0 ) ) {
217 printf("\n*** Wrong masses in kin2to2 ! *** \n\n");
218 return ;
219 }
220 if( ( Ecm < ma + mb ) || ( Ecm < mc + md ) ) {
221 printf("\n*** Wrong total energy in kin2to2 ! *** \n\n");
222 return ;
223 }
224
225 Ea = ( Ecm*Ecm + ma*ma - mb*mb )/2./Ecm;
226 Eb = ( Ecm*Ecm - ma*ma + mb*mb )/2./Ecm;
227 Ec = ( Ecm*Ecm + mc*mc - md*md )/2./Ecm;
228 Ed = ( Ecm*Ecm - mc*mc + md*md )/2./Ecm;
229
230 pi = sqrt(Ea*Ea - ma*ma);
231 pa[0] = Ea; pa[1] = 0; pa[2] = 0; pa[3] = +pi;
232 pb[0] = Eb; pb[1] = 0; pb[2] = 0; pb[3] = -pi;
233
234 pf = sqrt(Ec*Ec - mc*mc);
235 pc[0] = Ec; pc[1] = +pf*sin(theta); pc[2] = 0; pc[3] = +pf*cos(theta);
236 pd[0] = Ed; pd[1] = -pf*sin(theta); pd[2] = 0; pd[3] = -pf*cos(theta);
237
238 return;
239}
240
241// *********************************************************************************
242void kinematics(complex<double> s, complex<double> t, double mass[], struct Kin *var)
243{
244 double m12, m22, m32, m42; // masses squared
245 complex<double> t0, t1, u ;
246
247 var->s = s;
248 var->t = t;
249
250 var->m1 = mass[1];
251 var->m2 = mass[2];
252 var->m3 = mass[3];
253 var->m4 = mass[4];
254
255 m12 = mass[1] * mass[1] ;
256 m22 = mass[2] * mass[2] ;
257 m32 = mass[3] * mass[3] ;
258 m42 = mass[4] * mass[4] ;
259
260 var->ks = sqrt( lambda(s, m12, m22) / 4. / s );
261 var->qs = sqrt( lambda(s, m32, m42) / 4. / s );
262 var->kt = sqrt( lambda(t, m12, m32) / 4. / t );
263 var->pt = sqrt( lambda(t, m22, m42) / 4. / t );
264
265 // nuclei energies in s- and t-channel frames
266 var->E2s = ( s + m22 - m12 ) / 2. / sqrt(s);
267 var->E4s = ( s + m42 - m32 ) / 2. / sqrt(s);
268 var->E2t = ( t + m22 - m42 ) / 2. / sqrt(t);
269 var->E4t = ( t + m42 - m22 ) / 2. / sqrt(t);
270
271 t1 = pow(m12 - m32 - m22 + m42,2)/(4.*s) - (var->ks - var->qs) * (var->ks - var->qs);
272 t0 = t1 - 4. * var->ks *var->qs ;
273 u = -1.*s - t + m12 + m22 + m32 + m42 ; // Mandelstam s variable
274 var->phi = s * (t - t1) * (t0 - t) ; // Kibble function
275
276 var->coss = 1. + (t - t1)/(2. * var->qs * var->ks) ;
277 var->cost = (t*(s-u) + (m12-m32) * (m22-m42) )/(sqrt(lambda(var->t,m12,m32) * lambda(var->t,m22,m42) ) ) ;
278 var->sins = sqrt( var->phi / s ) / ( 2. * var->qs * var->ks) ;
279 var->sint = sqrt( var->phi / t ) / ( 2. * var->kt * var->pt) ;
280 var->cossh = sqrt( (1. + var->coss ) / 2. );
281 var->sinsh = sqrt( (1. - var->coss ) / 2. );
282 var->costh = sqrt( (1. + var->cost ) / 2. );
283 var->sinth = sqrt( (1. - var->cost ) / 2. );
284
285 return ;
286}
287
288// *********************************************************************************
289complex<double> lambda(complex<double> a, double b, double c){
290// triangle function
291 return a*a + b*b + c*c - 2.*(a*b + b*c + c*a);
292}
293
294// *********************************************************************************
295double snorm(double p[]){
296// Norm squared of a quadri-vector ; p[0] is the energy ; p[1-3] are x,y,z components
297 return p[0]*p[0] - ( p[1]*p[1] + p[2]*p[2] + p[3]*p[3] );
298}
299
300// *********************************************************************************
301complex<double> cgamma(complex<double> z,int OPT)
302{
303 complex<double> I(0,1);
304 complex<double> g, infini= 1e308+ 0.*I; // z0,z1
305 double x0,q1,q2,x,y,th,th1,th2,g0,gr,gi,gr1,gi1;
306 double na,t,x1,y1,sr,si;
307 int j,k;
308 x1=9e9;
309 na=9e9;
310
311 static double a[] = {
312 8.333333333333333e-02,
313 -2.777777777777778e-03,
314 7.936507936507937e-04,
315 -5.952380952380952e-04,
316 8.417508417508418e-04,
317 -1.917526917526918e-03,
318 6.410256410256410e-03,
319 -2.955065359477124e-02,
320 1.796443723688307e-01,
321 -1.39243221690590};
322
323 x = real(z);
324 y = imag(z);
325 if (x > 171) return infini;
326 if ((y == 0.0) && (x == (int)x) && (x <= 0.0))
327 return infini;
328 else if (x < 0.0) {
329 x1 = x;
330 y1 = y;
331 x = -x;
332 y = -y;
333 }
334 x0 = x;
335 if (x <= 7.0) {
336 na = (int)(7.0-x);
337 x0 = x+na;
338 }
339 q1 = sqrt(x0*x0+y*y);
340 th = atan(y/x0);
341 gr = (x0-0.5)*log(q1)-th*y-x0+0.5*log(2.0*M_PI3.14159265358979323846);
342 gi = th*(x0-0.5)+y*log(q1)-y;
343 for (k=0;k<10;k++){
344 t = pow(q1,-1.0-2.0*k);
345 gr += (a[k]*t*cos((2.0*k+1.0)*th));
346 gi -= (a[k]*t*sin((2.0*k+1.0)*th));
347 }
348 if (x <= 7.0) {
349 gr1 = 0.0;
350 gi1 = 0.0;
351 for (j=0;j<na;j++) {
352 gr1 += (0.5*log((x+j)*(x+j)+y*y));
353 gi1 += atan(y/(x+j));
354 }
355 gr -= gr1;
356 gi -= gi1;
357 }
358 if (x1 <= 0.0) {
359 q1 = sqrt(x*x+y*y);
360 th1 = atan(y/x);
361 sr = -sin(M_PI3.14159265358979323846*x)*cosh(M_PI3.14159265358979323846*y);
362 si = -cos(M_PI3.14159265358979323846*x)*sinh(M_PI3.14159265358979323846*y);
363 q2 = sqrt(sr*sr+si*si);
364 th2 = atan(si/sr);
365 if (sr < 0.0) th2 += M_PI3.14159265358979323846;
366 gr = log(M_PI3.14159265358979323846/(q1*q2))-gr;
367 gi = -th1-th2-gi;
368 x = x1;
Value stored to 'x' is never read
369 y = y1;
370 }
371 if (OPT == 0) {
372 g0 = exp(gr);
373 gr = g0*cos(gi);
374 gi = g0*sin(gi);
375 }
376 g = gr + I*gi;
377 return g;
378}