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