Bug Summary

File:generators/kkmc/photospp/forW-MEc.cc
Warning:line 957, column 7
Value stored to 'COSTHG' is never read

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -O3 -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name forW-MEc.cc -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/data/b2soft/buildbot/development/build -fcoverage-compilation-dir=/data/b2soft/buildbot/development/build -resource-dir /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++ -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/x86_64-redhat-linux -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/backward -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/python3.12 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/CLHEP -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/Geant4 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/root -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/belle_legacy -I include/ -D _PACKAGE_="generators" -D G4UI_USE_TCSH -D RaveDllExport= -D HAS_SQLITE -D HAS_CALLGRIND -I include -I /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/libxml2 -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++ -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/x86_64-redhat-linux -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/backward -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21/include -internal-isystem /usr/local/include -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../x86_64-redhat-linux/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -Wno-missing-braces -Wno-unused-command-line-argument -std=c++20 -fdeprecated-macro -ferror-limit 19 -fgnuc-version=4.2.1 -fno-implicit-modules -fskip-odr-check-in-gmf -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /scan_build/2026-05-31-004316-385593-1 -x c++ generators/kkmc/photospp/forW-MEc.cc
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// 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
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
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
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 {
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
120/////////////////////////////////////////////////////////////////////
121// //
122// this is small B_{s}(k,p) function when TrMartix is diaagonal!! //
123// //
124/////////////////////////////////////////////////////////////////////
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
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,
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
932 //! write(*,*) 'IDPHOs=',IDPHO(1),IDPHO(2),IDPHO(3),IDPHO(4),IDPHO(5)
933 //! write(*,*) 'IDPHOs=',pho.jdahep[1-i][1-i],npho
934 //! write(*,*) 'hep.IDPHOs=',hep.IDhep(1),hep.IDhep(2),hep.IDhep(3),hep.IDhep(4),hep.IDhep(5)
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]
Value stored to 'COSTHG' is never read
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;} //! position of W
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