]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/Photos/forW-MEc.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
making script executable
[u/mrichter/AliRoot.git] / TEvtGen / Photos / forW-MEc.cxx
CommitLineData
0ca57c2f 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>
7using std::cout;
8using std::endl;
9
10namespace Photospp
11{
12
13// COMMON /Kleiss_Stirling/spV,bet
14double PhotosMEforW::spV[4],PhotosMEforW::bet[4];
15
16// COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN
17double 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//////////////////////////////////////////////////////////////////
29complex<double> PhotosMEforW::InProd_zero(double p1[4],int l1,double p2[4],int l2){
30
31
32 double forSqrt1,forSqrt2,sqrt1,sqrt2;
33 complex<double> Dcmplx;
34 static complex<double> i_= complex<double>(0.0,1.0);
35 bool equal;
36
37
38
39 equal = true;
40 for (int i = 0; i < 4; i++){
41
42 if (p1[i]!=p2[i]) equal = equal && false ;
43 }
44
45
46 if ( (l1==l2) || equal ) return complex<double>(0.0,0.0);
47
48
49 else if ( (l1==+1) && (l2==-1) ){
50
51 forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
52 forSqrt2 = 1.0/forSqrt1;
53 sqrt1 = sqrt(forSqrt2);
54 sqrt2 = sqrt(forSqrt1);
55
56 return (p1[2]+i_*p1[3])*sqrt1 -
57 (p2[2]+i_*p2[3])*sqrt2 ;
58 }
59 else if ( (l1==-1) && (l2==+1) ){
60
61 forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
62 forSqrt2 = 1.0/forSqrt1;
63 sqrt1 = sqrt(forSqrt2);
64 sqrt2 = sqrt(forSqrt1);
65
66 return (p2[2]-i_*p2[3])*sqrt2 -
67 (p1[2]-i_*p1[3])*sqrt1 ;
68 }
69 else{
70
71
72 cout << " "<<endl;
73 cout << " ERROR IN InProd_zero:"<<endl;
74 cout << " WRONG VALUES FOR l1,l2: l1,l2 = -1,+1 "<<endl;
75 cout << " " <<endl;
76 exit(0);
77 }
78}
79
80double PhotosMEforW::InSqrt(double p[4],double q[4]){
81
82 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
91complex<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 if ((l1==+1)&&(l2==+1)) {
96 forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
97 sqrt1 = sqrt(forSqrt1);
98 sqrt2 = 1.0/sqrt1;
99 return complex<double>(m1*sqrt2+m2*sqrt1,0.0);
100 }
101 else if ((l1==+1)&&(l2==-1))
102 return InProd_zero(p1,+1,p2,-1);
103
104 else if ((l1==-1)&&(l2==+1))
105 return InProd_zero(p1,-1,p2,+1);
106
107 else if ((l1==-1)&&(l2==-1)){
108 forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
109 sqrt1 = sqrt(forSqrt1);
110 sqrt2 = 1.0/sqrt1;
111 return complex<double>(m1*sqrt2+m2*sqrt1,0.0);
112 }
113 else {
114 cout <<" " <<endl;
115 cout <<" ERROR IN InProd_mass.."<<endl;
116 cout <<" WRONG VALUES FOR l1,l2"<<endl;
117 cout <<" " <<endl;
118 exit(0);
119 }
120}
121
122/////////////////////////////////////////////////////////////////////
123// //
124// this is small B_{s}(k,p) function when TrMartix is diaagonal!! //
125// //
126/////////////////////////////////////////////////////////////////////
127complex<double> PhotosMEforW::BsFactor(int s,double k[4],double p[4],double m){
128 double forSqrt1,sqrt1;
129 complex<double> inPr1;
130
131 if ( s==1 ){
132
133 inPr1 = InProd_zero(k,+1,p,-1);
134 forSqrt1 = (p[0]-p[1])/(k[0]-k[1]);
135 sqrt1 = sqrt(2.0*forSqrt1);
136 //BsFactor =
137 return inPr1*sqrt1;
138 }
139
140 else if ( s==-1 ){
141
142 inPr1 = InProd_zero(k,-1,p,+1);
143 forSqrt1 = (p[0]-p[1])/(k[0]-k[1]);
144 sqrt1 = sqrt(2.0*forSqrt1);
145 //BsFactor =
146 return inPr1*sqrt1;
147 }
148 else{
149
150 cout << " "<<endl;
151 cout << " ERROR IN BsFactor: "<<endl;
152 cout << " WRONG VALUES FOR s : s = -1,+1"<<endl;
153 cout << " " <<endl;
154 exit(0);
155 }
156}
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//======================================================================
173complex<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 complex<double> wdecayeikonalks_1ph,BSoft1,BSoft2;
177
178 scalProd1 = p1[0]*k[0]-p1[1]*k[1]-p1[2]*k[2]-p1[3]*k[3];
179 scalProd2 = p2[0]*k[0]-p2[1]*k[1]-p2[2]*k[2]-p2[3]*k[3];
180 scalProd3 = p3[0]*k[0]-p3[1]*k[1]-p3[2]*k[2]-p3[3]*k[3];
181
182
183 BSoft1 = BsFactor(s,k,p1,mf1);
184 BSoft2 = BsFactor(s,k,p2,mf2);
185
186 //WDecayEikonalKS_1ph =
187 return sqrt(pi/alphaI)*(-(qf1/scalProd1+qb/scalProd3)*BSoft1
188 +(qf2/scalProd2-qb/scalProd3)*BSoft2);
189
190}
191
192//======================================================================
193//
194// Gauge invariant soft factor for decay!!
195// Gmass2 -- photon mass square
196//
197//======================================================================
198complex<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 complex<double> BsFactor2,BsFactor1;
202
203
204 ScalProd1 = k[0]*p1[0]-k[1]*p1[1]-k[2]*p1[2]-k[3]*p1[3];
205 ScalProd2 = k[0]*p2[0]-k[1]*p2[1]-k[2]*p2[2]-k[3]*p2[3];
206
207 BsFactor1 = BsFactor(s,k,p1,m1);
208 BsFactor2 = BsFactor(s,k,p2,m2);
209
210 return + BsFactor2/2.0/(ScalProd2-Gmass2)
211 - BsFactor1/2.0/(ScalProd1-Gmass2);
212}
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
234complex<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 double sqrt1,sqrt2; // ,scalProd1,scalProd2;
239 complex<double> inPr1,inPr2,inPr3;
240 bool equal;
241
242 equal = true;
243 for (int i = 0; i < 4; i++)
244 if (p1[i] != p2[i]) equal = equal&&false;
245
246
247
248 if ( (m1==m2)&&(equal) ){
249 //..
250 //.. when: p1=p2=p <=> m1=m2 TrMatrix_zero is diagonal
251 //..
252 if ( (l1==+1)&&(l2==+1) ){
253
254 inPr1 = InProd_zero(k,+s,p1,-s);
255 forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
256 sqrt1 = sqrt(2.0*forSqrt1);
257
258 return sqrt1*inPr1;
259 }
260
261 else if ( (l1==+1)&&(l2==-1) ){
262
263 return complex<double>(0.0,0.0);}
264
265
266 else if ( (l1==-1)&&(l2==+1) ){
267
268 return complex<double>(0.0,0.0);
269 }
270
271 else if ( (l1==-1)&&(l2==-1) ){
272
273 inPr1 = InProd_zero(k,+s,p1,-s);
274 forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
275 sqrt1 = sqrt(2.0*forSqrt1);
276
277 return sqrt1*inPr1;
278 }
279
280 else{
281
282 cout << "" <<endl;
283 cout << " ERROR IN TrMatrix_zero: " <<endl;
284 cout << " WRONG VALUES FOR l1,l2,s" <<endl;
285 cout << "" <<endl;
286 exit(0);
287
288 }
289
290 }
291
292 if ( (l1==+1)&&(l2==+1)&&(s==+1) ){
293
294 inPr1 = InProd_zero(k,+1,p1,-1);
295 forSqrt1 = (p2[0]-p2[1])/(k[0]-k[1]);
296 sqrt1 = sqrt(2.0*forSqrt1);
297
298 return sqrt1*inPr1;
299 }
300 else if ( (l1==+1)&&(l2==-1)&&(s==+1) ) {
301
302 return complex<double>(0.0,0.0);
303 }
304
305 else if( (l1==-1)&&(l2==+1)&&(s==+1) ){
306
307 forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
308 forSqrt2 = 1.0/forSqrt1;
309 sqrt1 = sqrt(2.0*forSqrt1);
310 sqrt2 = sqrt(2.0*forSqrt2);
311
312 return complex<double>(m2*sqrt1-m1*sqrt2,0.0);
313 }
314 else if ( (l1==-1)&&(l2==-1)&&(s==+1) ){
315
316 inPr1 = InProd_zero(k,+1,p2,-1);
317 forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
318 sqrt1 = sqrt(2.0*forSqrt1);
319
320 return inPr1*sqrt1;
321 }
322
323 else if ( (l1==+1)&&(l2==+1)&&(s==-1) ){
324
325 inPr1 = -InProd_zero(k,-1,p2,+1);
326 forSqrt1 = (p1[0]-p1[1])/(k[0]-k[1]);
327 sqrt1 = sqrt(2.0*forSqrt1);
328
329 return -sqrt1*inPr1;
330 }
331
332 else if ( (l1==+1)&&(l2==-1)&&(s==-1) ){
333
334 forSqrt1 = (p1[0]-p1[1])/(p2[0]-p2[1]);
335 forSqrt2 = 1.0/forSqrt1;
336 sqrt1 = sqrt(2.0*forSqrt1);
337 sqrt2 = sqrt(2.0*forSqrt2);
338
339 return complex<double>(m2*sqrt1-m1*sqrt2,0.0);
340 }
341
342 else if ( (l1==-1)&&(l2==+1)&&(s==-1) ){
343
344 return complex<double>(0.0,0.0);
345 }
346
347 else if( (l1==-1)&&(l2==-1)&&(s==-1) ){
348
349 inPr1 = -InProd_zero(k,-1,p1,+1);
350 forSqrt1 = (p2[0]-p2[1])/(k[0]-k[1]);
351 sqrt1 = sqrt(2.0*forSqrt1);
352
353 return -inPr1*sqrt1;
354 }
355 else {
356
357 cout << "" << endl;
358 cout << " ERROR IN TrMatrix_zero: " << endl;
359 cout << " WRONG VALUES FOR l1,l2,s" << endl;
360 cout << "" << endl;
361 exit(0);
362 }
363
364}
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////////////////////////////////////////////////////////////////
382complex<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 double k_1[4],k_2[4];
387 double forSqrt3,forSqrt4,sqrt3,sqrt1,sqrt2,sqrt4;
388 complex<double> inPr1,inPr2,inPr3,inPr4;
389
390 for (int i = 0; i < 4; i++) {
391 k_1[i] = 1.0/2.0*(k[i] - m*spV[i]);
392 k_2[i] = 1.0/2.0*(k[i] + m*spV[i]);
393 }
394
395 if ( (l1==+1)&&(l2==+1)&&(s==0) ){
396
397 inPr1 = InProd_zero(p1,+1,k_2,-1);
398 inPr2 = InProd_zero(p2,-1,k_2,+1);
399 inPr3 = InProd_zero(p1,+1,k_1,-1);
400 inPr4 = InProd_zero(p2,-1,k_1,+1);
401 sqrt1 = sqrt(p1[0]-p1[1]);
402 sqrt2 = sqrt(p2[0]-p2[1]);
403 sqrt3 = m1*m2/sqrt1/sqrt2;
404
405 return
406 (inPr1*inPr2-inPr3*inPr4)*(vf+af)/m
407 + (k_1[0]-k_2[0]-k_1[1]+k_2[1])*sqrt3*(vf-af)/m;
408 }
409
410 else if ( (l1==+1)&&(l2==-1)&&(s==0) ){
411
412 inPr1 = InProd_zero(p1,+1,k_1,-1);
413 inPr2 = InProd_zero(p1,+1,k_2,-1);
414 inPr3 = InProd_zero(p2,+1,k_2,-1);
415 inPr4 = InProd_zero(p2,+1,k_1,-1);
416
417 forSqrt1 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
418 forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
419 forSqrt3 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
420 forSqrt4 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
421 sqrt1 = sqrt(forSqrt1);
422 sqrt2 = sqrt(forSqrt2);
423 sqrt3 = sqrt(forSqrt3);
424 sqrt4 = sqrt(forSqrt4);
425
426 return
427 (inPr1*sqrt1 - inPr2*sqrt2)*(vf+af)*m2/m
428 + (inPr3*sqrt3 - inPr4*sqrt4)*(vf-af)*m1/m;
429 }
430 else if ( (l1==-1)&&(l2==+1)&&(s==0) ){
431
432 inPr1 = InProd_zero(p1,-1,k_1,+1);
433 inPr2 = InProd_zero(p1,-1,k_2,+1);
434 inPr3 = InProd_zero(p2,-1,k_2,+1);
435 inPr4 = InProd_zero(p2,-1,k_1,+1);
436
437 forSqrt1 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
438 forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
439 forSqrt3 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
440 forSqrt4 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
441 sqrt1 = sqrt(forSqrt1);
442 sqrt2 = sqrt(forSqrt2);
443 sqrt3 = sqrt(forSqrt3);
444 sqrt4 = sqrt(forSqrt4);
445
446 return
447 (inPr1*sqrt1 - inPr2*sqrt2)*(vf-af)*m2/m
448 + (inPr3*sqrt3 - inPr4*sqrt4)*(vf+af)*m1/m;
449 }
450 else if ( (l1==-1)&&(l2==-1)&&(s==0) ){
451
452 inPr1 = InProd_zero(p2,+1,k_2,-1);
453 inPr2 = InProd_zero(p1,-1,k_2,+1);
454 inPr3 = InProd_zero(p2,+1,k_1,-1);
455 inPr4 = InProd_zero(p1,-1,k_1,+1);
456 sqrt1 = sqrt(p1[0]-p1[1]);
457 sqrt2 = sqrt(p2[0]-p2[1]);
458 sqrt3 = m1*m2/sqrt1/sqrt2;
459
460 return
461 (inPr1*inPr2 - inPr3*inPr4)*(vf-af)/m
462 + (k_1[0]-k_2[0]-k_1[1]+k_2[1])*sqrt3*(vf+af)/m;
463 }
464 else if ( (l1==+1)&&(l2==+1)&&(s==+1) ){
465
466 inPr1 = InProd_zero(p1,+1,k_1,-1);
467 inPr2 = InProd_zero(k_2,-1,p2,+1);
468 inPr3 = inPr1*inPr2;
469
470 forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
471 forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
472 sqrt1 = sqrt(forSqrt1);
473 sqrt2 = sqrt(forSqrt2);
474 sqrt3 = m1*m2*sqrt1*sqrt2;
475
476 return
477 sqrt(2.0)/m*(inPr3*(vf+af)+sqrt3*(vf-af));
478 }
479
480 else if ( (l1==+1)&&(l2==-1)&&(s==+1) ){
481
482 inPr1 = InProd_zero(p1,+1,k_1,-1);
483 inPr2 = InProd_zero(p2,+1,k_1,-1);
484
485 forSqrt1 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
486 forSqrt2 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
487 sqrt1 = m2*sqrt(forSqrt1);
488 sqrt2 = m1*sqrt(forSqrt2);
489
490 return
491 sqrt(2.0)/m*( + inPr1*sqrt1*(vf+af)
492 - inPr2*sqrt2*(vf-af)
493 );
494 }
495 else if ( (l1==-1)&&(l2==+1)&&(s==+1) ){
496
497 inPr1 = InProd_zero(k_2,-1,p2,+1);
498 inPr2 = InProd_zero(k_2,-1,p1,+1);
499
500 forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
501 forSqrt2 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
502 sqrt1 = m1*sqrt(forSqrt1);
503 sqrt2 = m2*sqrt(forSqrt2);
504
505 return
506 sqrt(2.0)/m*( + inPr1*sqrt1*(vf+af)
507 - inPr2*sqrt2*(vf-af)
508 );
509 }
510 else if ( (l1==-1)&&(l2==-1)&&(s==+1) ){
511
512 inPr1 = InProd_zero(p2,+1,k_1,-1);
513 inPr2 = InProd_zero(k_2,-1,p1,+1);
514 inPr3 = inPr1*inPr2;
515
516 forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
517 forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
518 sqrt1 = sqrt(forSqrt1);
519 sqrt2 = sqrt(forSqrt2);
520 sqrt3 = m1*m2*sqrt1*sqrt2;
521
522 return
523 sqrt(2.0)/m*(inPr3*(vf-af)+sqrt3*(vf+af));
524 }
525
526 else if ( (l1==+1)&&(l2==+1)&&(s==-1) ){
527
528 inPr1 = InProd_zero(p2,-1,k_1,+1);
529 inPr2 = InProd_zero(k_2,+1,p1,-1);
530 inPr3 = inPr1*inPr2;
531
532 forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
533 forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
534 sqrt1 = sqrt(forSqrt1);
535 sqrt2 = sqrt(forSqrt2);
536 sqrt3 = m1*m2*sqrt1*sqrt2;
537
538 return
539 sqrt(2.0)/m*(inPr3*(vf+af)+sqrt3*(vf-af));
540 }
541 else if ( (l1==+1)&&(l2==-1)&&(s==-1) ){
542
543 inPr1 = InProd_zero(k_2,+1,p2,-1);
544 inPr2 = InProd_zero(k_2,+1,p1,-1);
545
546 forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
547 forSqrt2 = (k_1[0]-k_1[1])/(p2[0]-p2[1]);
548 sqrt1 = m1*sqrt(forSqrt1);
549 sqrt2 = m2*sqrt(forSqrt2);
550
551 return
552 sqrt(2.0)/m*(+ inPr1*sqrt1*(vf-af)
553 - inPr2*sqrt2*(vf+af)
554 );
555 }
556 else if ( (l1==-1)&&(l2==+1)&&(s==-1) ){
557
558 inPr1 = InProd_zero(p1,-1,k_1,+1);
559 inPr2 = InProd_zero(p2,-1,k_1,+1);
560
561 forSqrt1 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
562 forSqrt2 = (k_2[0]-k_2[1])/(p1[0]-p1[1]);
563 sqrt1 = m2*sqrt(forSqrt1);
564 sqrt2 = m1*sqrt(forSqrt2);
565
566 return
567 sqrt(2.0)/m*(+ inPr1*sqrt1*(vf-af)
568 - inPr2*sqrt2*(vf+af)
569 );
570 }
571 else if ( (l1==-1)&&(l2==-1)&&(s==-1) ){
572
573 inPr1 = InProd_zero(p1,-1,k_1,+1);
574 inPr2 = InProd_zero(k_2,+1,p2,-1);
575 inPr3 = inPr1*inPr2;
576
577 forSqrt1 = (k_1[0]-k_1[1])/(p1[0]-p1[1]);
578 forSqrt2 = (k_2[0]-k_2[1])/(p2[0]-p2[1]);
579 sqrt1 = sqrt(forSqrt1);
580 sqrt2 = sqrt(forSqrt2);
581 sqrt3 = m1*m2*sqrt1*sqrt2;
582
583 return
584 sqrt(2.0)/m*(inPr3*(vf-af)+sqrt3*(vf+af));
585 }
586
587 else{
588
589 cout << " "<< endl;
590 cout << " TrMatrix_mass: Wrong values for l1,l2,s:"<< endl;
591 cout << " l1,l2 = -1,+1; s = -1,0,1 "<< endl;
592 cout << " "<< endl;
593 exit(0);
594
595 }
596
597}
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//======================================================================
618complex<double> PhotosMEforW::WDecayBornAmpKS_1ph(double p3[4],int l3,double p1[4],int l1,double p2[4],int l2){
619
620 double coeff;
621
622
623 coeff = sqrt(pi/alphaI/2.0)/sw; // vertex: g/2/sqrt(2)
624
625 return coeff*TrMatrix_mass(p2,mf2,l2,p3,mb,l3,p1,-mf1,-l1);
626}
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
676complex<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 double scalProd1,scalProd2,scalProd3,coeff; //,theta3,ph3;
679 complex<double> bornAmp,TrMx1,TrMx2;
680 complex<double> BSoft1,BSoft2;
681
682 coeff = sqrt(2.0)*pi/sw/alphaI; // vertex: g/2/sqrt[2] * e
683
684 scalProd1 = p1[0]*k[0]-p1[1]*k[1]-p1[2]*k[2]-p1[3]*k[3];
685 scalProd2 = p2[0]*k[0]-p2[1]*k[1]-p2[2]*k[2]-p2[3]*k[3];
686 scalProd3 = p3[0]*k[0]-p3[1]*k[1]-p3[2]*k[2]-p3[3]*k[3];
687
688 BSoft1 = BsFactor(s,k,p1,mf1);
689 BSoft2 = BsFactor(s,k,p2,mf2);
690 bornAmp = TrMatrix_mass(p2,mf2,l2,p3,mb,l3,p1,-mf1,-l1);
691 TrMx1 = complex<double>(0.0,0.0);
692 TrMx2 = complex<double>(0.0,0.0);
693
694 for (int la1 = -1; la1< 3 ; la1+=2) {
695 // DO la1=-1,1,2
696 TrMx1 = TrMx1 + TrMatrix_zero(k,0.0,-la1,k,s,p1,-mf1,-l1)*
697 TrMatrix_mass(p2,mf2,l2,p3,mb,l3,k,0.0,-la1);
698 TrMx2 = TrMx2 + TrMatrix_zero(p2,mf2,l2,k,s,k,0.0,la1)*
699 TrMatrix_mass(k,0.0,la1,p3,mb,l3,p1,-mf1,-l1);
700 }
701
702 return coeff * (
703 + (-(qf1/scalProd1+qb/scalProd3)*BSoft1 // IR-divergent part of amplitude
704 +(qf2/scalProd2-qb/scalProd3)*BSoft2)/2.0*bornAmp
705 //
706 - (qf1/scalProd1+qb/scalProd3)*TrMx1/2.0 // IR-finite part of amplitude
707 + (qf2/scalProd2-qb/scalProd3)*TrMx2/2.0
708 );
709}
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
721double PhotosMEforW::WDecayEikonalSqrKS_1ph(double p3[4],double p1[4],double p2[4],double k[4]){
722 double spinSumAvrg;
723 complex<double> wDecAmp;
724
725 spinSumAvrg = 0.0;
726 for (int s = -1; s< 3 ; s+=2) {
727 wDecAmp = WDecayEikonalKS_1ph(p3,p1,p2,k,s);
728 spinSumAvrg = spinSumAvrg + real(wDecAmp*conj(wDecAmp));
729 }
730 return spinSumAvrg;
731}
732
733//========================================================
734// The squared eikonal factor for W decay =
735// into fermion pair and one photon =
736// INPUT : =
737// =
738// OUTPUT: =
739//========================================================
740
741double PhotosMEforW::WDecayBornAmpSqrKS_1ph(double p3[4],double p1[4],double p2[4]){
742 double spinSumAvrg;
743 complex<double> wDecAmp;
744
745 spinSumAvrg = 0.0;
746 for (int l3 = -1; l3< 2 ; l3++) {
747 for (int l1 = -1; l1< 3 ; l1+=2) {
748 for (int l2 = -1; l2< 3 ; l2+=2) {
749 wDecAmp = WDecayBornAmpKS_1ph(p3,l3,p1,l1,p2,l2);
750 spinSumAvrg = spinSumAvrg + real(wDecAmp*conj(wDecAmp));
751 }
752 }
753 }
754 return spinSumAvrg;
755}
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
767double PhotosMEforW::WDecayAmplitudeSqrKS_1ph(double p3[4],double p1[4],double p2[4], double k[4]){
768
769 double spinSumAvrg;
770 complex<double> wDecAmp;
771
772 spinSumAvrg = 0.0;
773 for (int l3 = -1; l3< 2 ; l3++) {
774 for (int l1 = -1; l1< 3 ; l1+=2) {
775 for (int l2 = -1; l2< 3 ; l2+=2) {
776 for (int s = -1; s < 3 ; s+=2) {
777 wDecAmp = WDecayAmplitudeKS_1ph(p3,l3,p1,l1,p2,l2,k,s);
778 spinSumAvrg = spinSumAvrg + real(wDecAmp*conj(wDecAmp));
779 }
780 }
781 }
782 }
783 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}
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
817double 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 double AMPSQR=WDecayAmplitudeSqrKS_1ph(PW,PNE,PMU,PPHOT);
822
823 double EIKONALFACTOR=WDecayBornAmpSqrKS_1ph(B_PW,B_PNE,B_PMU)
824 *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 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
851void PhotosMEforW::SANC_INIT1(double QB0,double QF20,double MF10,double MF20,double MB0){
852 qb =QB0;
853 qf2=QF20;
854 mf1=MF10;
855 mf2=MF20;
856 mb =MB0;
857}
858
859void 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 if (SANC_MC_INIT==-123456789){
867 SANC_MC_INIT=1;
868
869 pi=4*atan(1.0);
870 qf1=0.0; // neutrino charge
871 mf1=1.0e-10; // newutrino mass
872 vf=1.0; // V&A couplings
873 af=1.0;
874 alphaI=1.0/ALPHA;
875 cw=0.881731727; // Weak Weinberg angle
876 sw=0.471751166;
877
878
879 //... An auxilary K&S vectors
880 bet[0]= 1.0;
881 bet[1]= 0.0722794881816159;
882 bet[2]=-0.994200045099866;
883 bet[3]= 0.0796363353729248;
884
885 spV[0]= 0.0;
886 spV[1]= 7.22794881816159e-2;
887 spV[2]=-0.994200045099866;
888 spV[3]= 7.96363353729248e-2;
889
890 mcLUN = PHLUN;
891 }
892}
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//----------------------------------------------------------------------
915void 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 double PW[4],PMU[4],PPHOT[4],PNE[4];
921 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 if(abs(pho.idhep[1-i])==24&&
934 abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])>=11&&
935 abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])<=16&&
936 abs(pho.idhep[pho.jdahep[1-i][1-i]-i+1])>=11&&
937 abs(pho.idhep[pho.jdahep[1-i][1-i]-i+1])<=16 ){
938
939 if(
940 abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])==11||
941 abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])==13||
942 abs(pho.idhep[pho.jdahep[1-i][1-i]-i ])==15 ){
943 I=pho.jdahep[1-i][1-i];
944 }
945 else{
946 I=pho.jdahep[1-i][1-i]+1;
947 }
948 //.. muon energy
949 EMU=pho.phep[I-i][4-i];
950 //.. muon mass square
951 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 -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 BETA=sqrt(1- MCHREN/ pho.phep[I-i][4-i]*pho.phep[I-i][4-i]);
954 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 +pho.phep[I-i][1-i]*pho.phep[pho.nhep-i][1-i])/
956 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 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 MPASQR=pho.phep[1-i][4-i]*pho.phep[1-i][4-i];
959 XPH=pho.phep[pho.nhep-i][4-i];
960
961 //... Initialization of the W->l\nu\gamma
962 //... decay Matrix Element parameters
963 SANC_INIT(phocop_.alpha,phlun);
964
965
966 MB=pho.phep[1-i][4-i];// ! W boson mass
967 MF2=sqrt(MCHREN);// ! muon mass
968 I3=-1;
969 for(IJ=1;IJ<=hep.nhep;IJ++){
970 if(abs(hep.idhep[IJ-i])==24){ I3=IJ;} //! position of W
971 }
972 if(I3==-1) {cout << " ERROR IN PHOBWnlo of PHOTS W-ME: I3= &2i"<<I3<<endl;}
973 if(
974 abs(hep.idhep[hep.jdahep[I3-i][1-i]-i ])==11||
975 abs(hep.idhep[hep.jdahep[I3-i][1-i]-i ])==13||
976 abs(hep.idhep[hep.jdahep[I3-i][1-i]-i ])==15 ){
977 I4=hep.jdahep[I3-i][1-i];} // ! position of lepton
978 else{
979 I4=hep.jdahep[I3-i][1-i]+1 ; // ! position of lepton
980 }
981
982 if (hep.idhep[I3-i]==-24) QB=-1.0;// ! W boson charge
983 if (hep.idhep[I3-i]==+24) QB=+1.0;//
984 if (hep.idhep[I4-i]>0.0) QF2=-1.0; // ! lepton charge
985 if (hep.idhep[I4-i]<0.0) QF2=+1.0;
986
987
988 //... Particle momenta before foton radiation; effective Born level
989 for( JJ=1; JJ<=4;JJ++){
990 B_PW [(JJ % 4)]=hep.phep[I3-i][JJ-i];// ! W boson
991 B_PNE[(JJ % 4)]=hep.phep[I3-i][JJ-i]-hep.phep[I4-i][JJ-i];// ! neutrino
992 B_PMU[(JJ % 4)]=hep.phep[I4-i][JJ-i]; // ! muon
993 }
994
995 //.. Particle monenta after photon radiation
996 for( JJ=1; JJ<=4;JJ++){
997 PW [(JJ % 4)]=pho.phep[1-i][JJ-i];
998 PMU [(JJ % 4)]=pho.phep[I-i][JJ-i];
999 PPHOT[(JJ % 4)]=pho.phep[pho.nhep-i][JJ-i];
1000 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 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 MF1=sqrt(fabs( PNE[0]*PNE[0]- PNE[1]*PNE[1]- PNE[2]*PNE[2]- PNE[3]*PNE[3]));
1006
1007 SANC_INIT1(QB,QF2,MF1,MF2,MB);
1008 *WT=(*WT)*SANC_WT(PW,PNE,PMU,PPHOT,B_PW,B_PNE,B_PMU);
1009 }
1010 // write(*,*) 'AMPSQR/EIKONALFACTOR= ', AMPSQR/EIKONALFACTOR
1011}
1012
1013} // namespace Photospp
1014