Line data Source code
1 : #include "forW-MEc.h"
2 : #include "Photos.h"
3 : #include "f_Init.h"
4 : #include "PH_HEPEVT_Interface.h"
5 : #include <cstdlib>
6 : #include<iostream>
7 : using std::cout;
8 : using std::endl;
9 :
10 : namespace Photospp
11 : {
12 :
13 : // COMMON /Kleiss_Stirling/spV,bet
14 : double PhotosMEforW::spV[4],PhotosMEforW::bet[4];
15 :
16 : // COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN
17 : double PhotosMEforW::pi,PhotosMEforW::sw,PhotosMEforW::cw,PhotosMEforW::alphaI,PhotosMEforW::qb,PhotosMEforW::mb,PhotosMEforW::mf1,PhotosMEforW::mf2,PhotosMEforW::qf1,PhotosMEforW::qf2,PhotosMEforW::vf,PhotosMEforW::af,PhotosMEforW::mcLUN;
18 :
19 : //////////////////////////////////////////////////////////////////
20 : // small s_{+,-}(p1,p2) for massless case: //
21 : // p1^2 = p2^2 = 0 //
22 : // //
23 : // k0(0) = 1.d0 //
24 : // k0(1) = 1.d0 //
25 : // k0(2) = 0.d0 Kleisse_Stirling k0 points to X-axis //
26 : // k0(3) = 0.d0 //
27 : // //
28 : //////////////////////////////////////////////////////////////////
29 : complex<double> PhotosMEforW::InProd_zero(double p1[4],int l1,double p2[4],int l2){
30 :
31 :
32 0 : double forSqrt1,forSqrt2,sqrt1,sqrt2;
33 0 : complex<double> Dcmplx;
34 : static complex<double> i_= complex<double>(0.0,1.0);
35 : bool equal;
36 :
37 :
38 :
39 : equal = true;
40 0 : for (int i = 0; i < 4; i++){
41 :
42 0 : if (p1[i]!=p2[i]) equal = equal && false ;
43 : }
44 :
45 :
46 0 : if ( (l1==l2) || equal ) return complex<double>(0.0,0.0);
47 :
48 :
49 0 : else if ( (l1==+1) && (l2==-1) ){
50 :
51 0 : forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
52 0 : forSqrt2 = 1.0/forSqrt1;
53 0 : sqrt1 = sqrt(forSqrt2);
54 0 : sqrt2 = sqrt(forSqrt1);
55 :
56 0 : return (p1[2]+i_*p1[3])*sqrt1 -
57 0 : (p2[2]+i_*p2[3])*sqrt2 ;
58 : }
59 0 : else if ( (l1==-1) && (l2==+1) ){
60 :
61 0 : forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
62 0 : forSqrt2 = 1.0/forSqrt1;
63 0 : sqrt1 = sqrt(forSqrt2);
64 0 : sqrt2 = sqrt(forSqrt1);
65 :
66 0 : return (p2[2]-i_*p2[3])*sqrt2 -
67 0 : (p1[2]-i_*p1[3])*sqrt1 ;
68 : }
69 : else{
70 :
71 :
72 0 : cout << " "<<endl;
73 0 : cout << " ERROR IN InProd_zero:"<<endl;
74 0 : cout << " WRONG VALUES FOR l1,l2: l1,l2 = -1,+1 "<<endl;
75 0 : cout << " " <<endl;
76 0 : exit(-1);
77 : }
78 0 : }
79 :
80 : double PhotosMEforW::InSqrt(double p[4],double q[4]){
81 :
82 0 : return sqrt( (p[0]-p[1]) / (q[0]-q[1]) );
83 : }
84 :
85 : //////////////////////////////////////////////////////////////////
86 : // //
87 : // Inner product for massive spinors: Ub(p1,m1,l1)*U(p2,m2,l2) //
88 : // //
89 : //////////////////////////////////////////////////////////////////
90 :
91 : complex<double> PhotosMEforW::InProd_mass(double p1[4],double m1,int l1,double p2[4],double m2,int l2){
92 : double sqrt1,sqrt2,forSqrt1;
93 :
94 :
95 0 : if ((l1==+1)&&(l2==+1)) {
96 0 : forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
97 0 : sqrt1 = sqrt(forSqrt1);
98 0 : sqrt2 = 1.0/sqrt1;
99 0 : return complex<double>(m1*sqrt2+m2*sqrt1,0.0);
100 : }
101 0 : else if ((l1==+1)&&(l2==-1))
102 0 : return InProd_zero(p1,+1,p2,-1);
103 :
104 0 : else if ((l1==-1)&&(l2==+1))
105 0 : return InProd_zero(p1,-1,p2,+1);
106 :
107 0 : else if ((l1==-1)&&(l2==-1)){
108 0 : forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
109 0 : sqrt1 = sqrt(forSqrt1);
110 0 : sqrt2 = 1.0/sqrt1;
111 0 : return complex<double>(m1*sqrt2+m2*sqrt1,0.0);
112 : }
113 : else {
114 0 : cout <<" " <<endl;
115 0 : cout <<" ERROR IN InProd_mass.."<<endl;
116 0 : cout <<" WRONG VALUES FOR l1,l2"<<endl;
117 0 : cout <<" " <<endl;
118 0 : exit(-1);
119 : }
120 0 : }
121 :
122 : /////////////////////////////////////////////////////////////////////
123 : // //
124 : // this is small B_{s}(k,p) function when TrMartix is diaagonal!! //
125 : // //
126 : /////////////////////////////////////////////////////////////////////
127 : complex<double> PhotosMEforW::BsFactor(int s,double k[4],double p[4],double m){
128 0 : double forSqrt1,sqrt1;
129 0 : complex<double> inPr1;
130 :
131 0 : if ( s==1 ){
132 :
133 0 : inPr1 = InProd_zero(k,+1,p,-1);
134 0 : forSqrt1 = (p[0]-p[1])/(k[0]-k[1]);
135 0 : sqrt1 = sqrt(2.0*forSqrt1);
136 : //BsFactor =
137 0 : return inPr1*sqrt1;
138 : }
139 :
140 0 : else if ( s==-1 ){
141 :
142 0 : inPr1 = InProd_zero(k,-1,p,+1);
143 0 : forSqrt1 = (p[0]-p[1])/(k[0]-k[1]);
144 0 : sqrt1 = sqrt(2.0*forSqrt1);
145 : //BsFactor =
146 0 : return inPr1*sqrt1;
147 : }
148 : else{
149 :
150 0 : cout << " "<<endl;
151 0 : cout << " ERROR IN BsFactor: "<<endl;
152 0 : cout << " WRONG VALUES FOR s : s = -1,+1"<<endl;
153 0 : cout << " " <<endl;
154 0 : exit(-1);
155 : }
156 0 : }
157 :
158 :
159 :
160 :
161 : //======================================================================
162 : //
163 : // Eikonal factor of decay W->l_1+l_2+\gamma in terms of K&S objects !
164 : //
165 : // EikFactor = q1*eps.p1/k.p1 + q2*eps.p2/k.p2 - q3*eps.p3/k.p3
166 : //
167 : // indices 1,2 are for charged decay products
168 : // index 3 is for W
169 : //
170 : // q - charge
171 : //
172 : //======================================================================
173 : complex<double> PhotosMEforW::WDecayEikonalKS_1ph(double p3[4],double p1[4],double p2[4],double k[4],int s){
174 :
175 : double scalProd1,scalProd2,scalProd3;
176 0 : complex<double> wdecayeikonalks_1ph,BSoft1,BSoft2;
177 :
178 0 : scalProd1 = p1[0]*k[0]-p1[1]*k[1]-p1[2]*k[2]-p1[3]*k[3];
179 0 : scalProd2 = p2[0]*k[0]-p2[1]*k[1]-p2[2]*k[2]-p2[3]*k[3];
180 0 : scalProd3 = p3[0]*k[0]-p3[1]*k[1]-p3[2]*k[2]-p3[3]*k[3];
181 :
182 :
183 0 : BSoft1 = BsFactor(s,k,p1,mf1);
184 0 : BSoft2 = BsFactor(s,k,p2,mf2);
185 :
186 : //WDecayEikonalKS_1ph =
187 0 : return sqrt(pi/alphaI)*(-(qf1/scalProd1+qb/scalProd3)*BSoft1
188 0 : +(qf2/scalProd2-qb/scalProd3)*BSoft2);
189 :
190 0 : }
191 :
192 : //======================================================================
193 : //
194 : // Gauge invariant soft factor for decay!!
195 : // Gmass2 -- photon mass square
196 : //
197 : //======================================================================
198 : complex<double> PhotosMEforW::SoftFactor(int s,double k[4],double p1[4],double m1,double p2[4],double m2,double Gmass2){
199 :
200 : double ScalProd1,ScalProd2;
201 0 : complex<double> BsFactor2,BsFactor1;
202 :
203 :
204 0 : ScalProd1 = k[0]*p1[0]-k[1]*p1[1]-k[2]*p1[2]-k[3]*p1[3];
205 0 : ScalProd2 = k[0]*p2[0]-k[1]*p2[1]-k[2]*p2[2]-k[3]*p2[3];
206 :
207 0 : BsFactor1 = BsFactor(s,k,p1,m1);
208 0 : BsFactor2 = BsFactor(s,k,p2,m2);
209 :
210 0 : return + BsFactor2/2.0/(ScalProd2-Gmass2)
211 0 : - BsFactor1/2.0/(ScalProd1-Gmass2);
212 0 : }
213 :
214 : //#############################################################################
215 : // #
216 : // \ eps(k,0,s) #
217 : // / #
218 : // _\ #
219 : // /\ #
220 : // \ #
221 : // / #
222 : // ---<----------\-------------<--- #
223 : // Ub(p1,m1,l1) U(p2,m2,l2) #
224 : // #
225 : // #
226 : // definition of arbitrary light-like vector beta!! #
227 : // #
228 : // bet[0] = 1.d0 #
229 : // bet[1] = 1.d0 #
230 : // bet[2] = 0.d0 <==> bet == k0 expression becomes easy!! #
231 : // bet[3] = 0.d0 #
232 : //#############################################################################
233 :
234 : complex<double> PhotosMEforW::TrMatrix_zero(double p1[4],double m1,int l1,double k[4],int s,double p2[4],double m2,int l2){
235 :
236 : double forSqrt1,forSqrt2;
237 : // double p1_1[4],p2_1[4];
238 0 : double sqrt1,sqrt2; // ,scalProd1,scalProd2;
239 0 : complex<double> inPr1,inPr2,inPr3;
240 : bool equal;
241 :
242 : equal = true;
243 0 : for (int i = 0; i < 4; i++)
244 0 : if (p1[i] != p2[i]) equal = equal&&false;
245 :
246 :
247 :
248 0 : if ( (m1==m2)&&(equal) ){
249 : //..
250 : //.. when: p1=p2=p <=> m1=m2 TrMatrix_zero is diagonal
251 : //..
252 0 : if ( (l1==+1)&&(l2==+1) ){
253 :
254 0 : inPr1 = InProd_zero(k,+s,p1,-s);
255 0 : forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
256 0 : sqrt1 = sqrt(2.0*forSqrt1);
257 :
258 0 : return sqrt1*inPr1;
259 : }
260 :
261 0 : else if ( (l1==+1)&&(l2==-1) ){
262 :
263 0 : return complex<double>(0.0,0.0);}
264 :
265 :
266 0 : else if ( (l1==-1)&&(l2==+1) ){
267 :
268 0 : return complex<double>(0.0,0.0);
269 : }
270 :
271 0 : else if ( (l1==-1)&&(l2==-1) ){
272 :
273 0 : inPr1 = InProd_zero(k,+s,p1,-s);
274 0 : forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
275 0 : sqrt1 = sqrt(2.0*forSqrt1);
276 :
277 0 : return sqrt1*inPr1;
278 : }
279 :
280 : else{
281 :
282 0 : cout << "" <<endl;
283 0 : cout << " ERROR IN TrMatrix_zero: " <<endl;
284 0 : cout << " WRONG VALUES FOR l1,l2,s" <<endl;
285 0 : cout << "" <<endl;
286 0 : exit(-1);
287 :
288 : }
289 :
290 : }
291 :
292 0 : if ( (l1==+1)&&(l2==+1)&&(s==+1) ){
293 :
294 0 : inPr1 = InProd_zero(k,+1,p1,-1);
295 0 : forSqrt1 = (p2[0]-p2[1])/(k[0]-k[1]);
296 0 : sqrt1 = sqrt(2.0*forSqrt1);
297 :
298 0 : return sqrt1*inPr1;
299 : }
300 0 : else if ( (l1==+1)&&(l2==-1)&&(s==+1) ) {
301 :
302 0 : return complex<double>(0.0,0.0);
303 : }
304 :
305 0 : else if( (l1==-1)&&(l2==+1)&&(s==+1) ){
306 :
307 0 : forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
308 0 : forSqrt2 = 1.0/forSqrt1;
309 0 : sqrt1 = sqrt(2.0*forSqrt1);
310 0 : sqrt2 = sqrt(2.0*forSqrt2);
311 :
312 0 : return complex<double>(m2*sqrt1-m1*sqrt2,0.0);
313 : }
314 0 : else if ( (l1==-1)&&(l2==-1)&&(s==+1) ){
315 :
316 0 : inPr1 = InProd_zero(k,+1,p2,-1);
317 0 : forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
318 0 : sqrt1 = sqrt(2.0*forSqrt1);
319 :
320 0 : return inPr1*sqrt1;
321 : }
322 :
323 0 : else if ( (l1==+1)&&(l2==+1)&&(s==-1) ){
324 :
325 0 : inPr1 = -InProd_zero(k,-1,p2,+1);
326 0 : forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
327 0 : sqrt1 = sqrt(2.0*forSqrt1);
328 :
329 0 : return -sqrt1*inPr1;
330 : }
331 :
332 0 : else if ( (l1==+1)&&(l2==-1)&&(s==-1) ){
333 :
334 0 : forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
335 0 : forSqrt2 = 1.0/forSqrt1;
336 0 : sqrt1 = sqrt(2.0*forSqrt1);
337 0 : sqrt2 = sqrt(2.0*forSqrt2);
338 :
339 0 : return complex<double>(m2*sqrt1-m1*sqrt2,0.0);
340 : }
341 :
342 0 : else if ( (l1==-1)&&(l2==+1)&&(s==-1) ){
343 :
344 0 : return complex<double>(0.0,0.0);
345 : }
346 :
347 0 : else if( (l1==-1)&&(l2==-1)&&(s==-1) ){
348 :
349 0 : inPr1 = -InProd_zero(k,-1,p1,+1);
350 0 : forSqrt1 = (p2[0]-p2[1])/(k[0]-k[1]);
351 0 : sqrt1 = sqrt(2.0*forSqrt1);
352 :
353 0 : return -inPr1*sqrt1;
354 : }
355 : else {
356 :
357 0 : cout << "" << endl;
358 0 : cout << " ERROR IN TrMatrix_zero: " << endl;
359 0 : cout << " WRONG VALUES FOR l1,l2,s" << endl;
360 0 : cout << "" << endl;
361 0 : exit(-1);
362 : }
363 :
364 0 : }
365 :
366 :
367 :
368 : ////////////////////////////////////////////////////////////////
369 : // transition matrix for massive boson //
370 : // //
371 : // //
372 : // \ eps(k,m,s) //
373 : // / //
374 : // _\ //
375 : // /\ k //
376 : // \ //
377 : // <-- p1 / <-- p2 //
378 : // ---<----------\----------<--- //
379 : // Ub(p1,m1,l1) U(p2,m2,l2) //
380 : // //
381 : ////////////////////////////////////////////////////////////////
382 : complex<double> PhotosMEforW::TrMatrix_mass(double p1[4],double m1,int l1,double k[4],double m,int s,double p2[4],double m2,int l2){
383 :
384 :
385 : double forSqrt1,forSqrt2;
386 0 : double k_1[4],k_2[4];
387 0 : double forSqrt3,forSqrt4,sqrt3,sqrt1,sqrt2,sqrt4;
388 0 : complex<double> inPr1,inPr2,inPr3,inPr4;
389 :
390 0 : for (int i = 0; i < 4; i++) {
391 0 : k_1[i] = 1.0/2.0*(k[i] - m*spV[i]);
392 0 : k_2[i] = 1.0/2.0*(k[i] + m*spV[i]);
393 : }
394 :
395 0 : if ( (l1==+1)&&(l2==+1)&&(s==0) ){
396 :
397 0 : inPr1 = InProd_zero(p1,+1,k_2,-1);
398 0 : inPr2 = InProd_zero(p2,-1,k_2,+1);
399 0 : inPr3 = InProd_zero(p1,+1,k_1,-1);
400 0 : inPr4 = InProd_zero(p2,-1,k_1,+1);
401 0 : sqrt1 = sqrt(p1[0]-p1[1]);
402 0 : sqrt2 = sqrt(p2[0]-p2[1]);
403 0 : sqrt3 = m1*m2/sqrt1/sqrt2;
404 :
405 0 : return
406 0 : (inPr1*inPr2-inPr3*inPr4)*(vf+af)/m
407 0 : + (k_1[0]-k_2[0]-k_1[1]+k_2[1])*sqrt3*(vf-af)/m;
408 : }
409 :
410 0 : else if ( (l1==+1)&&(l2==-1)&&(s==0) ){
411 :
412 0 : inPr1 = InProd_zero(p1,+1,k_1,-1);
413 0 : inPr2 = InProd_zero(p1,+1,k_2,-1);
414 0 : inPr3 = InProd_zero(p2,+1,k_2,-1);
415 0 : inPr4 = InProd_zero(p2,+1,k_1,-1);
416 :
417 0 : forSqrt1 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
418 0 : forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
419 0 : forSqrt3 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
420 0 : forSqrt4 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
421 0 : sqrt1 = sqrt(forSqrt1);
422 0 : sqrt2 = sqrt(forSqrt2);
423 0 : sqrt3 = sqrt(forSqrt3);
424 0 : sqrt4 = sqrt(forSqrt4);
425 :
426 0 : return
427 0 : (inPr1*sqrt1 - inPr2*sqrt2)*(vf+af)*m2/m
428 0 : + (inPr3*sqrt3 - inPr4*sqrt4)*(vf-af)*m1/m;
429 : }
430 0 : else if ( (l1==-1)&&(l2==+1)&&(s==0) ){
431 :
432 0 : inPr1 = InProd_zero(p1,-1,k_1,+1);
433 0 : inPr2 = InProd_zero(p1,-1,k_2,+1);
434 0 : inPr3 = InProd_zero(p2,-1,k_2,+1);
435 0 : inPr4 = InProd_zero(p2,-1,k_1,+1);
436 :
437 0 : forSqrt1 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
438 0 : forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
439 0 : forSqrt3 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
440 0 : forSqrt4 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
441 0 : sqrt1 = sqrt(forSqrt1);
442 0 : sqrt2 = sqrt(forSqrt2);
443 0 : sqrt3 = sqrt(forSqrt3);
444 0 : sqrt4 = sqrt(forSqrt4);
445 :
446 0 : return
447 0 : (inPr1*sqrt1 - inPr2*sqrt2)*(vf-af)*m2/m
448 0 : + (inPr3*sqrt3 - inPr4*sqrt4)*(vf+af)*m1/m;
449 : }
450 0 : else if ( (l1==-1)&&(l2==-1)&&(s==0) ){
451 :
452 0 : inPr1 = InProd_zero(p2,+1,k_2,-1);
453 0 : inPr2 = InProd_zero(p1,-1,k_2,+1);
454 0 : inPr3 = InProd_zero(p2,+1,k_1,-1);
455 0 : inPr4 = InProd_zero(p1,-1,k_1,+1);
456 0 : sqrt1 = sqrt(p1[0]-p1[1]);
457 0 : sqrt2 = sqrt(p2[0]-p2[1]);
458 0 : sqrt3 = m1*m2/sqrt1/sqrt2;
459 :
460 0 : return
461 0 : (inPr1*inPr2 - inPr3*inPr4)*(vf-af)/m
462 0 : + (k_1[0]-k_2[0]-k_1[1]+k_2[1])*sqrt3*(vf+af)/m;
463 : }
464 0 : else if ( (l1==+1)&&(l2==+1)&&(s==+1) ){
465 :
466 0 : inPr1 = InProd_zero(p1,+1,k_1,-1);
467 0 : inPr2 = InProd_zero(k_2,-1,p2,+1);
468 0 : inPr3 = inPr1*inPr2;
469 :
470 0 : forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
471 0 : forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
472 0 : sqrt1 = sqrt(forSqrt1);
473 0 : sqrt2 = sqrt(forSqrt2);
474 0 : sqrt3 = m1*m2*sqrt1*sqrt2;
475 :
476 0 : return
477 0 : sqrt(2.0)/m*(inPr3*(vf+af)+sqrt3*(vf-af));
478 : }
479 :
480 0 : else if ( (l1==+1)&&(l2==-1)&&(s==+1) ){
481 :
482 0 : inPr1 = InProd_zero(p1,+1,k_1,-1);
483 0 : inPr2 = InProd_zero(p2,+1,k_1,-1);
484 :
485 0 : forSqrt1 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
486 0 : forSqrt2 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
487 0 : sqrt1 = m2*sqrt(forSqrt1);
488 0 : sqrt2 = m1*sqrt(forSqrt2);
489 :
490 0 : return
491 0 : sqrt(2.0)/m*( + inPr1*sqrt1*(vf+af)
492 0 : - inPr2*sqrt2*(vf-af)
493 : );
494 : }
495 0 : else if ( (l1==-1)&&(l2==+1)&&(s==+1) ){
496 :
497 0 : inPr1 = InProd_zero(k_2,-1,p2,+1);
498 0 : inPr2 = InProd_zero(k_2,-1,p1,+1);
499 :
500 0 : forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
501 0 : forSqrt2 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
502 0 : sqrt1 = m1*sqrt(forSqrt1);
503 0 : sqrt2 = m2*sqrt(forSqrt2);
504 :
505 0 : return
506 0 : sqrt(2.0)/m*( + inPr1*sqrt1*(vf+af)
507 0 : - inPr2*sqrt2*(vf-af)
508 : );
509 : }
510 0 : else if ( (l1==-1)&&(l2==-1)&&(s==+1) ){
511 :
512 0 : inPr1 = InProd_zero(p2,+1,k_1,-1);
513 0 : inPr2 = InProd_zero(k_2,-1,p1,+1);
514 0 : inPr3 = inPr1*inPr2;
515 :
516 0 : forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
517 0 : forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
518 0 : sqrt1 = sqrt(forSqrt1);
519 0 : sqrt2 = sqrt(forSqrt2);
520 0 : sqrt3 = m1*m2*sqrt1*sqrt2;
521 :
522 0 : return
523 0 : sqrt(2.0)/m*(inPr3*(vf-af)+sqrt3*(vf+af));
524 : }
525 :
526 0 : else if ( (l1==+1)&&(l2==+1)&&(s==-1) ){
527 :
528 0 : inPr1 = InProd_zero(p2,-1,k_1,+1);
529 0 : inPr2 = InProd_zero(k_2,+1,p1,-1);
530 0 : inPr3 = inPr1*inPr2;
531 :
532 0 : forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
533 0 : forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
534 0 : sqrt1 = sqrt(forSqrt1);
535 0 : sqrt2 = sqrt(forSqrt2);
536 0 : sqrt3 = m1*m2*sqrt1*sqrt2;
537 :
538 0 : return
539 0 : sqrt(2.0)/m*(inPr3*(vf+af)+sqrt3*(vf-af));
540 : }
541 0 : else if ( (l1==+1)&&(l2==-1)&&(s==-1) ){
542 :
543 0 : inPr1 = InProd_zero(k_2,+1,p2,-1);
544 0 : inPr2 = InProd_zero(k_2,+1,p1,-1);
545 :
546 0 : forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
547 0 : forSqrt2 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
548 0 : sqrt1 = m1*sqrt(forSqrt1);
549 0 : sqrt2 = m2*sqrt(forSqrt2);
550 :
551 0 : return
552 0 : sqrt(2.0)/m*(+ inPr1*sqrt1*(vf-af)
553 0 : - inPr2*sqrt2*(vf+af)
554 : );
555 : }
556 0 : else if ( (l1==-1)&&(l2==+1)&&(s==-1) ){
557 :
558 0 : inPr1 = InProd_zero(p1,-1,k_1,+1);
559 0 : inPr2 = InProd_zero(p2,-1,k_1,+1);
560 :
561 0 : forSqrt1 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
562 0 : forSqrt2 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
563 0 : sqrt1 = m2*sqrt(forSqrt1);
564 0 : sqrt2 = m1*sqrt(forSqrt2);
565 :
566 0 : return
567 0 : sqrt(2.0)/m*(+ inPr1*sqrt1*(vf-af)
568 0 : - inPr2*sqrt2*(vf+af)
569 : );
570 : }
571 0 : else if ( (l1==-1)&&(l2==-1)&&(s==-1) ){
572 :
573 0 : inPr1 = InProd_zero(p1,-1,k_1,+1);
574 0 : inPr2 = InProd_zero(k_2,+1,p2,-1);
575 0 : inPr3 = inPr1*inPr2;
576 :
577 0 : forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
578 0 : forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
579 0 : sqrt1 = sqrt(forSqrt1);
580 0 : sqrt2 = sqrt(forSqrt2);
581 0 : sqrt3 = m1*m2*sqrt1*sqrt2;
582 :
583 0 : return
584 0 : sqrt(2.0)/m*(inPr3*(vf-af)+sqrt3*(vf+af));
585 : }
586 :
587 : else{
588 :
589 0 : cout << " "<< endl;
590 0 : cout << " TrMatrix_mass: Wrong values for l1,l2,s:"<< endl;
591 0 : cout << " l1,l2 = -1,+1; s = -1,0,1 "<< endl;
592 0 : cout << " "<< endl;
593 0 : exit(-1);
594 :
595 : }
596 :
597 0 : }
598 :
599 :
600 :
601 : //======================================================================
602 : // =
603 : // p1,mf1,l1 =
604 : // / =
605 : // \/_ =
606 : // / =
607 : // p3,mb,l3 / =
608 : // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
609 : // \ =
610 : // _\/ =
611 : // \ =
612 : // p2,mf2,l2 =
613 : // INPUT : p1,m1,l1; p2,m2,l2; p3,m3,l3 -- momenta,mass and helicity =
614 : // =
615 : // OUTPUT: value of functions -- decay amplitude =
616 : // =
617 : //======================================================================
618 : complex<double> PhotosMEforW::WDecayBornAmpKS_1ph(double p3[4],int l3,double p1[4],int l1,double p2[4],int l2){
619 :
620 0 : double coeff;
621 :
622 :
623 0 : coeff = sqrt(pi/alphaI/2.0)/sw; // vertex: g/2/sqrt(2)
624 :
625 0 : return coeff*TrMatrix_mass(p2,mf2,l2,p3,mb,l3,p1,-mf1,-l1);
626 0 : }
627 :
628 :
629 : //======================================================================
630 : // k,0,l =
631 : // \ p1,mf1,l1 =
632 : // / / =
633 : // \ \/_ =
634 : // / / =
635 : // p3,mb,l3 \ / =
636 : // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
637 : // \ =
638 : // _\/ =
639 : // \ =
640 : // p2,mf2,l2 =
641 : // { + } =
642 : // p1,mf1,l1 =
643 : // / =
644 : // \/_~~~~~~~ k,0,s =
645 : // / =
646 : // p3,mb,l3 / =
647 : // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
648 : // \ =
649 : // _\/ =
650 : // \ =
651 : // p2,mf2,l2 =
652 : // { + } =
653 : // p1,mf1,l1 =
654 : // / =
655 : // \/_ =
656 : // / =
657 : // p3,mb,l3 / =
658 : // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
659 : // \ =
660 : // _\/ ~~~~~~~ k,0,s =
661 : // \ =
662 : // p2,mf2,l2 =
663 : // =
664 : // all momentas, exept k are incoming !!! =
665 : // =
666 : // This function culculates The W-ff\gamma decay amplitude into permion=
667 : // pair and one photon using Kleisse&Stirling method for helicity =
668 : // amplitudes, which includes three above feynman diagramms.. =
669 : // =
670 : // INPUT : p1,m1,l1; p2,m2,l2; p3,m3,l3 -- momenta,mass and helicity =
671 : // =
672 : // OUTPUT: value of functions -- decay amplitude =
673 : // =
674 : //======================================================================
675 :
676 : complex<double> PhotosMEforW::WDecayAmplitudeKS_1ph(double p3[4],int l3,double p1[4],int l1,double p2[4],int l2,double k[4],int s){
677 :
678 0 : double scalProd1,scalProd2,scalProd3,coeff; //,theta3,ph3;
679 0 : complex<double> bornAmp,TrMx1,TrMx2;
680 0 : complex<double> BSoft1,BSoft2;
681 :
682 0 : coeff = sqrt(2.0)*pi/sw/alphaI; // vertex: g/2/sqrt[2] * e
683 :
684 0 : scalProd1 = p1[0]*k[0]-p1[1]*k[1]-p1[2]*k[2]-p1[3]*k[3];
685 0 : scalProd2 = p2[0]*k[0]-p2[1]*k[1]-p2[2]*k[2]-p2[3]*k[3];
686 0 : scalProd3 = p3[0]*k[0]-p3[1]*k[1]-p3[2]*k[2]-p3[3]*k[3];
687 :
688 0 : BSoft1 = BsFactor(s,k,p1,mf1);
689 0 : BSoft2 = BsFactor(s,k,p2,mf2);
690 0 : bornAmp = TrMatrix_mass(p2,mf2,l2,p3,mb,l3,p1,-mf1,-l1);
691 0 : TrMx1 = complex<double>(0.0,0.0);
692 0 : TrMx2 = complex<double>(0.0,0.0);
693 :
694 0 : for (int la1 = -1; la1< 3 ; la1+=2) {
695 : // DO la1=-1,1,2
696 0 : TrMx1 = TrMx1 + TrMatrix_zero(k,0.0,-la1,k,s,p1,-mf1,-l1)*
697 0 : TrMatrix_mass(p2,mf2,l2,p3,mb,l3,k,0.0,-la1);
698 0 : TrMx2 = TrMx2 + TrMatrix_zero(p2,mf2,l2,k,s,k,0.0,la1)*
699 0 : TrMatrix_mass(k,0.0,la1,p3,mb,l3,p1,-mf1,-l1);
700 : }
701 :
702 0 : return coeff * (
703 0 : + (-(qf1/scalProd1+qb/scalProd3)*BSoft1 // IR-divergent part of amplitude
704 0 : +(qf2/scalProd2-qb/scalProd3)*BSoft2)/2.0*bornAmp
705 : //
706 0 : - (qf1/scalProd1+qb/scalProd3)*TrMx1/2.0 // IR-finite part of amplitude
707 0 : + (qf2/scalProd2-qb/scalProd3)*TrMx2/2.0
708 : );
709 0 : }
710 :
711 :
712 :
713 : //========================================================
714 : // The squared eikonal factor for W decay =
715 : // into fermion pair and one photon =
716 : // INPUT : =
717 : // =
718 : // OUTPUT: =
719 : //========================================================
720 :
721 : double PhotosMEforW::WDecayEikonalSqrKS_1ph(double p3[4],double p1[4],double p2[4],double k[4]){
722 : double spinSumAvrg;
723 0 : complex<double> wDecAmp;
724 :
725 : spinSumAvrg = 0.0;
726 0 : for (int s = -1; s< 3 ; s+=2) {
727 0 : wDecAmp = WDecayEikonalKS_1ph(p3,p1,p2,k,s);
728 0 : spinSumAvrg = spinSumAvrg + real(wDecAmp*conj(wDecAmp));
729 : }
730 0 : return spinSumAvrg;
731 0 : }
732 :
733 : //========================================================
734 : // The squared eikonal factor for W decay =
735 : // into fermion pair and one photon =
736 : // INPUT : =
737 : // =
738 : // OUTPUT: =
739 : //========================================================
740 :
741 : double PhotosMEforW::WDecayBornAmpSqrKS_1ph(double p3[4],double p1[4],double p2[4]){
742 : double spinSumAvrg;
743 0 : complex<double> wDecAmp;
744 :
745 : spinSumAvrg = 0.0;
746 0 : for (int l3 = -1; l3< 2 ; l3++) {
747 0 : for (int l1 = -1; l1< 3 ; l1+=2) {
748 0 : for (int l2 = -1; l2< 3 ; l2+=2) {
749 0 : wDecAmp = WDecayBornAmpKS_1ph(p3,l3,p1,l1,p2,l2);
750 0 : spinSumAvrg = spinSumAvrg + real(wDecAmp*conj(wDecAmp));
751 : }
752 : }
753 : }
754 0 : return spinSumAvrg;
755 0 : }
756 :
757 :
758 :
759 : //========================================================
760 : // The squared amplitude for W decay =
761 : // into fermion pair and one photon =
762 : // INPUT : =
763 : // =
764 : // OUTPUT: =
765 : //========================================================
766 :
767 : double PhotosMEforW::WDecayAmplitudeSqrKS_1ph(double p3[4],double p1[4],double p2[4], double k[4]){
768 :
769 : double spinSumAvrg;
770 0 : complex<double> wDecAmp;
771 :
772 : spinSumAvrg = 0.0;
773 0 : for (int l3 = -1; l3< 2 ; l3++) {
774 0 : for (int l1 = -1; l1< 3 ; l1+=2) {
775 0 : for (int l2 = -1; l2< 3 ; l2+=2) {
776 0 : for (int s = -1; s < 3 ; s+=2) {
777 0 : wDecAmp = WDecayAmplitudeKS_1ph(p3,l3,p1,l1,p2,l2,k,s);
778 0 : spinSumAvrg = spinSumAvrg + real(wDecAmp*conj(wDecAmp));
779 : }
780 : }
781 : }
782 : }
783 0 : return spinSumAvrg;
784 :
785 :
786 :
787 : //$$$
788 : //$$$
789 : //$$$
790 : //$$$
791 : //$$$
792 : //$$$
793 : //$$$
794 : //$$$
795 : //$$$
796 : //$$$
797 : //$$$
798 : //$$$
799 : //$$$
800 : //$$$
801 : //$$$
802 : //$$$ WffGammaME.f ends above:
803 : //$$$
804 : //$$$
805 : //$$$
806 : //$$$
807 0 : }
808 :
809 :
810 :
811 : //C========================================================== ==
812 : //C========================================================== ==
813 : //C these will be public for PHOTOS functions of W_ME class ==
814 : //C========================================================== ==
815 : //C========================================================== ==
816 :
817 : double PhotosMEforW::SANC_WT(double PW[4],double PNE[4],double PMU[4],double PPHOT[4],double B_PW[4],double B_PNE[4],double B_PMU[4]){
818 :
819 :
820 : //.. Exact amplitude square
821 0 : double AMPSQR=WDecayAmplitudeSqrKS_1ph(PW,PNE,PMU,PPHOT);
822 :
823 0 : double EIKONALFACTOR=WDecayBornAmpSqrKS_1ph(B_PW,B_PNE,B_PMU)
824 0 : *WDecayEikonalSqrKS_1ph(PW,PNE,PMU,PPHOT);
825 :
826 : //.. New weight
827 :
828 : // cout << 'B_pne=',B_PNE << endl;
829 : // cout << 'B_PMU=',B_PMU << endl;
830 : // cout << 'bornie=',WDecayBornAmpSqrKS_1ph(B_PW,B_PNE,B_PMU) << endl;
831 :
832 : // cout << ' ' << endl;
833 : // cout << ' pne=',pne << endl;
834 : // cout << ' pmu=',pmu << endl;
835 : // cout << 'pphot=',pphot << endl;
836 : // cout << ' ' << endl;
837 : // cout << ' b_pw=',B_PW << endl;
838 : // cout << ' b_pne=',B_PNE << endl;
839 : // cout << 'b_pmu=',B_PMU << endl;
840 :
841 : // cout << 'cori=',AMPSQR/EIKONALFACTOR,AMPSQR,EIKONALFACTOR << endl;
842 :
843 0 : return AMPSQR/EIKONALFACTOR;
844 : //
845 : // return (1-8*EMU*XPH*(1-COSTHG*BETA)*
846 : // (MCHREN+2*XPH*SQRT(MPASQR))/
847 : // MPASQR**2/(1-MCHREN/MPASQR)/(4-MCHREN/MPASQR))
848 : }
849 :
850 :
851 : void PhotosMEforW::SANC_INIT1(double QB0,double QF20,double MF10,double MF20,double MB0){
852 0 : qb =QB0;
853 0 : qf2=QF20;
854 0 : mf1=MF10;
855 0 : mf2=MF20;
856 0 : mb =MB0;
857 0 : }
858 :
859 : void PhotosMEforW::SANC_INIT(double ALPHA,int PHLUN){
860 :
861 :
862 : static int SANC_MC_INIT=-123456789;
863 :
864 : //... Initialization of the W->l\nu\gamma
865 : //... decay Matrix Element parameters
866 0 : if (SANC_MC_INIT==-123456789){
867 0 : SANC_MC_INIT=1;
868 :
869 0 : pi=4*atan(1.0);
870 0 : qf1=0.0; // neutrino charge
871 0 : mf1=1.0e-10; // newutrino mass
872 0 : vf=1.0; // V&A couplings
873 0 : af=1.0;
874 0 : alphaI=1.0/ALPHA;
875 0 : cw=0.881731727; // Weak Weinberg angle
876 0 : sw=0.471751166;
877 :
878 :
879 : //... An auxilary K&S vectors
880 0 : bet[0]= 1.0;
881 0 : bet[1]= 0.0722794881816159;
882 0 : bet[2]=-0.994200045099866;
883 0 : bet[3]= 0.0796363353729248;
884 :
885 0 : spV[0]= 0.0;
886 0 : spV[1]= 7.22794881816159e-2;
887 0 : spV[2]=-0.994200045099866;
888 0 : spV[3]= 7.96363353729248e-2;
889 :
890 0 : mcLUN = PHLUN;
891 0 : }
892 0 : }
893 : //----------------------------------------------------------------------
894 : //
895 : // PHOTOS: PHOtos Boson W correction weight
896 : //
897 : // Purpose: calculates correction weight due to amplitudes of
898 : // emission from W boson. It is ecact, but not verified
899 : // for exponentiation yet.
900 : //
901 : //
902 : //
903 : //
904 : // Input Parameters: Common /PHOEVT/, with photon added.
905 : // wt to be corrected
906 : //
907 : //
908 : //
909 : // Output Parameters: wt
910 : //
911 : // Author(s): G. Nanava, Z. Was Created at: 13/03/03
912 : // Last Update: 22/06/13
913 : //
914 : //----------------------------------------------------------------------
915 : void PhotosMEforW::PHOBWnlo(double *WT){
916 : // FILE *PHLUN = stdout; // printouts from matrix element calculations
917 : // directed with phlun still
918 : int phlun=6;
919 : double EMU,MCHREN,BETA,COSTHG,MPASQR,XPH;
920 0 : double PW[4],PMU[4],PPHOT[4],PNE[4];
921 0 : double B_PW[4],B_PNE[4],B_PMU[4]; //,AMPSQR;
922 : static int i=1;
923 : int I,IJ,I3,I4,JJ;
924 : double MB,MF1,MF2,QB,QF2;
925 : // double pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af;
926 :
927 :
928 : //! write(*,*) 'IDPHOs=',IDPHO(1),IDPHO(2),IDPHO(3),IDPHO(4),IDPHO(5)
929 : //! write(*,*) 'IDPHOs=',pho.jdahep[1-i][1-i],npho
930 : //! write(*,*) 'hep.IDPHOs=',hep.IDhep(1),hep.IDhep(2),hep.IDhep(3),hep.IDhep(4),hep.IDhep(5)
931 :
932 : //--
933 0 : if(abs(pho.idhep[1-i])==24&&
934 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])>=11&&
935 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])<=16&&
936 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i+1])>=11&&
937 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i+1])<=16 ){
938 :
939 : if(
940 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])==11||
941 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])==13||
942 0 : abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])==15 ){
943 0 : I=pho.jdahep[1-i][1-i];
944 0 : }
945 : else{
946 0 : I=pho.jdahep[1-i][1-i]+1;
947 : }
948 : //.. muon energy
949 0 : EMU=pho.phep[I-i][4-i];
950 : //.. muon mass square
951 0 : MCHREN=fabs(pho.phep[I-i][4-i]*pho.phep[I-i][4-i]-pho.phep[I-i][3-i]*pho.phep[I-i][3-i]
952 0 : -pho.phep[I-i][2-i]*pho.phep[I-i][2-i]-pho.phep[I-i][1-i]*pho.phep[I-i][1-i]);
953 0 : BETA=sqrt(1- MCHREN/ pho.phep[I-i][4-i]*pho.phep[I-i][4-i]);
954 0 : COSTHG=((pho.phep[I-i][3-i]*pho.phep[pho.nhep-i][3-i]+pho.phep[I-i][2-i]*pho.phep[pho.nhep-i][2-i]
955 0 : +pho.phep[I-i][1-i]*pho.phep[pho.nhep-i][1-i])/
956 0 : sqrt(pho.phep[I-i][3-i]*pho.phep[I-i][3-i]+pho.phep[I-i][2-i]*pho.phep[I-i][2-i]+pho.phep[I-i][1-i]*pho.phep[I-i][1-i]) /
957 0 : sqrt(pho.phep[pho.nhep-i][3-i]*pho.phep[pho.nhep-i][3-i]+pho.phep[pho.nhep-i][2-i]*pho.phep[pho.nhep-i][2-i]+pho.phep[pho.nhep-i][1-i]*pho.phep[pho.nhep-i][1-i]));
958 0 : MPASQR=pho.phep[1-i][4-i]*pho.phep[1-i][4-i];
959 0 : XPH=pho.phep[pho.nhep-i][4-i];
960 :
961 : //... Initialization of the W->l\nu\gamma
962 : //... decay Matrix Element parameters
963 0 : SANC_INIT(phocop_.alpha,phlun);
964 :
965 :
966 0 : MB=pho.phep[1-i][4-i];// ! W boson mass
967 0 : MF2=sqrt(MCHREN);// ! muon mass
968 : I3=-1;
969 0 : for(IJ=1;IJ<=hep.nhep;IJ++){
970 0 : if(abs(hep.idhep[IJ-i])==24){ I3=IJ;} //! position of W
971 : }
972 0 : if(I3==-1) {cout << " ERROR IN PHOBWnlo of PHOTS W-ME: I3= &2i"<<I3<<endl;}
973 : if(
974 0 : abs(hep.idhep[hep.jdahep[I3-i][1-i]-i ])==11||
975 0 : abs(hep.idhep[hep.jdahep[I3-i][1-i]-i ])==13||
976 0 : abs(hep.idhep[hep.jdahep[I3-i][1-i]-i ])==15 ){
977 0 : I4=hep.jdahep[I3-i][1-i];} // ! position of lepton
978 : else{
979 0 : I4=hep.jdahep[I3-i][1-i]+1 ; // ! position of lepton
980 : }
981 :
982 0 : if (hep.idhep[I3-i]==-24) QB=-1.0;// ! W boson charge
983 0 : if (hep.idhep[I3-i]==+24) QB=+1.0;//
984 0 : if (hep.idhep[I4-i]>0.0) QF2=-1.0; // ! lepton charge
985 0 : if (hep.idhep[I4-i]<0.0) QF2=+1.0;
986 :
987 :
988 : //... Particle momenta before foton radiation; effective Born level
989 0 : for( JJ=1; JJ<=4;JJ++){
990 0 : B_PW [(JJ % 4)]=hep.phep[I3-i][JJ-i];// ! W boson
991 0 : B_PNE[(JJ % 4)]=hep.phep[I3-i][JJ-i]-hep.phep[I4-i][JJ-i];// ! neutrino
992 0 : B_PMU[(JJ % 4)]=hep.phep[I4-i][JJ-i]; // ! muon
993 : }
994 :
995 : //.. Particle monenta after photon radiation
996 0 : for( JJ=1; JJ<=4;JJ++){
997 0 : PW [(JJ % 4)]=pho.phep[1-i][JJ-i];
998 0 : PMU [(JJ % 4)]=pho.phep[I-i][JJ-i];
999 0 : PPHOT[(JJ % 4)]=pho.phep[pho.nhep-i][JJ-i];
1000 0 : PNE [(JJ % 4)]=pho.phep[1-i][JJ-i]-pho.phep[I-i][JJ-i]-pho.phep[pho.nhep-i][JJ-i];
1001 : }
1002 :
1003 : // two options of calculating neutrino (spectator) mass
1004 0 : MF1=sqrt(fabs(B_PNE[0]*B_PNE[0]*-B_PNE[1]*B_PNE[1]-B_PNE[2]*B_PNE[2]-B_PNE[3]*B_PNE[3]));
1005 0 : MF1=sqrt(fabs( PNE[0]*PNE[0]- PNE[1]*PNE[1]- PNE[2]*PNE[2]- PNE[3]*PNE[3]));
1006 :
1007 0 : SANC_INIT1(QB,QF2,MF1,MF2,MB);
1008 0 : *WT=(*WT)*SANC_WT(PW,PNE,PMU,PPHOT,B_PW,B_PNE,B_PMU);
1009 0 : }
1010 : // write(*,*) 'AMPSQR/EIKONALFACTOR= ', AMPSQR/EIKONALFACTOR
1011 0 : }
1012 :
1013 : } // namespace Photospp
1014 :
|
