1 |
|
2 |
|
3 |
|
4 |
|
5 |
|
6 |
|
7 |
|
8 |
|
9 |
|
10 |
|
11 |
|
12 |
|
13 |
|
14 | #include "Pi0ReggeModel.h"
|
15 |
|
16 |
|
17 | double Pi0PhotCS_S(double E, double theta, double &BeamSigma){
|
18 |
|
19 |
|
20 |
|
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};
|
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);
|
33 |
|
34 | hel[1] = -1;
|
35 | res[1] = Pi0PhotAmpS(pa, pb, pc, hel);
|
36 |
|
37 | hel[3] = -1;
|
38 | res[2] = Pi0PhotAmpS(pa, pb, pc, hel);
|
39 |
|
40 | hel[1] = +1;
|
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 |
|
46 |
|
47 |
|
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 |
|
56 | complex<double> Pi0PhotAmpS(double pa[],double pb[],double pc[], int hel[]){
|
57 |
|
58 |
|
59 |
|
60 |
|
61 |
|
62 |
|
63 |
|
64 |
|
65 |
|
66 |
|
67 | complex<double> amp = 0.0;
|
68 | complex<double> Ai[5] = {0}, Fi[5] = {0};
|
69 |
|
70 | double pab[4], pca[4], pd[4];
|
71 | struct Kin var;
|
72 | double mass[5]={0};
|
73 |
|
74 |
|
75 | if ( abs(hel[0]) != 2 ) return 0.0;
|
76 |
|
77 | int i;
|
78 | for(i=0;i<4;i++){
|
79 | pab[i] = pa[i] + pb[i];
|
80 |
|
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 |
|
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);
|
96 | CGLN_Ai(var.s,var.t,Ai);
|
97 |
|
98 | CGLNA2F(var, Ai, Fi);
|
99 |
|
100 |
|
101 | if ( hel[1] == 1 && hel[3] == 1 ) {
|
102 |
|
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 |
|
108 | amp = -1/sqrt(2) * var.sins * var.cossh * (Fi[3] + Fi[4]);
|
109 | }
|
110 | else if ( hel[1] == -1 && hel[3] == 1 ){
|
111 |
|
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 |
|
117 | amp = 1/sqrt(2) * var.sins * var.sinsh * (Fi[3] - Fi[4]);
|
118 | }
|
119 | else return 0.0 ;
|
120 |
|
121 |
|
122 | if (hel[0] == -2 && hel[1] != hel[3] ) amp = -1.*amp;
|
123 |
|
124 |
|
125 | return amp * 8. * M_PI3.14159265358979323846 * sqrt(var.s);
|
126 | }
|
127 |
|
128 |
|
129 | void CGLN_Ai(complex<double> s, complex<double> t, complex<double> CGLNA[]){
|
130 |
|
131 |
|
132 |
|
133 |
|
134 |
|
135 |
|
136 | double alpV[3] = { 0.442457, 1.099910, 0.0};
|
137 | double alpC[3] = { 0.461811, 0.166543, 0.0};
|
138 | double alpA[3] = {-0.193332, 1.021300, 0.0};
|
139 | double g1 = 3.8873, g4 = -10.1403, g1c = -1.76187, g4c = -3.58089, g2 = -8.887;
|
140 |
|
141 | complex<double> Rv, Rc, Ra;
|
142 | complex<double> avec, acut, aaxi;
|
143 |
|
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 |
|
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 |
|
155 |
|
156 |
|
157 |
|
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 |
|
168 | void CGLNA2F(struct Kin var, complex<double> Ai[], complex<double> Fi[]){
|
169 |
|
170 |
|
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 |
|
199 | void kin2to2(double Ecm, double theta, double mass[], double pa[],double pb[],double pc[], double pd[]){
|
200 |
|
201 |
|
202 |
|
203 |
|
204 |
|
205 |
|
206 |
|
207 |
|
208 |
|
209 | double Ea, Eb, Ec, Ed;
|
210 | double pi, pf;
|
211 | double ma, mb, mc, md;
|
212 |
|
213 | ma = mass[1]; mb = mass[2]; mc = mass[3]; md = mass[4];
|
214 |
|
215 |
|
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 |
|
242 | void kinematics(complex<double> s, complex<double> t, double mass[], struct Kin *var)
|
243 | {
|
244 | double m12, m22, m32, m42;
|
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 |
|
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 ;
|
274 | var->phi = s * (t - t1) * (t0 - t) ;
|
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 |
|
289 | complex<double> lambda(complex<double> a, double b, double c){
|
290 |
|
291 | return a*a + b*b + c*c - 2.*(a*b + b*c + c*a);
|
292 | }
|
293 |
|
294 |
|
295 | double snorm(double p[]){
|
296 |
|
297 | return p[0]*p[0] - ( p[1]*p[1] + p[2]*p[2] + p[3]*p[3] );
|
298 | }
|
299 |
|
300 |
|
301 | complex<double> cgamma(complex<double> z,int OPT)
|
302 | {
|
303 | complex<double> I(0,1);
|
304 | complex<double> g, infini= 1e308+ 0.*I;
|
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;
|
369 | y = y1;
|
| Value stored to 'y' is never read |
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 | }
|