forW-MEc.cxx
1 #include "forW-MEc.h"
2 #include "Photos.h"
3 #include "photosC.h"
4 #include <cstdlib>
5 #include<iostream>
6 using std::cout;
7 using std::endl;
8 
9 namespace Photospp
10 {
11 
12 // COMMON /Kleiss_Stirling/spV,bet
13 double PhotosMEforW::spV[4],PhotosMEforW::bet[4];
14 
15 // COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN
16 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;
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 //////////////////////////////////////////////////////////////////
28 complex<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 
79 double 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 
90 complex<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 /////////////////////////////////////////////////////////////////////
126 complex<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 //======================================================================
172 complex<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 //======================================================================
197 complex<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 
233 complex<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 ////////////////////////////////////////////////////////////////
381 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){
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 //======================================================================
617 complex<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 
675 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){
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 
720 double 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 
740 double 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 
766 double 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 
816 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]){
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 
850 void 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 
858 void 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 //----------------------------------------------------------------------
914 void 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 
Photospp::PhotosMEforW::PHOBWnlo
static void PHOBWnlo(double *WT)
Definition: forW-MEc.cxx:914
Generated on Tue Mar 12 2024 23:11:40 for C++InterfacetoPHOTOS by   doxygen 1.8.7