Belle II Software  release-06-00-14
forW-MEc.cc
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 // COMMON /Kleiss_Stirling/spV,bet
12  double PhotosMEforW::spV[4], PhotosMEforW::bet[4];
13 
14 // COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN
15  double PhotosMEforW::pi, PhotosMEforW::sw, PhotosMEforW::cw, PhotosMEforW::alphaI, PhotosMEforW::qb, PhotosMEforW::mb,
16  PhotosMEforW::mf1, PhotosMEforW::mf2, PhotosMEforW::qf1, PhotosMEforW::qf2, PhotosMEforW::vf, PhotosMEforW::af, PhotosMEforW::mcLUN;
17 
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 // //
28  complex<double> PhotosMEforW::InProd_zero(double p1[4], int l1, double p2[4], int l2)
29  {
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  } 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  } else {
68 
69 
70  cout << " " << endl;
71  cout << " ERROR IN InProd_zero:" << endl;
72  cout << " WRONG VALUES FOR l1,l2: l1,l2 = -1,+1 " << endl;
73  cout << " " << endl;
74  exit(-1);
75  }
76  }
77 
78  double PhotosMEforW::InSqrt(double p[4], double q[4])
79  {
80 
81  return sqrt((p[0] - p[1]) / (q[0] - q[1]));
82  }
83 
85 // //
86 // Inner product for massive spinors: Ub(p1,m1,l1)*U(p2,m2,l2) //
87 // //
89 
90  complex<double> PhotosMEforW::InProd_mass(double p1[4], double m1, int l1, double p2[4], double m2, int l2)
91  {
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  } 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  } else {
112  cout << " " << endl;
113  cout << " ERROR IN InProd_mass.." << endl;
114  cout << " WRONG VALUES FOR l1,l2" << endl;
115  cout << " " << endl;
116  exit(-1);
117  }
118  }
119 
121 // //
122 // this is small B_{s}(k,p) function when TrMartix is diaagonal!! //
123 // //
125  complex<double> PhotosMEforW::BsFactor(int s, double k[4], double p[4], double m)
126  {
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  } else {
147 
148  cout << " " << endl;
149  cout << " ERROR IN BsFactor: " << endl;
150  cout << " WRONG VALUES FOR s : s = -1,+1" << endl;
151  cout << " " << endl;
152  exit(-1);
153  }
154  }
155 
156 
157 
158 
159 //======================================================================
160 //
161 // Eikonal factor of decay W->l_1+l_2+\gamma in terms of K&S objects !
162 //
163 // EikFactor = q1*eps.p1/k.p1 + q2*eps.p2/k.p2 - q3*eps.p3/k.p3
164 //
165 // indices 1,2 are for charged decay products
166 // index 3 is for W
167 //
168 // q - charge
169 //
170 //======================================================================
171  complex<double> PhotosMEforW::WDecayEikonalKS_1ph(double p3[4], double p1[4], double p2[4], double k[4], int s)
172  {
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 
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 
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 
237  double forSqrt1, forSqrt2;
238  // double p1_1[4],p2_1[4];
239  double sqrt1, sqrt2; // ,scalProd1,scalProd2;
240  complex<double> inPr1, inPr2, inPr3;
241  bool equal;
242 
243  equal = true;
244  for (int i = 0; i < 4; i++)
245  if (p1[i] != p2[i]) equal = equal && false;
246 
247 
248 
249  if ((m1 == m2) && (equal)) {
250  //..
251  //.. when: p1=p2=p <=> m1=m2 TrMatrix_zero is diagonal
252  //..
253  if ((l1 == +1) && (l2 == +1)) {
254 
255  inPr1 = InProd_zero(k, +s, p1, -s);
256  forSqrt1 = (p1[0] - p1[1]) / (k[0] - k[1]);
257  sqrt1 = sqrt(2.0 * forSqrt1);
258 
259  return sqrt1 * inPr1;
260  }
261 
262  else if ((l1 == +1) && (l2 == -1)) {
263 
264  return complex<double>(0.0, 0.0);
265  }
266 
267 
268  else if ((l1 == -1) && (l2 == +1)) {
269 
270  return complex<double>(0.0, 0.0);
271  }
272 
273  else if ((l1 == -1) && (l2 == -1)) {
274 
275  inPr1 = InProd_zero(k, +s, p1, -s);
276  forSqrt1 = (p1[0] - p1[1]) / (k[0] - k[1]);
277  sqrt1 = sqrt(2.0 * forSqrt1);
278 
279  return sqrt1 * inPr1;
280  }
281 
282  else {
283 
284  cout << "" << endl;
285  cout << " ERROR IN TrMatrix_zero: " << endl;
286  cout << " WRONG VALUES FOR l1,l2,s" << endl;
287  cout << "" << endl;
288  exit(-1);
289 
290  }
291 
292  }
293 
294  if ((l1 == +1) && (l2 == +1) && (s == +1)) {
295 
296  inPr1 = InProd_zero(k, +1, p1, -1);
297  forSqrt1 = (p2[0] - p2[1]) / (k[0] - k[1]);
298  sqrt1 = sqrt(2.0 * forSqrt1);
299 
300  return sqrt1 * inPr1;
301  } else if ((l1 == +1) && (l2 == -1) && (s == +1)) {
302 
303  return complex<double>(0.0, 0.0);
304  }
305 
306  else if ((l1 == -1) && (l2 == +1) && (s == +1)) {
307 
308  forSqrt1 = (p1[0] - p1[1]) / (p2[0] - p2[1]);
309  forSqrt2 = 1.0 / forSqrt1;
310  sqrt1 = sqrt(2.0 * forSqrt1);
311  sqrt2 = sqrt(2.0 * forSqrt2);
312 
313  return complex<double>(m2 * sqrt1 - m1 * sqrt2, 0.0);
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  } 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 
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 // //
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,
382  int l2)
383  {
384 
385 
386  double forSqrt1, forSqrt2;
387  double k_1[4], k_2[4];
388  double forSqrt3, forSqrt4, sqrt3, sqrt1, sqrt2, sqrt4;
389  complex<double> inPr1, inPr2, inPr3, inPr4;
390 
391  for (int i = 0; i < 4; i++) {
392  k_1[i] = 1.0 / 2.0 * (k[i] - m * spV[i]);
393  k_2[i] = 1.0 / 2.0 * (k[i] + m * spV[i]);
394  }
395 
396  if ((l1 == +1) && (l2 == +1) && (s == 0)) {
397 
398  inPr1 = InProd_zero(p1, +1, k_2, -1);
399  inPr2 = InProd_zero(p2, -1, k_2, +1);
400  inPr3 = InProd_zero(p1, +1, k_1, -1);
401  inPr4 = InProd_zero(p2, -1, k_1, +1);
402  sqrt1 = sqrt(p1[0] - p1[1]);
403  sqrt2 = sqrt(p2[0] - p2[1]);
404  sqrt3 = m1 * m2 / sqrt1 / sqrt2;
405 
406  return
407  (inPr1 * inPr2 - inPr3 * inPr4) * (vf + af) / m
408  + (k_1[0] - k_2[0] - k_1[1] + k_2[1]) * sqrt3 * (vf - af) / m;
409  }
410 
411  else if ((l1 == +1) && (l2 == -1) && (s == 0)) {
412 
413  inPr1 = InProd_zero(p1, +1, k_1, -1);
414  inPr2 = InProd_zero(p1, +1, k_2, -1);
415  inPr3 = InProd_zero(p2, +1, k_2, -1);
416  inPr4 = InProd_zero(p2, +1, k_1, -1);
417 
418  forSqrt1 = (k_1[0] - k_1[1]) / (p2[0] - p2[1]);
419  forSqrt2 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
420  forSqrt3 = (k_2[0] - k_2[1]) / (p1[0] - p1[1]);
421  forSqrt4 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
422  sqrt1 = sqrt(forSqrt1);
423  sqrt2 = sqrt(forSqrt2);
424  sqrt3 = sqrt(forSqrt3);
425  sqrt4 = sqrt(forSqrt4);
426 
427  return
428  (inPr1 * sqrt1 - inPr2 * sqrt2) * (vf + af) * m2 / m
429  + (inPr3 * sqrt3 - inPr4 * sqrt4) * (vf - af) * m1 / m;
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  } 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  } else if ((l1 == +1) && (l2 == +1) && (s == +1)) {
463 
464  inPr1 = InProd_zero(p1, +1, k_1, -1);
465  inPr2 = InProd_zero(k_2, -1, p2, +1);
466  inPr3 = inPr1 * inPr2;
467 
468  forSqrt1 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
469  forSqrt2 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
470  sqrt1 = sqrt(forSqrt1);
471  sqrt2 = sqrt(forSqrt2);
472  sqrt3 = m1 * m2 * sqrt1 * sqrt2;
473 
474  return
475  sqrt(2.0) / m * (inPr3 * (vf + af) + sqrt3 * (vf - af));
476  }
477 
478  else if ((l1 == +1) && (l2 == -1) && (s == +1)) {
479 
480  inPr1 = InProd_zero(p1, +1, k_1, -1);
481  inPr2 = InProd_zero(p2, +1, k_1, -1);
482 
483  forSqrt1 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
484  forSqrt2 = (k_2[0] - k_2[1]) / (p1[0] - p1[1]);
485  sqrt1 = m2 * sqrt(forSqrt1);
486  sqrt2 = m1 * sqrt(forSqrt2);
487 
488  return
489  sqrt(2.0) / m * (+ inPr1 * sqrt1 * (vf + af)
490  - inPr2 * sqrt2 * (vf - af)
491  );
492  } else if ((l1 == -1) && (l2 == +1) && (s == +1)) {
493 
494  inPr1 = InProd_zero(k_2, -1, p2, +1);
495  inPr2 = InProd_zero(k_2, -1, p1, +1);
496 
497  forSqrt1 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
498  forSqrt2 = (k_1[0] - k_1[1]) / (p2[0] - p2[1]);
499  sqrt1 = m1 * sqrt(forSqrt1);
500  sqrt2 = m2 * sqrt(forSqrt2);
501 
502  return
503  sqrt(2.0) / m * (+ inPr1 * sqrt1 * (vf + af)
504  - inPr2 * sqrt2 * (vf - af)
505  );
506  } else if ((l1 == -1) && (l2 == -1) && (s == +1)) {
507 
508  inPr1 = InProd_zero(p2, +1, k_1, -1);
509  inPr2 = InProd_zero(k_2, -1, p1, +1);
510  inPr3 = inPr1 * inPr2;
511 
512  forSqrt1 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
513  forSqrt2 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
514  sqrt1 = sqrt(forSqrt1);
515  sqrt2 = sqrt(forSqrt2);
516  sqrt3 = m1 * m2 * sqrt1 * sqrt2;
517 
518  return
519  sqrt(2.0) / m * (inPr3 * (vf - af) + sqrt3 * (vf + af));
520  }
521 
522  else if ((l1 == +1) && (l2 == +1) && (s == -1)) {
523 
524  inPr1 = InProd_zero(p2, -1, k_1, +1);
525  inPr2 = InProd_zero(k_2, +1, p1, -1);
526  inPr3 = inPr1 * inPr2;
527 
528  forSqrt1 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
529  forSqrt2 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
530  sqrt1 = sqrt(forSqrt1);
531  sqrt2 = sqrt(forSqrt2);
532  sqrt3 = m1 * m2 * sqrt1 * sqrt2;
533 
534  return
535  sqrt(2.0) / m * (inPr3 * (vf + af) + sqrt3 * (vf - af));
536  } else if ((l1 == +1) && (l2 == -1) && (s == -1)) {
537 
538  inPr1 = InProd_zero(k_2, +1, p2, -1);
539  inPr2 = InProd_zero(k_2, +1, p1, -1);
540 
541  forSqrt1 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
542  forSqrt2 = (k_1[0] - k_1[1]) / (p2[0] - p2[1]);
543  sqrt1 = m1 * sqrt(forSqrt1);
544  sqrt2 = m2 * sqrt(forSqrt2);
545 
546  return
547  sqrt(2.0) / m * (+ inPr1 * sqrt1 * (vf - af)
548  - inPr2 * sqrt2 * (vf + af)
549  );
550  } else if ((l1 == -1) && (l2 == +1) && (s == -1)) {
551 
552  inPr1 = InProd_zero(p1, -1, k_1, +1);
553  inPr2 = InProd_zero(p2, -1, k_1, +1);
554 
555  forSqrt1 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
556  forSqrt2 = (k_2[0] - k_2[1]) / (p1[0] - p1[1]);
557  sqrt1 = m2 * sqrt(forSqrt1);
558  sqrt2 = m1 * sqrt(forSqrt2);
559 
560  return
561  sqrt(2.0) / m * (+ inPr1 * sqrt1 * (vf - af)
562  - inPr2 * sqrt2 * (vf + af)
563  );
564  } else if ((l1 == -1) && (l2 == -1) && (s == -1)) {
565 
566  inPr1 = InProd_zero(p1, -1, k_1, +1);
567  inPr2 = InProd_zero(k_2, +1, p2, -1);
568  inPr3 = inPr1 * inPr2;
569 
570  forSqrt1 = (k_1[0] - k_1[1]) / (p1[0] - p1[1]);
571  forSqrt2 = (k_2[0] - k_2[1]) / (p2[0] - p2[1]);
572  sqrt1 = sqrt(forSqrt1);
573  sqrt2 = sqrt(forSqrt2);
574  sqrt3 = m1 * m2 * sqrt1 * sqrt2;
575 
576  return
577  sqrt(2.0) / m * (inPr3 * (vf - af) + sqrt3 * (vf + af));
578  }
579 
580  else {
581 
582  cout << " " << endl;
583  cout << " TrMatrix_mass: Wrong values for l1,l2,s:" << endl;
584  cout << " l1,l2 = -1,+1; s = -1,0,1 " << endl;
585  cout << " " << endl;
586  exit(-1);
587 
588  }
589 
590  }
591 
592 
593 
594 //======================================================================
595 // =
596 // p1,mf1,l1 =
597 // / =
598 // \/_ =
599 // / =
600 // p3,mb,l3 / =
601 // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
602 // \ =
603 // _\/ =
604 // \ =
605 // p2,mf2,l2 =
606 // INPUT : p1,m1,l1; p2,m2,l2; p3,m3,l3 -- momenta,mass and helicity =
607 // =
608 // OUTPUT: value of functions -- decay amplitude =
609 // =
610 //======================================================================
611  complex<double> PhotosMEforW::WDecayBornAmpKS_1ph(double p3[4], int l3, double p1[4], int l1, double p2[4], int l2)
612  {
613 
614  double coeff;
615 
616 
617  coeff = sqrt(pi / alphaI / 2.0) / sw; // vertex: g/2/sqrt(2)
618 
619  return coeff * TrMatrix_mass(p2, mf2, l2, p3, mb, l3, p1, -mf1, -l1);
620  }
621 
622 
623 //======================================================================
624 // k,0,l =
625 // \ p1,mf1,l1 =
626 // / / =
627 // \ \/_ =
628 // / / =
629 // p3,mb,l3 \ / =
630 // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
631 // \ =
632 // _\/ =
633 // \ =
634 // p2,mf2,l2 =
635 // { + } =
636 // p1,mf1,l1 =
637 // / =
638 // \/_~~~~~~~ k,0,s =
639 // / =
640 // p3,mb,l3 / =
641 // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
642 // \ =
643 // _\/ =
644 // \ =
645 // p2,mf2,l2 =
646 // { + } =
647 // p1,mf1,l1 =
648 // / =
649 // \/_ =
650 // / =
651 // p3,mb,l3 / =
652 // \/\/\/\/\/\ ------> g_(mu,1)*(1+g5_(1)) =
653 // \ =
654 // _\/ ~~~~~~~ k,0,s =
655 // \ =
656 // p2,mf2,l2 =
657 // =
658 // all momentas, exept k are incoming !!! =
659 // =
660 // This function culculates The W-ff\gamma decay amplitude into permion=
661 // pair and one photon using Kleisse&Stirling method for helicity =
662 // amplitudes, which includes three above feynman diagramms.. =
663 // =
664 // INPUT : p1,m1,l1; p2,m2,l2; p3,m3,l3 -- momenta,mass and helicity =
665 // =
666 // OUTPUT: value of functions -- decay amplitude =
667 // =
668 //======================================================================
669 
670  complex<double> PhotosMEforW::WDecayAmplitudeKS_1ph(double p3[4], int l3, double p1[4], int l1, double p2[4], int l2, double k[4],
671  int s)
672  {
673 
674  double scalProd1, scalProd2, scalProd3, coeff; //,theta3,ph3;
675  complex<double> bornAmp, TrMx1, TrMx2;
676  complex<double> BSoft1, BSoft2;
677 
678  coeff = sqrt(2.0) * pi / sw / alphaI; // vertex: g/2/sqrt[2] * e
679 
680  scalProd1 = p1[0] * k[0] - p1[1] * k[1] - p1[2] * k[2] - p1[3] * k[3];
681  scalProd2 = p2[0] * k[0] - p2[1] * k[1] - p2[2] * k[2] - p2[3] * k[3];
682  scalProd3 = p3[0] * k[0] - p3[1] * k[1] - p3[2] * k[2] - p3[3] * k[3];
683 
684  BSoft1 = BsFactor(s, k, p1, mf1);
685  BSoft2 = BsFactor(s, k, p2, mf2);
686  bornAmp = TrMatrix_mass(p2, mf2, l2, p3, mb, l3, p1, -mf1, -l1);
687  TrMx1 = complex<double>(0.0, 0.0);
688  TrMx2 = complex<double>(0.0, 0.0);
689 
690  for (int la1 = -1; la1 < 3 ; la1 += 2) {
691  // DO la1=-1,1,2
692  TrMx1 = TrMx1 + TrMatrix_zero(k, 0.0, -la1, k, s, p1, -mf1, -l1) *
693  TrMatrix_mass(p2, mf2, l2, p3, mb, l3, k, 0.0, -la1);
694  TrMx2 = TrMx2 + TrMatrix_zero(p2, mf2, l2, k, s, k, 0.0, la1) *
695  TrMatrix_mass(k, 0.0, la1, p3, mb, l3, p1, -mf1, -l1);
696  }
697 
698  return coeff * (
699  + (-(qf1 / scalProd1 + qb / scalProd3) * BSoft1 // IR-divergent part of amplitude
700  + (qf2 / scalProd2 - qb / scalProd3) * BSoft2) / 2.0 * bornAmp
701  //
702  - (qf1 / scalProd1 + qb / scalProd3) * TrMx1 / 2.0 // IR-finite part of amplitude
703  + (qf2 / scalProd2 - qb / scalProd3) * TrMx2 / 2.0
704  );
705  }
706 
707 
708 
709 //========================================================
710 // The squared eikonal factor for W decay =
711 // into fermion pair and one photon =
712 // INPUT : =
713 // =
714 // OUTPUT: =
715 //========================================================
716 
717  double PhotosMEforW::WDecayEikonalSqrKS_1ph(double p3[4], double p1[4], double p2[4], double k[4])
718  {
719  double spinSumAvrg;
720  complex<double> wDecAmp;
721 
722  spinSumAvrg = 0.0;
723  for (int s = -1; s < 3 ; s += 2) {
724  wDecAmp = WDecayEikonalKS_1ph(p3, p1, p2, k, s);
725  spinSumAvrg = spinSumAvrg + real(wDecAmp * conj(wDecAmp));
726  }
727  return spinSumAvrg;
728  }
729 
730 //========================================================
731 // The squared eikonal factor for W decay =
732 // into fermion pair and one photon =
733 // INPUT : =
734 // =
735 // OUTPUT: =
736 //========================================================
737 
738  double PhotosMEforW::WDecayBornAmpSqrKS_1ph(double p3[4], double p1[4], double p2[4])
739  {
740  double spinSumAvrg;
741  complex<double> wDecAmp;
742 
743  spinSumAvrg = 0.0;
744  for (int l3 = -1; l3 < 2 ; l3++) {
745  for (int l1 = -1; l1 < 3 ; l1 += 2) {
746  for (int l2 = -1; l2 < 3 ; l2 += 2) {
747  wDecAmp = WDecayBornAmpKS_1ph(p3, l3, p1, l1, p2, l2);
748  spinSumAvrg = spinSumAvrg + real(wDecAmp * conj(wDecAmp));
749  }
750  }
751  }
752  return spinSumAvrg;
753  }
754 
755 
756 
757 //========================================================
758 // The squared amplitude for W decay =
759 // into fermion pair and one photon =
760 // INPUT : =
761 // =
762 // OUTPUT: =
763 //========================================================
764 
765  double PhotosMEforW::WDecayAmplitudeSqrKS_1ph(double p3[4], double p1[4], double p2[4], double k[4])
766  {
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],
817  double B_PMU[4])
818  {
819 
820 
821  //.. Exact amplitude square
822  double AMPSQR = WDecayAmplitudeSqrKS_1ph(PW, PNE, PMU, PPHOT);
823 
824  double EIKONALFACTOR = WDecayBornAmpSqrKS_1ph(B_PW, B_PNE, B_PMU)
825  * WDecayEikonalSqrKS_1ph(PW, PNE, PMU, PPHOT);
826 
827  //.. New weight
828 
829  // cout << 'B_pne=',B_PNE << endl;
830  // cout << 'B_PMU=',B_PMU << endl;
831  // cout << 'bornie=',WDecayBornAmpSqrKS_1ph(B_PW,B_PNE,B_PMU) << endl;
832 
833  // cout << ' ' << endl;
834  // cout << ' pne=',pne << endl;
835  // cout << ' pmu=',pmu << endl;
836  // cout << 'pphot=',pphot << endl;
837  // cout << ' ' << endl;
838  // cout << ' b_pw=',B_PW << endl;
839  // cout << ' b_pne=',B_PNE << endl;
840  // cout << 'b_pmu=',B_PMU << endl;
841 
842  // cout << 'cori=',AMPSQR/EIKONALFACTOR,AMPSQR,EIKONALFACTOR << endl;
843 
844  return AMPSQR / EIKONALFACTOR;
845  //
846  // return (1-8*EMU*XPH*(1-COSTHG*BETA)*
847  // (MCHREN+2*XPH*SQRT(MPASQR))/
848  // MPASQR**2/(1-MCHREN/MPASQR)/(4-MCHREN/MPASQR))
849  }
850 
851 
852  void PhotosMEforW::SANC_INIT1(double QB0, double QF20, double MF10, double MF20, double MB0)
853  {
854  qb = QB0;
855  qf2 = QF20;
856  mf1 = MF10;
857  mf2 = MF20;
858  mb = MB0;
859  }
860 
861  void PhotosMEforW::SANC_INIT(double ALPHA, int PHLUN)
862  {
863 
864 
865  static int SANC_MC_INIT = -123456789;
866 
867  //... Initialization of the W->l\nu\gamma
868  //... decay Matrix Element parameters
869  if (SANC_MC_INIT == -123456789) {
870  SANC_MC_INIT = 1;
871 
872  pi = 4 * atan(1.0);
873  qf1 = 0.0; // neutrino charge
874  mf1 = 1.0e-10; // newutrino mass
875  vf = 1.0; // V&A couplings
876  af = 1.0;
877  alphaI = 1.0 / ALPHA;
878  cw = 0.881731727; // Weak Weinberg angle
879  sw = 0.471751166;
880 
881 
882  //... An auxilary K&S vectors
883  bet[0] = 1.0;
884  bet[1] = 0.0722794881816159;
885  bet[2] = -0.994200045099866;
886  bet[3] = 0.0796363353729248;
887 
888  spV[0] = 0.0;
889  spV[1] = 7.22794881816159e-2;
890  spV[2] = -0.994200045099866;
891  spV[3] = 7.96363353729248e-2;
892 
893  mcLUN = PHLUN;
894  }
895  }
896 //----------------------------------------------------------------------
897 //
898 // PHOTOS: PHOtos Boson W correction weight
899 //
900 // Purpose: calculates correction weight due to amplitudes of
901 // emission from W boson. It is ecact, but not verified
902 // for exponentiation yet.
903 //
904 //
905 //
906 //
907 // Input Parameters: Common /PHOEVT/, with photon added.
908 // wt to be corrected
909 //
910 //
911 //
912 // Output Parameters: wt
913 //
914 // Author(s): G. Nanava, Z. Was Created at: 13/03/03
915 // Last Update: 22/06/13
916 //
917 //----------------------------------------------------------------------
918  void PhotosMEforW::PHOBWnlo(double* WT)
919  {
920  // FILE *PHLUN = stdout; // printouts from matrix element calculations
921  // directed with phlun still
922  int phlun = 6;
923  double EMU, MCHREN, BETA, COSTHG, MPASQR, XPH;
924  double PW[4], PMU[4], PPHOT[4], PNE[4];
925  double B_PW[4], B_PNE[4], B_PMU[4]; //,AMPSQR;
926  static int i = 1;
927  int I, IJ, I3, I4, JJ;
928  double MB, MF1, MF2, QB, QF2;
929  // double pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af;
930 
931 
935 
936  //--
937  if (abs(pho.idhep[1 - i]) == 24 &&
938  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i ]) >= 11 &&
939  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i ]) <= 16 &&
940  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i + 1]) >= 11 &&
941  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i + 1]) <= 16) {
942 
943  if (
944  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i ]) == 11 ||
945  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i ]) == 13 ||
946  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i ]) == 15) {
947  I = pho.jdahep[1 - i][1 - i];
948  } else {
949  I = pho.jdahep[1 - i][1 - i] + 1;
950  }
951  //.. muon energy
952  EMU = pho.phep[I - i][4 - i];
953  //.. muon mass square
954  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]
955  - 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  BETA = sqrt(1 - MCHREN / pho.phep[I - i][4 - i] * pho.phep[I - i][4 - i]);
957  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]
958  + pho.phep[I - i][1 - i] * pho.phep[pho.nhep - i][1 - i]) /
959  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] *
960  pho.phep[I - i][1 - i]) /
961  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] +
962  pho.phep[pho.nhep - i][1 - i] * pho.phep[pho.nhep - i][1 - i]));
963  MPASQR = pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i];
964  XPH = pho.phep[pho.nhep - i][4 - i];
965 
966  //... Initialization of the W->l\nu\gamma
967  //... decay Matrix Element parameters
968  SANC_INIT(phocop.alpha, phlun);
969 
970 
971  MB = pho.phep[1 - i][4 - i]; // ! W boson mass
972  MF2 = sqrt(MCHREN); // ! muon mass
973  I3 = -1;
974  for (IJ = 1; IJ <= hep.nhep; IJ++) {
975  if (abs(hep.idhep[IJ - i]) == 24) { I3 = IJ;}
976  }
977  if (I3 == -1) {cout << " ERROR IN PHOBWnlo of PHOTS W-ME: I3= &2i" << I3 << endl;}
978  if (
979  abs(hep.idhep[hep.jdahep[I3 - i][1 - i] - i ]) == 11 ||
980  abs(hep.idhep[hep.jdahep[I3 - i][1 - i] - i ]) == 13 ||
981  abs(hep.idhep[hep.jdahep[I3 - i][1 - i] - i ]) == 15) {
982  I4 = hep.jdahep[I3 - i][1 - i];
983  } // ! position of lepton
984  else {
985  I4 = hep.jdahep[I3 - i][1 - i] + 1 ; // ! position of lepton
986  }
987 
988  if (hep.idhep[I3 - i] == -24) QB = -1.0; // ! W boson charge
989  if (hep.idhep[I3 - i] == +24) QB = +1.0; //
990  if (hep.idhep[I4 - i] > 0.0) QF2 = -1.0; // ! lepton charge
991  if (hep.idhep[I4 - i] < 0.0) QF2 = +1.0;
992 
993 
994  //... Particle momenta before foton radiation; effective Born level
995  for (JJ = 1; JJ <= 4; JJ++) {
996  B_PW [(JJ % 4)] = hep.phep[I3 - i][JJ - i]; // ! W boson
997  B_PNE[(JJ % 4)] = hep.phep[I3 - i][JJ - i] - hep.phep[I4 - i][JJ - i]; // ! neutrino
998  B_PMU[(JJ % 4)] = hep.phep[I4 - i][JJ - i]; // ! muon
999  }
1000 
1001  //.. Particle monenta after photon radiation
1002  for (JJ = 1; JJ <= 4; JJ++) {
1003  PW [(JJ % 4)] = pho.phep[1 - i][JJ - i];
1004  PMU [(JJ % 4)] = pho.phep[I - i][JJ - i];
1005  PPHOT[(JJ % 4)] = pho.phep[pho.nhep - i][JJ - i];
1006  PNE [(JJ % 4)] = pho.phep[1 - i][JJ - i] - pho.phep[I - i][JJ - i] - pho.phep[pho.nhep - i][JJ - i];
1007  }
1008 
1009  // two options of calculating neutrino (spectator) mass
1010  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]));
1011  MF1 = sqrt(fabs(PNE[0] * PNE[0] - PNE[1] * PNE[1] - PNE[2] * PNE[2] - PNE[3] * PNE[3]));
1012 
1013  SANC_INIT1(QB, QF2, MF1, MF2, MB);
1014  *WT = (*WT) * SANC_WT(PW, PNE, PMU, PPHOT, B_PW, B_PNE, B_PMU);
1015  }
1016  // write(*,*) 'AMPSQR/EIKONALFACTOR= ', AMPSQR/EIKONALFACTOR
1017  }
1018 
1019 } // namespace Photospp
1020 
Photospp::PhotosMEforW::PHOBWnlo
static void PHOBWnlo(double *WT)
Definition: forW-MEc.cc:918
scripts.bitstring.equal
def equal(a, b)
Definition: bitstring.py:292