]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8130/src/SigmaEW.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaEW.cxx
CommitLineData
5ad4eb21 1// SigmaEW.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for the
7// electroweak simulation classes (including photon processes).
8
9#include "SigmaEW.h"
10
11namespace Pythia8 {
12
13//**************************************************************************
14
15// Sigma2qg2qgamma class.
16// Cross section for q g -> q gamma.
17
18//*********
19
20// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
21
22void Sigma2qg2qgamma::sigmaKin() {
23
24 // Calculate kinematics dependence.
25 sigUS = (1./3.) * (sH2 + uH2) / (-sH * uH);
26
27 // Answer.
28 sigma0 = (M_PI/sH2) * alpS * alpEM * sigUS;
29
30 }
31
32//*********
33
34// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
35
36double Sigma2qg2qgamma::sigmaHat() {
37
38 // Incoming flavour gives charge factor.
39 int idNow = (id2 == 21) ? id1 : id2;
40 double eNow = CoupEW::ef( abs(idNow) );
41 return sigma0 * pow2(eNow);
42
43}
44
45//*********
46
47// Select identity, colour and anticolour.
48
49void Sigma2qg2qgamma::setIdColAcol() {
50
51 // Construct outgoing flavours.
52 id3 = (id1 == 21) ? 22 : id1;
53 id4 = (id2 == 21) ? 22 : id2;
54 setId( id1, id2, id3, id4);
55
56 // Colour flow topology. Swap if first is gluon, or when antiquark.
57 setColAcol( 1, 0, 2, 1, 2, 0, 0, 0);
58 if (id1 == 21) swapCol1234();
59 if (id1 < 0 || id2 < 0) swapColAcol();
60
61}
62
63//**************************************************************************
64
65// Sigma2qqbar2ggamma class.
66// Cross section for q qbar -> g gamma.
67
68//*********
69
70// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
71
72void Sigma2qqbar2ggamma::sigmaKin() {
73
74 // Calculate kinematics dependence.
75 double sigTU = (8./9.) * (tH2 + uH2) / (tH * uH);
76
77 // Answer.
78 sigma0 = (M_PI/sH2) * alpS * alpEM * sigTU;
79
80}
81
82//*********
83
84// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
85
86double Sigma2qqbar2ggamma::sigmaHat() {
87
88 // Incoming flavour gives charge factor.
89 double eNow = CoupEW::ef( abs(id1) );
90 return sigma0 * pow2(eNow);
91
92}
93
94//*********
95
96// Select identity, colour and anticolour.
97
98void Sigma2qqbar2ggamma::setIdColAcol() {
99
100 // Outgoing flavours trivial.
101 setId( id1, id2, 21, 22);
102
103 // One colour flow topology. Swap if first is antiquark.
104 setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
105 if (id1 < 0) swapColAcol();
106
107}
108
109//**************************************************************************
110
111// Sigma2gg2ggamma class.
112// Cross section for g g -> g gamma.
113// Proceeds through a quark box, by default using 5 massless quarks.
114
115//*********
116
117// Initialize process, especially parton-flux object.
118
119void Sigma2gg2ggamma::initProc() {
120
121 // Maximum quark flavour in loop.
122 int nQuarkLoop = Settings::mode("PromptPhoton:nQuarkLoop");
123
124 // Calculate charge factor from the allowed quarks in the box.
125 chargeSum = - 1./3. + 2./3. - 1./3.;
126 if (nQuarkLoop >= 4) chargeSum += 2./3.;
127 if (nQuarkLoop >= 5) chargeSum -= 1./3.;
128 if (nQuarkLoop >= 6) chargeSum += 2./3.;
129
130}
131
132//*********
133
134// Evaluate d(sigmaHat)/d(tHat).
135
136void Sigma2gg2ggamma::sigmaKin() {
137
138 // Logarithms of Mandelstam variable ratios.
139 double logST = log( -sH / tH );
140 double logSU = log( -sH / uH );
141 double logTU = log( tH / uH );
142
143 // Real and imaginary parts of separate amplitudes.
144 double b0stuRe = 1. + (tH - uH) / sH * logTU
145 + 0.5 * (tH2 + uH2) / sH2 * (pow2(logTU) + pow2(M_PI));
146 double b0stuIm = 0.;
147 double b0tsuRe = 1. + (sH - uH) / tH * logSU
148 + 0.5 * (sH2 + uH2) / tH2 * pow2(logSU);
149 double b0tsuIm = -M_PI * ( (sH - uH) / tH + (sH2 + uH2) / tH2 * logSU);
150 double b0utsRe = 1. + (sH - tH) / uH * logST
151 + 0.5 * (sH2 + tH2) / uH2 * pow2(logST);
152 double b0utsIm = -M_PI * ( (sH - tH) / uH + (sH2 + tH2) / uH2 * logST);
153 double b1stuRe = -1.;
154 double b1stuIm = 0.;
155 double b2stuRe = -1.;
156 double b2stuIm = 0.;
157
158 // Calculate kinematics dependence.
159 double sigBox = pow2(b0stuRe) + pow2(b0stuIm) + pow2(b0tsuRe)
160 + pow2(b0tsuIm) + pow2(b0utsRe) + pow2(b0utsIm) + 4. * pow2(b1stuRe)
161 + 4. * pow2(b1stuIm) + pow2(b2stuRe) + pow2(b2stuIm);
162
163 // Answer.
164 sigma = (5. / (192. * M_PI * sH2)) * pow2(chargeSum)
165 * pow3(alpS) * alpEM * sigBox;
166
167}
168
169//*********
170
171// Select identity, colour and anticolour.
172
173void Sigma2gg2ggamma::setIdColAcol() {
174
175 // Flavours and colours are trivial.
176 setId( id1, id2, 21, 22);
177 setColAcol( 1, 2, 2, 3, 1, 3, 0, 0);
178 if (Rndm::flat() > 0.5) swapColAcol();
179
180}
181
182//**************************************************************************
183
184// Sigma2ffbar2gammagamma class.
185// Cross section for q qbar -> gamma gamma.
186
187//*********
188
189// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
190
191void Sigma2ffbar2gammagamma::sigmaKin() {
192
193 // Calculate kinematics dependence.
194 sigTU = 2. * (tH2 + uH2) / (tH * uH);
195
196 // Answer contains factor 1/2 from identical photons.
197 sigma0 = (M_PI/sH2) * pow2(alpEM) * 0.5 * sigTU;
198
199}
200
201//*********
202
203// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
204
205double Sigma2ffbar2gammagamma::sigmaHat() {
206
207 // Incoming flavour gives charge and colour factors.
208 double eNow = CoupEW::ef( abs(id1) );
209 double colFac = (abs(id1) < 9) ? 1. / 3. : 1.;
210 return sigma0 * pow4(eNow) * colFac;
211
212}
213
214//*********
215
216// Select identity, colour and anticolour.
217
218void Sigma2ffbar2gammagamma::setIdColAcol() {
219
220 // Outgoing flavours trivial.
221 setId( id1, id2, 22, 22);
222
223 // No colours at all or one flow topology. Swap if first is antiquark.
224 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
225 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
226 if (id1 < 0) swapColAcol();
227
228}
229
230//**************************************************************************
231
232// Sigma2gg2gammagamma class.
233// Cross section for g g -> gamma gamma.
234// Proceeds through a quark box, by default using 5 massless quarks.
235
236//*********
237
238// Initialize process.
239
240void Sigma2gg2gammagamma::initProc() {
241
242 // Maximum quark flavour in loop.
243 int nQuarkLoop = Settings::mode("PromptPhoton:nQuarkLoop");
244
245 // Calculate charge factor from the allowed quarks in the box.
246 charge2Sum = 1./9. + 4./9. + 1./9.;
247 if (nQuarkLoop >= 4) charge2Sum += 4./9.;
248 if (nQuarkLoop >= 5) charge2Sum += 1./9.;
249 if (nQuarkLoop >= 6) charge2Sum += 4./9.;
250
251}
252
253//*********
254
255// Evaluate d(sigmaHat)/d(tHat).
256
257void Sigma2gg2gammagamma::sigmaKin() {
258
259 // Logarithms of Mandelstam variable ratios.
260 double logST = log( -sH / tH );
261 double logSU = log( -sH / uH );
262 double logTU = log( tH / uH );
263
264 // Real and imaginary parts of separate amplitudes.
265 double b0stuRe = 1. + (tH - uH) / sH * logTU
266 + 0.5 * (tH2 + uH2) / sH2 * (pow2(logTU) + pow2(M_PI));
267 double b0stuIm = 0.;
268 double b0tsuRe = 1. + (sH - uH) / tH * logSU
269 + 0.5 * (sH2 + uH2) / tH2 * pow2(logSU);
270 double b0tsuIm = -M_PI * ( (sH - uH) / tH + (sH2 + uH2) / tH2 * logSU);
271 double b0utsRe = 1. + (sH - tH) / uH * logST
272 + 0.5 * (sH2 + tH2) / uH2 * pow2(logST);
273 double b0utsIm = -M_PI * ( (sH - tH) / uH + (sH2 + tH2) / uH2 * logST);
274 double b1stuRe = -1.;
275 double b1stuIm = 0.;
276 double b2stuRe = -1.;
277 double b2stuIm = 0.;
278
279 // Calculate kinematics dependence.
280 double sigBox = pow2(b0stuRe) + pow2(b0stuIm) + pow2(b0tsuRe)
281 + pow2(b0tsuIm) + pow2(b0utsRe) + pow2(b0utsIm) + 4. * pow2(b1stuRe)
282 + 4. * pow2(b1stuIm) + pow2(b2stuRe) + pow2(b2stuIm);
283
284 // Answer contains factor 1/2 from identical photons.
285 sigma = (0.5 / (16. * M_PI * sH2)) * pow2(charge2Sum)
286 * pow2(alpS) * pow2(alpEM) * sigBox;
287
288}
289
290//*********
291
292// Select identity, colour and anticolour.
293
294void Sigma2gg2gammagamma::setIdColAcol() {
295
296 // Flavours and colours are trivial.
297 setId( id1, id2, 22, 22);
298 setColAcol( 1, 2, 2, 1, 0, 0, 0, 0);
299
300}
301
302//**************************************************************************
303
304// Sigma2ff2fftgmZ class.
305// Cross section for f f' -> f f' via t-channel gamma*/Z0 exchange
306// (f is quark or lepton).
307
308//*********
309
310// Initialize process.
311
312void Sigma2ff2fftgmZ::initProc() {
313
314 // Store Z0 mass for propagator. Common coupling factor.
315 gmZmode = Settings::mode("WeakZ0:gmZmode");
316 mZ = ParticleDataTable::m0(23);
317 mZS = mZ*mZ;
318 thetaWRat = 1. / (16. * CoupEW::sin2thetaW() * CoupEW::cos2thetaW());
319
320}
321
322//*********
323
324// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
325
326void Sigma2ff2fftgmZ::sigmaKin() {
327
328 // Cross section part common for all incoming flavours.
329 double sigma0 = (M_PI / sH2) * pow2(alpEM);
330
331 // Kinematical functions for gamma-gamma, gamma-Z and Z-Z parts.
332 sigmagmgm = sigma0 * 2. * (sH2 + uH2) / tH2;
333 sigmagmZ = sigma0 * 4. * thetaWRat * sH2 / (tH * (tH - mZS));
334 sigmaZZ = sigma0 * 2. * pow2(thetaWRat) * sH2 / pow2(tH - mZS);
335 if (gmZmode == 1) {sigmagmZ = 0.; sigmaZZ = 0.;}
336 if (gmZmode == 2) {sigmagmgm = 0.; sigmagmZ = 0.;}
337
338}
339
340//*********
341
342// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
343
344double Sigma2ff2fftgmZ::sigmaHat() {
345
346 // Couplings for current flavour combination.
347 int id1Abs = abs(id1);
348 double e1 = CoupEW::ef(id1Abs);
349 double v1 = CoupEW::vf(id1Abs);
350 double a1 = CoupEW::af(id1Abs);
351 int id2Abs = abs(id2);
352 double e2 = CoupEW::ef(id2Abs);
353 double v2 = CoupEW::vf(id2Abs);
354 double a2 = CoupEW::af(id2Abs);
355
356 // Distinguish same-sign and opposite-sign fermions.
357 double epsi = (id1 * id2 > 0) ? 1. : -1.;
358
359 // Flavour-dependent cross section.
360 double sigma = sigmagmgm * pow2(e1 * e2)
361 + sigmagmZ * e1 * e2 * (v1 * v2 * (1. + uH2 / sH2)
362 + a1 * a2 * epsi * (1. - uH2 / sH2))
363 + sigmaZZ * ((v1*v1 + a1*a1) * (v2*v2 + a2*a2) * (1. + uH2 / sH2)
364 + 4. * v1 * a1 * v2 * a2 * epsi * (1. - uH2 / sH2));
365
366 // Spin-state extra factor 2 per incoming neutrino.
367 if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
368 if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
369
370 // Answer.
371 return sigma;
372
373}
374
375//*********
376
377// Select identity, colour and anticolour.
378
379void Sigma2ff2fftgmZ::setIdColAcol() {
380
381 // Trivial flavours: out = in.
382 setId( id1, id2, id1, id2);
383
384 // Colour flow topologies. Swap when antiquarks.
385 if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
386 setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
387 else if (abs(id1) < 9 && abs(id2) < 9)
388 setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
389 else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
390 else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
391 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
392 if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
393 swapColAcol();
394
395}
396
397//**************************************************************************
398
399// Sigma2ff2fftW class.
400// Cross section for f_1 f_2 -> f_3 f_4 via t-channel W+- exchange
401// (f is quark or lepton).
402
403//*********
404
405// Initialize process.
406
407void Sigma2ff2fftW::initProc() {
408
409 // Store W+- mass for propagator. Common coupling factor.
410 mW = ParticleDataTable::m0(24);
411 mWS = mW*mW;
412 thetaWRat = 1. / (4. * CoupEW::sin2thetaW());
413
414}
415
416//*********
417
418// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
419
420void Sigma2ff2fftW::sigmaKin() {
421
422 // Cross section part common for all incoming flavours.
423 sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat)
424 * 4. * sH2 / pow2(tH - mWS);
425
426}
427
428//*********
429
430// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
431
432double Sigma2ff2fftW::sigmaHat() {
433
434 // Some flavour combinations not possible.
435 int id1Abs = abs(id1);
436 int id2Abs = abs(id2);
437 if ( (id1Abs%2 == id2Abs%2 && id1 * id2 > 0)
438 || (id1Abs%2 != id2Abs%2 && id1 * id2 < 0) ) return 0.;
439
440 // Basic cross section.
441 double sigma = sigma0;
442 if (id1 * id2 < 0) sigma *= uH2 / sH2;
443
444 // CKM factors for final states.
445 sigma *= VCKM::V2sum(id1Abs) * VCKM::V2sum(id2Abs);
446
447 // Spin-state extra factor 2 per incoming neutrino.
448 if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
449 if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
450
451 // Answer.
452 return sigma;
453
454}
455
456//*********
457
458// Select identity, colour and anticolour.
459
460void Sigma2ff2fftW::setIdColAcol() {
461
462 // Pick out-flavours by relative CKM weights.
463 id3 = VCKM::V2pick(id1);
464 id4 = VCKM::V2pick(id2);
465 setId( id1, id2, id3, id4);
466
467 // Colour flow topologies. Swap when antiquarks.
468 if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
469 setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
470 else if (abs(id1) < 9 && abs(id2) < 9)
471 setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
472 else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
473 else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
474 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
475 if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
476 swapColAcol();
477
478}
479
480
481//**************************************************************************
482
483// Sigma2qq2QqtW class.
484// Cross section for q q' -> Q q" via t-channel W+- exchange.
485// Related to Sigma2ff2ffViaW class, but with massive matrix elements.
486
487//*********
488
489// Initialize process.
490
491void Sigma2qq2QqtW::initProc() {
492
493 // Process name.
494 nameSave = "q q -> Q q (t-channel W+-)";
495 if (idNew == 4) nameSave = "q q -> c q (t-channel W+-)";
496 if (idNew == 5) nameSave = "q q -> b q (t-channel W+-)";
497 if (idNew == 6) nameSave = "q q -> t q (t-channel W+-)";
498 if (idNew == 7) nameSave = "q q -> b' q (t-channel W+-)";
499 if (idNew == 8) nameSave = "q q -> t' q (t-channel W+-)";
500
501 // Store W+- mass for propagator. Common coupling factor.
502 mW = ParticleDataTable::m0(24);
503 mWS = mW*mW;
504 thetaWRat = 1. / (4. * CoupEW::sin2thetaW());
505
506 // Secondary open width fractions, relevant for top (or heavier).
507 openFracPos = ParticleDataTable::resOpenFrac(idNew);
508 openFracNeg = ParticleDataTable::resOpenFrac(-idNew);
509
510}
511
512//*********
513
514// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
515
516void Sigma2qq2QqtW::sigmaKin() {
517
518 // Cross section part common for all incoming flavours.
519 sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat) * 4. / pow2(tH - mWS);
520
521}
522
523//*********
524
525// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
526
527double Sigma2qq2QqtW::sigmaHat() {
528
529 // Some flavour combinations not possible.
530 int id1Abs = abs(id1);
531 int id2Abs = abs(id2);
532 bool diff12 = (id1Abs%2 != id2Abs%2);
533 if ( (!diff12 && id1 * id2 > 0)
534 || ( diff12 && id1 * id2 < 0) ) return 0.;
535
536 // Basic cross section.
537 double sigma = sigma0;
538 sigma *= (id1 * id2 > 0) ? sH * (sH - s3) : uH * (uH - s3);
539
540 // Secondary width if t or tbar produced on either side.
541 double openFrac1 = (id1 > 0) ? openFracPos : openFracNeg;
542 double openFrac2 = (id2 > 0) ? openFracPos : openFracNeg;
543
544 // CKM factors for final states; further impossible case.
545 bool diff1N = (id1Abs%2 != idNew%2);
546 bool diff2N = (id2Abs%2 != idNew%2);
547 if (diff1N && diff2N)
548 sigma *= ( VCKM::V2id(id1Abs, idNew) * openFrac1 * VCKM::V2sum(id2Abs)
549 + VCKM::V2sum(id1Abs) * VCKM::V2id(id2Abs, idNew) * openFrac2 );
550 else if (diff1N)
551 sigma *= VCKM::V2id(id1Abs, idNew) * openFrac1 * VCKM::V2sum(id2Abs);
552 else if (diff2N)
553 sigma *= VCKM::V2sum(id1Abs) * VCKM::V2id(id2Abs, idNew) * openFrac2;
554 else sigma = 0.;
555
556 // Spin-state extra factor 2 per incoming neutrino.
557 if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
558 if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
559
560 // Answer.
561 return sigma;
562
563}
564
565//*********
566
567// Select identity, colour and anticolour.
568
569void Sigma2qq2QqtW::setIdColAcol() {
570
571 // For topologies like d dbar -> (t/c/u) (t/c/u)bar pick side.
572 int id1Abs = abs(id1);
573 int id2Abs = abs(id2);
574 int side = 1;
575 if ( (id1Abs + idNew)%2 == 1 && (id2Abs + idNew)%2 == 1 ) {
576 double prob1 = VCKM::V2id(id1Abs, idNew) * VCKM::V2sum(id2Abs);
577 prob1 *= (id1 > 0) ? openFracPos : openFracNeg;
578 double prob2 = VCKM::V2id(id2Abs, idNew) * VCKM::V2sum(id1Abs);
579 prob2 *= (id2 > 0) ? openFracPos : openFracNeg;
580 if (prob2 > Rndm::flat() * (prob1 + prob2)) side = 2;
581 }
582 else if ((id2Abs + idNew)%2 == 1) side = 2;
583
584 // Pick out-flavours by relative CKM weights.
585 if (side == 1) {
586 // q q' -> t q" : correct order from start.
587 id3 = (id1 > 0) ? idNew : -idNew;
588 id4 = VCKM::V2pick(id2);
589 setId( id1, id2, id3, id4);
590 } else {
591 // q q' -> q" t : stored as t q" so swap tHat <-> uHat.
592 swapTU = true;
593 id3 = VCKM::V2pick(id1);
594 id4 = (id2 > 0) ? idNew : -idNew;
595 setId( id1, id2, id4, id3);
596 }
597
598 // Colour flow topologies. Swap when antiquarks on side 1.
599 if (side == 1 && id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
600 else if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
601 else if (side == 1) setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
602 else setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
603 if (id1 < 0) swapColAcol();
604
605}
606
607//*********
608
609// Evaluate weight for decay angles of W in top decay.
610
611double Sigma2qq2QqtW::weightDecay( Event& process, int iResBeg,
612 int iResEnd) {
613
614 // For top decay hand over to standard routine, else done.
615 if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
616 return weightTopDecay( process, iResBeg, iResEnd);
617 else return 1.;
618
619}
620
621//**************************************************************************
622
623// Sigma1ffbar2gmZ class.
624// Cross section for f fbar -> gamma*/Z0 (f is quark or lepton).
625
626//*********
627
628// Initialize process.
629
630void Sigma1ffbar2gmZ::initProc() {
631
632 // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
633 gmZmode = Settings::mode("WeakZ0:gmZmode");
634
635 // Store Z0 mass and width for propagator.
636 mRes = ParticleDataTable::m0(23);
637 GammaRes = ParticleDataTable::mWidth(23);
638 m2Res = mRes*mRes;
639 GamMRat = GammaRes / mRes;
640 thetaWRat = 1. / (16. * CoupEW::sin2thetaW() * CoupEW::cos2thetaW());
641
642 // Set pointer to particle properties and decay table.
643 particlePtr = ParticleDataTable::particleDataPtr(23);
644
645}
646
647//*********
648
649// Evaluate sigmaHat(sHat), part independent of incoming flavour.
650
651void Sigma1ffbar2gmZ::sigmaKin() {
652
653 // Common coupling factors.
654 double colQ = 3. * (1. + alpS / M_PI);
655
656 // Reset quantities to sum. Declare variables in loop.
657 gamSum = 0.;
658 intSum = 0.;
659 resSum = 0.;
660 int idAbs, onMode;
661 double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
662
663 // Loop over all Z0 decay channels.
664 for (int i = 0; i < particlePtr->decay.size(); ++i) {
665 idAbs = abs( particlePtr->decay[i].product(0) );
666
667 // Only contributions from three fermion generations, except top.
668 if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
669 mf = ParticleDataTable::m0(idAbs);
670
671 // Check that above threshold. Phase space.
672 if (mH > 2. * mf + MASSMARGIN) {
673 mr = pow2(mf / mH);
674 betaf = sqrtpos(1. - 4. * mr);
675 psvec = betaf * (1. + 2. * mr);
676 psaxi = pow3(betaf);
677
678 // Combine phase space with couplings.
679 ef2 = CoupEW::ef2(idAbs) * psvec;
680 efvf = CoupEW::efvf(idAbs) * psvec;
681 vf2af2 = CoupEW::vf2(idAbs) * psvec + CoupEW::af2(idAbs) * psaxi;
682 colf = (idAbs < 6) ? colQ : 1.;
683
684 // Store sum of combinations. For outstate only open channels.
685 onMode = particlePtr->decay[i].onMode();
686 if (onMode == 1 || onMode == 2) {
687 gamSum += colf * ef2;
688 intSum += colf * efvf;
689 resSum += colf * vf2af2;
690 }
691
692 // End loop over fermions.
693 }
694 }
695 }
696
697 // Calculate prefactors for gamma/interference/Z0 cross section terms.
698 gamProp = 4. * M_PI * pow2(alpEM) / (3. * sH);
699 intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
700 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
701 resProp = gamProp * pow2(thetaWRat * sH)
702 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
703
704 // Optionally only keep gamma* or Z0 term.
705 if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
706 if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
707
708}
709
710//*********
711
712// Evaluate sigmaHat(sHat), including incoming flavour dependence.
713
714double Sigma1ffbar2gmZ::sigmaHat() {
715
716 // Combine gamma, interference and Z0 parts.
717 int idAbs = abs(id1);
718 double sigma = CoupEW::ef2(idAbs) * gamProp * gamSum
719 + CoupEW::efvf(idAbs) * intProp * intSum
720 + CoupEW::vf2af2(idAbs) * resProp * resSum;
721
722 // Colour factor. Answer.
723 if (idAbs < 9) sigma /= 3.;
724 return sigma;
725
726}
727
728//*********
729
730// Select identity, colour and anticolour.
731
732void Sigma1ffbar2gmZ::setIdColAcol() {
733
734 // Flavours trivial.
735 setId( id1, id2, 23);
736
737 // Colour flow topologies. Swap when antiquarks.
738 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
739 else setColAcol( 0, 0, 0, 0, 0, 0);
740 if (id1 < 0) swapColAcol();
741
742}
743
744//*********
745
746// Evaluate weight for gamma*/Z0 decay angle.
747
748double Sigma1ffbar2gmZ::weightDecay( Event& process, int iResBeg,
749 int iResEnd) {
750
751 // Z should sit in entry 5.
752 if (iResBeg != 5 || iResEnd != 5) return 1.;
753
754 // Couplings for in- and out-flavours.
755 int idInAbs = process[3].idAbs();
756 double ei = CoupEW::ef(idInAbs);
757 double vi = CoupEW::vf(idInAbs);
758 double ai = CoupEW::af(idInAbs);
759 int idOutAbs = process[6].idAbs();
760 double ef = CoupEW::ef(idOutAbs);
761 double vf = CoupEW::vf(idOutAbs);
762 double af = CoupEW::af(idOutAbs);
763
764 // Phase space factors. (One power of beta left out in formulae.)
765 double mf = process[6].m();
766 double mr = mf*mf / sH;
767 double betaf = sqrtpos(1. - 4. * mr);
768
769 // Coefficients of angular expression.
770 double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
771 + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
772 double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
773 + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
774 double coefAsym = betaf * ( ei * ai * intProp * ef * af
775 + 4. * vi * ai * resProp * vf * af );
776
777 // Flip asymmetry for in-fermion + out-antifermion.
778 if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
779
780 // Reconstruct decay angle and weight for it.
781 double cosThe = (process[3].p() - process[4].p())
782 * (process[7].p() - process[6].p()) / (sH * betaf);
783 double wtMax = 2. * (coefTran + abs(coefAsym));
784 double wt = coefTran * (1. + pow2(cosThe))
785 + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
786
787 // Done.
788 return (wt / wtMax);
789
790}
791
792//**************************************************************************
793
794// Sigma1ffbar2W class.
795// Cross section for f fbar' -> W+- (f is quark or lepton).
796
797//*********
798
799// Initialize process.
800
801void Sigma1ffbar2W::initProc() {
802
803 // Store W+- mass and width for propagator.
804 mRes = ParticleDataTable::m0(24);
805 GammaRes = ParticleDataTable::mWidth(24);
806 m2Res = mRes*mRes;
807 GamMRat = GammaRes / mRes;
808 thetaWRat = 1. / (12. * CoupEW::sin2thetaW());
809
810 // Set pointer to particle properties and decay table.
811 particlePtr = ParticleDataTable::particleDataPtr(24);
812
813}
814
815//*********
816
817// Evaluate sigmaHat(sHat), part independent of incoming flavour.
818
819void Sigma1ffbar2W::sigmaKin() {
820
821 // Set up Breit-Wigner. Cross section for W+ and W- separately.
822 double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
823 double preFac = alpEM * thetaWRat * mH;
824 sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(34, mH);
825 sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-34, mH);
826
827}
828
829//*********
830
831// Evaluate sigmaHat(sHat), including incoming flavour dependence.
832
833double Sigma1ffbar2W::sigmaHat() {
834
835 // Secondary width for W+ or W-. CKM and colour factors.
836 int idUp = (abs(id1)%2 == 0) ? id1 : id2;
837 double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
838 if (abs(id1) < 9) sigma *= VCKM::V2id(abs(id1), abs(id2)) / 3.;
839
840 // Answer.
841 return sigma;
842
843}
844
845//*********
846
847// Select identity, colour and anticolour.
848
849void Sigma1ffbar2W::setIdColAcol() {
850
851 // Sign of outgoing W.
852 int sign = 1 - 2 * (abs(id1)%2);
853 if (id1 < 0) sign = -sign;
854 setId( id1, id2, 24 * sign);
855
856 // Colour flow topologies. Swap when antiquarks.
857 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
858 else setColAcol( 0, 0, 0, 0, 0, 0);
859 if (id1 < 0) swapColAcol();
860
861}
862
863//*********
864
865// Evaluate weight for W decay angle.
866
867double Sigma1ffbar2W::weightDecay( Event& process, int iResBeg,
868 int iResEnd) {
869
870 // W should sit in entry 5.
871 if (iResBeg != 5 || iResEnd != 5) return 1.;
872
873 // Phase space factors.
874 double mr1 = pow2(process[6].m()) / sH;
875 double mr2 = pow2(process[7].m()) / sH;
876 double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
877
878 // Sign of asymmetry.
879 double eps = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
880
881 // Reconstruct decay angle and weight for it.
882 double cosThe = (process[3].p() - process[4].p())
883 * (process[7].p() - process[6].p()) / (sH * betaf);
884 double wtMax = 4.;
885 double wt = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
886
887 // Done.
888 return (wt / wtMax);
889
890}
891
892//**************************************************************************
893
894// Sigma2ffbar2ffbarsgm class.
895// Cross section f fbar -> gamma* -> f' fbar', for multiple interactions.
896
897
898//*********
899
900// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
901
902void Sigma2ffbar2ffbarsgm::sigmaKin() {
903
904 // Pick new flavour. Allow three leptons and five quarks.
905 double colQ = 1. + (alpS / M_PI);
906 double flavWt = 3. + colQ * 11. / 3.;
907 double flavRndm = Rndm::flat() * flavWt;
908 if (flavRndm < 3.) {
909 if (flavRndm < 1.) idNew = 11;
910 else if (flavRndm < 2.) idNew = 13;
911 else idNew = 15;
912 } else {
913 flavRndm = 3. * (flavWt - 3.) / colQ;
914 if (flavRndm < 4.) idNew = 2;
915 else if (flavRndm < 8.) idNew = 4;
916 else if (flavRndm < 9.) idNew = 1;
917 else if (flavRndm < 10.) idNew = 3;
918 else idNew = 5;
919 }
920 double mNew = ParticleDataTable::m0(idNew);
921 double m2New = mNew*mNew;
922
923 // Calculate kinematics dependence. Give correct mass factors for
924 // tHat, uHat defined as if massless kinematics, d(sigma)/d(Omega)
925 // = beta (1 + cos^2(theta) + (1 - beta^2) sin^2(theta)).
926 // Special case related to phase space form in multiple interactions.
927 double sigS = 0.;
928 if (sH > 4. * m2New) {
929 double beta = sqrt(1. - 4. * m2New / sH);
930 sigS = beta * (2.* (tH2 + uH2) + 4. * (1. - beta * beta) * tH * uH)
931 / sH2;
932 }
933
934 // Answer is proportional to number of outgoing flavours.
935 sigma0 = (M_PI/sH2) * pow2(alpEM) * sigS * flavWt;
936
937}
938
939//*********
940
941// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
942
943double Sigma2ffbar2ffbarsgm::sigmaHat() {
944
945 // Charge and colour factors.
946 double eNow = CoupEW::ef( abs(id1) );
947 double sigma = sigma0 * pow2(eNow);
948 if (abs(id1) < 9) sigma /= 3.;
949
950 // Answer.
951 return sigma;
952
953}
954
955//*********
956
957// Select identity, colour and anticolour.
958
959void Sigma2ffbar2ffbarsgm::setIdColAcol() {
960
961 // Set outgoing flavours.
962 id3 = (id1 > 0) ? idNew : -idNew;
963 setId( id1, id2, id3, -id3);
964
965 // Colour flow topologies. Swap when antiquarks.
966 if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
967 else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
968 else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
969 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
970 if (id1 < 0) swapColAcol();
971
972}
973
974//**************************************************************************
975
976// Sigma2ffbar2FFbarsgmZ class.
977// Cross section f fbar -> gamma*/Z0 -> F Fbar.
978
979//*********
980
981// Initialize process.
982
983void Sigma2ffbar2FFbarsgmZ::initProc() {
984
985 // Process name.
986 nameSave = "f fbar -> F Fbar (s-channel gamma*/Z0)";
987 if (idNew == 4) nameSave = "f fbar -> c cbar (s-channel gamma*/Z0)";
988 if (idNew == 5) nameSave = "f fbar -> b bbar (s-channel gamma*/Z0)";
989 if (idNew == 6) nameSave = "f fbar -> t tbar (s-channel gamma*/Z0)";
990 if (idNew == 7) nameSave = "f fbar -> b' b'bar (s-channel gamma*/Z0)";
991 if (idNew == 8) nameSave = "f fbar -> t' t'bar (s-channel gamma*/Z0)";
992 if (idNew == 15) nameSave = "f fbar -> tau+ tau- (s-channel gamma*/Z0)";
993 if (idNew == 17) nameSave = "f fbar -> tau'+ tau'- (s-channel gamma*/Z0)";
994 if (idNew == 18)
995 nameSave = "f fbar -> nu'_tau nu'bar_tau (s-channel gamma*/Z0)";
996
997 // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
998 gmZmode = Settings::mode("WeakZ0:gmZmode");
999
1000 // Store Z0 mass and width for propagator.
1001 mRes = ParticleDataTable::m0(23);
1002 GammaRes = ParticleDataTable::mWidth(23);
1003 m2Res = mRes*mRes;
1004 GamMRat = GammaRes / mRes;
1005 thetaWRat = 1. / (16. * CoupEW::sin2thetaW() * CoupEW::cos2thetaW());
1006
1007 // Store couplings of F.
1008 ef = CoupEW::ef(idNew);
1009 vf = CoupEW::vf(idNew);
1010 af = CoupEW::af(idNew);
1011
1012 // Secondary open width fraction, relevant for top (or heavier).
1013 openFracPair = ParticleDataTable::resOpenFrac(idNew, -idNew);
1014
1015}
1016
1017//*********
1018
1019// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1020
1021void Sigma2ffbar2FFbarsgmZ::sigmaKin() {
1022
1023 // Check that above threshold.
1024 isPhysical = true;
1025 if (mH < m3 + m4 + MASSMARGIN) {
1026 isPhysical = false;
1027 return;
1028 }
1029
1030 // Define average F, Fbar mass so same beta. Phase space.
1031 double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1032 mr = s34Avg / sH;
1033 betaf = sqrtpos(1. - 4. * mr);
1034
1035 // Final-state colour factor.
1036 double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) : 1.;
1037
1038 // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1039 cosThe = (tH - uH) / (betaf * sH);
1040
1041 // Calculate prefactors for gamma/interference/Z0 cross section terms.
1042 gamProp = colF * M_PI * pow2(alpEM) / sH2;
1043 intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1044 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1045 resProp = gamProp * pow2(thetaWRat * sH)
1046 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1047
1048 // Optionally only keep gamma* or Z0 term.
1049 if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1050 if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1051
1052}
1053
1054//*********
1055
1056// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1057
1058double Sigma2ffbar2FFbarsgmZ::sigmaHat() {
1059
1060 // Fail if below threshold.
1061 if (!isPhysical) return 0.;
1062
1063 // Couplings for in-flavours.
1064 int idAbs = abs(id1);
1065 double ei = CoupEW::ef(idAbs);
1066 double vi = CoupEW::vf(idAbs);
1067 double ai = CoupEW::af(idAbs);
1068
1069 // Coefficients of angular expression.
1070 double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
1071 + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
1072 double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
1073 + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
1074 double coefAsym = betaf * ( ei * ai * intProp * ef * af
1075 + 4. * vi * ai * resProp * vf * af );
1076
1077 // Combine gamma, interference and Z0 parts.
1078 double sigma = coefTran * (1. + pow2(cosThe))
1079 + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
1080
1081 // Top: corrections for closed decay channels.
1082 sigma *= openFracPair;
1083
1084 // Initial-state colour factor. Answer.
1085 if (idAbs < 9) sigma /= 3.;
1086 return sigma;
1087
1088}
1089
1090//*********
1091
1092// Select identity, colour and anticolour.
1093
1094void Sigma2ffbar2FFbarsgmZ::setIdColAcol() {
1095
1096 // Set outgoing flavours.
1097 id3 = (id1 > 0) ? idNew : -idNew;
1098 setId( id1, id2, id3, -id3);
1099
1100 // Colour flow topologies. Swap when antiquarks.
1101 if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1102 else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1103 else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1104 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1105 if (id1 < 0) swapColAcol();
1106
1107}
1108
1109//*********
1110
1111// Evaluate weight for decay angles of W in top decay.
1112
1113double Sigma2ffbar2FFbarsgmZ::weightDecay( Event& process, int iResBeg,
1114 int iResEnd) {
1115
1116 // For top decay hand over to standard routine, else done.
1117 if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1118 return weightTopDecay( process, iResBeg, iResEnd);
1119 else return 1.;
1120
1121}
1122
1123//**************************************************************************
1124
1125// Sigma2ffbar2FfbarsW class.
1126// Cross section f fbar' -> W+- -> F fbar".
1127
1128//*********
1129
1130// Initialize process.
1131
1132void Sigma2ffbar2FfbarsW::initProc() {
1133
1134 // Process name.
1135 nameSave = "f fbar -> F fbar (s-channel W+-)";
1136 if (idNew == 4) nameSave = "f fbar -> c qbar (s-channel W+-)";
1137 if (idNew == 5) nameSave = "f fbar -> b qbar (s-channel W+-)";
1138 if (idNew == 6) nameSave = "f fbar -> t qbar (s-channel W+-)";
1139 if (idNew == 7) nameSave = "f fbar -> b' qbar (s-channel W+-)";
1140 if (idNew == 8) nameSave = "f fbar -> t' qbar (s-channel W+-)";
1141 if (idNew == 7 && idNew2 == 6)
1142 nameSave = "f fbar -> b' tbar (s-channel W+-)";
1143 if (idNew == 8 && idNew2 == 7)
1144 nameSave = "f fbar -> t' b'bar (s-channel W+-)";
1145 if (idNew == 15 || idNew == 16)
1146 nameSave = "f fbar -> tau nu_taubar (s-channel W+-)";
1147 if (idNew == 17 || idNew == 18)
1148 nameSave = "f fbar -> tau' nu'_taubar (s-channel W+-)";
1149
1150 // Store W+- mass and width for propagator.
1151 mRes = ParticleDataTable::m0(24);
1152 GammaRes = ParticleDataTable::mWidth(24);
1153 m2Res = mRes*mRes;
1154 GamMRat = GammaRes / mRes;
1155 thetaWRat = 1. / (12. * CoupEW::sin2thetaW());
1156
1157 // For t/t' want to use at least b mass.
1158 idPartner = idNew2;
1159 if ( (idNew == 6 || idNew == 8) && idNew2 == 0 ) idPartner = 5;
1160
1161 // Sum of CKM weights for quarks.
1162 V2New = (idNew < 9) ? VCKM::V2sum(idNew) : 1.;
1163 if (idNew2 != 0) V2New = VCKM::V2id(idNew, idNew2);
1164
1165 // Secondary open width fractions, relevant for top or heavier.
1166 openFracPos = ParticleDataTable::resOpenFrac( idNew, -idNew2);
1167 openFracNeg = ParticleDataTable::resOpenFrac(-idNew, idNew2);
1168
1169}
1170
1171//*********
1172
1173// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1174
1175void Sigma2ffbar2FfbarsW::sigmaKin() {
1176
1177 // Check that above threshold.
1178 isPhysical = true;
1179 if (mH < m3 + m4 + MASSMARGIN) {
1180 isPhysical = false;
1181 return;
1182 }
1183
1184 // Phase space factors.
1185 double mr1 = s3 / sH;
1186 double mr2 = s4 / sH;
1187 double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
1188
1189 // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1190 double cosThe = (tH - uH) / (betaf * sH);
1191
1192 // Set up Breit-Wigner and in- and out-widths.
1193 double sigBW = 9. * M_PI * pow2(alpEM * thetaWRat)
1194 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1195
1196 // Initial-state colour factor.
1197 double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) * V2New : 1.;
1198
1199 // Angular dependence.
1200 double wt = pow2(1. + betaf * cosThe) - pow2(mr1 - mr2);
1201
1202 // Temporary answer.
1203 sigma0 = sigBW * colF * wt;
1204
1205}
1206
1207//*********
1208
1209// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1210
1211double Sigma2ffbar2FfbarsW::sigmaHat() {
1212
1213 // Fail if below threshold.
1214 if (!isPhysical) return 0.;
1215
1216 // CKM and colour factors.
1217 double sigma = sigma0;
1218 if (abs(id1) < 9) sigma *= VCKM::V2id(abs(id1), abs(id2)) / 3.;
1219
1220 // Correction for secondary width in top (or heavier) decay.
1221 int idSame = ((abs(id1) + idNew)%2 == 0) ? id1 : id2;
1222 sigma *= (idSame > 0) ? openFracPos : openFracNeg;
1223
1224 // Answer.
1225 return sigma;
1226
1227}
1228
1229//*********
1230
1231// Select identity, colour and anticolour.
1232
1233void Sigma2ffbar2FfbarsW::setIdColAcol() {
1234
1235 // Set outgoing flavours.
1236 id3 = idNew;
1237 id4 = (idNew2 != 0) ? idNew2 : VCKM::V2pick(idNew);
1238 if (idNew%2 == 0) {
1239 int idInUp = (abs(id1)%2 == 0) ? id1 : id2;
1240 if (idInUp > 0) id4 = -id4;
1241 else id3 = -id3;
1242 } else {
1243 int idInDn = (abs(id1)%2 == 1) ? id1 : id2;
1244 if (idInDn > 0) id4 = -id4;
1245 else id3 = -id3;
1246 }
1247 setId( id1, id2, id3, id4);
1248
1249 // Swap tHat and uHat for fbar' f -> F f".
1250 if (id1 * id3 < 0) swapTU = true;
1251
1252 // Colour flow topologies. Swap when antiquarks.
1253 if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1254 else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1255 else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1256 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1257 if (id1 < 0) swapCol12();
1258 if (id3 < 0) swapCol34();
1259
1260}
1261
1262//*********
1263
1264// Evaluate weight for decay angles of W in top decay.
1265
1266double Sigma2ffbar2FfbarsW::weightDecay( Event& process, int iResBeg,
1267 int iResEnd) {
1268
1269 // For top decay hand over to standard routine, else done.
1270 if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1271 return weightTopDecay( process, iResBeg, iResEnd);
1272 else return 1.;
1273
1274}
1275
1276//**************************************************************************
1277
1278// Sigma2ffbargmZWgmZW class.
1279// Collects common methods for f fbar -> gamma*/Z0/W+- gamma*/Z0/W-+.
1280
1281//*********
1282
1283// Calculate and store internal products.
1284
1285void Sigma2ffbargmZWgmZW::setupProd( Event& process, int i1, int i2,
1286 int i3, int i4, int i5, int i6) {
1287
1288 // Store incoming and outgoing momenta,
1289 pRot[1] = process[i1].p();
1290 pRot[2] = process[i2].p();
1291 pRot[3] = process[i3].p();
1292 pRot[4] = process[i4].p();
1293 pRot[5] = process[i5].p();
1294 pRot[6] = process[i6].p();
1295
1296 // Do random rotation to avoid accidental zeroes in HA expressions.
1297 bool smallPT = false;
1298 do {
1299 smallPT = false;
1300 double thetaNow = acos(2. * Rndm::flat() - 1.);
1301 double phiNow = 2. * M_PI * Rndm::flat();
1302 for (int i = 1; i <= 6; ++i) {
1303 pRot[i].rot( thetaNow, phiNow);
1304 if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
1305 }
1306 } while (smallPT);
1307
1308 // Calculate internal products.
1309 for (int i = 1; i < 6; ++i) {
1310 for (int j = i + 1; j <= 6; ++j) {
1311 hA[i][j] =
1312 sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
1313 / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
1314 - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
1315 / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
1316 hC[i][j] = conj( hA[i][j] );
1317 if (i <= 2) {
1318 hA[i][j] *= complex( 0., 1.);
1319 hC[i][j] *= complex( 0., 1.);
1320 }
1321 hA[j][i] = - hA[i][j];
1322 hC[j][i] = - hC[i][j];
1323 }
1324 }
1325
1326}
1327
1328//*********
1329
1330// Evaluate the F function of Gunion and Kunszt.
1331
1332complex Sigma2ffbargmZWgmZW::fGK(int j1, int j2, int j3, int j4, int j5,
1333 int j6) {
1334
1335 return 4. * hA[j1][j3] * hC[j2][j6]
1336 * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
1337
1338}
1339
1340//*********
1341
1342// Evaluate the Xi function of Gunion and Kunszt.
1343
1344double Sigma2ffbargmZWgmZW::xiGK( double tHnow, double uHnow) {
1345
1346 return - 4. * s3 * s4 + tHnow * (3. * tHnow + 4. * uHnow)
1347 + tHnow * tHnow * ( tHnow * uHnow / (s3 * s4)
1348 - 2. * (1. / s3 + 1./s4) * (tHnow + uHnow)
1349 + 2. * (s3 / s4 + s4 / s3) );
1350
1351}
1352
1353//*********
1354
1355// Evaluate the Xj function of Gunion and Kunszt.
1356
1357double Sigma2ffbargmZWgmZW::xjGK( double tHnow, double uHnow) {
1358
1359 return 8. * pow2(s3 + s4) - 8. * (s3 + s4) * (tHnow + uHnow)
1360 - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
1361 / (s3 * s4) - 2. * (1. / s3 + 1. / s4) * (tHnow + uHnow)
1362 + 2. * (s3 / s4 + s4 / s3) );
1363
1364}
1365
1366//**************************************************************************
1367
1368// Sigma2ffbar2gmZgmZ class.
1369// Cross section for f fbar -> gamma*/Z0 gamma*/Z0 (f is quark or lepton).
1370
1371//*********
1372
1373// Initialize process.
1374
1375void Sigma2ffbar2gmZgmZ::initProc() {
1376
1377 // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1378 gmZmode = Settings::mode("WeakZ0:gmZmode");
1379
1380 // Store Z0 mass and width for propagator.
1381 mRes = ParticleDataTable::m0(23);
1382 GammaRes = ParticleDataTable::mWidth(23);
1383 m2Res = mRes*mRes;
1384 GamMRat = GammaRes / mRes;
1385 thetaWRat = 1. / (16. * CoupEW::sin2thetaW() * CoupEW::cos2thetaW());
1386
1387 // Set pointer to particle properties and decay table.
1388 particlePtr = ParticleDataTable::particleDataPtr(23);
1389
1390}
1391
1392//*********
1393
1394// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1395
1396void Sigma2ffbar2gmZgmZ::sigmaKin() {
1397
1398 // Cross section part common for all incoming flavours.
1399 sigma0 = (M_PI / sH2) * pow2(alpEM) * 0.5
1400 * ( (tH2 + uH2 + 2. * (s3 + s4) * sH) / (tH * uH)
1401 - s3 * s4 * (1./tH2 + 1./uH2) );
1402
1403 // Common coupling factors at the resonance masses
1404 double alpEM3 = alphaEMPtr->alphaEM(s3);
1405 double alpS3 = alphaSPtr->alphaS(s3);
1406 double colQ3 = 3. * (1. + alpS3 / M_PI);
1407 double alpEM4 = alphaEMPtr->alphaEM(s4);
1408 double alpS4 = alphaSPtr->alphaS(s4);
1409 double colQ4 = 3. * (1. + alpS4 / M_PI);
1410
1411 // Reset quantities to sum. Declare variables in loop.
1412 gamSum3 = 0.;
1413 intSum3 = 0.;
1414 resSum3 = 0.;
1415 gamSum4 = 0.;
1416 intSum4 = 0.;
1417 resSum4 = 0.;
1418 int onMode;
1419 double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
1420
1421 // Loop over all Z0 decay channels.
1422 for (int i = 0; i < particlePtr->decay.size(); ++i) {
1423 int idAbs = abs( particlePtr->decay[i].product(0) );
1424
1425 // Only contributions from three fermion generations, except top.
1426 if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
1427 mf = ParticleDataTable::m0(idAbs);
1428 onMode = particlePtr->decay[i].onMode();
1429
1430 // First Z0: check that above threshold. Phase space.
1431 if (m3 > 2. * mf + MASSMARGIN) {
1432 mr = pow2(mf / m3);
1433 betaf = sqrtpos(1. - 4. * mr);
1434 psvec = betaf * (1. + 2. * mr);
1435 psaxi = pow3(betaf);
1436
1437 // First Z0: combine phase space with couplings.
1438 ef2 = CoupEW::ef2(idAbs) * psvec;
1439 efvf = CoupEW::efvf(idAbs) * psvec;
1440 vf2af2 = CoupEW::vf2(idAbs) * psvec + CoupEW::af2(idAbs) * psaxi;
1441 colf = (idAbs < 6) ? colQ3 : 1.;
1442
1443 // First Z0: store sum of combinations for open outstate channels.
1444 if (onMode == 1 || onMode == 2) {
1445 gamSum3 += colf * ef2;
1446 intSum3 += colf * efvf;
1447 resSum3 += colf * vf2af2;
1448 }
1449 }
1450
1451 // Second Z0: check that above threshold. Phase space.
1452 if (m4 > 2. * mf + MASSMARGIN) {
1453 mr = pow2(mf / m4);
1454 betaf = sqrtpos(1. - 4. * mr);
1455 psvec = betaf * (1. + 2. * mr);
1456 psaxi = pow3(betaf);
1457
1458 // Second Z0: combine phase space with couplings.
1459 ef2 = CoupEW::ef2(idAbs) * psvec;
1460 efvf = CoupEW::efvf(idAbs) * psvec;
1461 vf2af2 = CoupEW::vf2(idAbs) * psvec + CoupEW::af2(idAbs) * psaxi;
1462 colf = (idAbs < 6) ? colQ4 : 1.;
1463
1464 // Second Z0: store sum of combinations for open outstate channels.
1465 if (onMode == 1 || onMode == 2) {
1466 gamSum4 += colf * ef2;
1467 intSum4 += colf * efvf;
1468 resSum4 += colf * vf2af2;
1469 }
1470 }
1471
1472 // End loop over fermions.
1473 }
1474 }
1475
1476 // First Z0: calculate prefactors for gamma/interference/Z0 terms.
1477 gamProp3 = 4. * alpEM3 / (3. * M_PI * s3);
1478 intProp3 = gamProp3 * 2. * thetaWRat * s3 * (s3 - m2Res)
1479 / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1480 resProp3 = gamProp3 * pow2(thetaWRat * s3)
1481 / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1482
1483 // First Z0: optionally only keep gamma* or Z0 term.
1484 if (gmZmode == 1) {intProp3 = 0.; resProp3 = 0.;}
1485 if (gmZmode == 2) {gamProp3 = 0.; intProp3 = 0.;}
1486
1487 // Second Z0: calculate prefactors for gamma/interference/Z0 terms.
1488 gamProp4 = 4. * alpEM4 / (3. * M_PI * s4);
1489 intProp4 = gamProp4 * 2. * thetaWRat * s4 * (s4 - m2Res)
1490 / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1491 resProp4 = gamProp4 * pow2(thetaWRat * s4)
1492 / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1493
1494 // Second Z0: optionally only keep gamma* or Z0 term.
1495 if (gmZmode == 1) {intProp4 = 0.; resProp4 = 0.;}
1496 if (gmZmode == 2) {gamProp4 = 0.; intProp4 = 0.;}
1497
1498}
1499
1500//*********
1501
1502// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1503
1504double Sigma2ffbar2gmZgmZ::sigmaHat() {
1505
1506 // Charge/2, left- and righthanded couplings for in-fermion.
1507 int idAbs = abs(id1);
1508 double ei = 0.5 * CoupEW::ef(idAbs);
1509 double li = CoupEW::lf(idAbs);
1510 double ri = CoupEW::rf(idAbs);
1511
1512 // Combine left/right gamma, interference and Z0 parts for each Z0.
1513 double left3 = ei * ei * gamProp3 * gamSum3
1514 + ei * li * intProp3 * intSum3
1515 + li * li * resProp3 * resSum3;
1516 double right3 = ei * ei * gamProp3 * gamSum3
1517 + ei * ri * intProp3 * intSum3
1518 + ri * ri * resProp3 * resSum3;
1519 double left4 = ei * ei * gamProp4 * gamSum4
1520 + ei * li * intProp4 * intSum4
1521 + li * li * resProp4 * resSum4;
1522 double right4 = ei * ei * gamProp4 * gamSum4
1523 + ei * ri * intProp4 * intSum4
1524 + ri * ri * resProp4 * resSum4;
1525
1526 // Combine left- and right-handed couplings for the two Z0's.
1527 double sigma = sigma0 * (left3 * left4 + right3 * right4);
1528
1529 // Correct for the running-width Z0 propagators weight in PhaseSpace.
1530 sigma /= (runBW3 * runBW4);
1531
1532 // Initial-state colour factor. Answer.
1533 if (idAbs < 9) sigma /= 3.;
1534 return sigma;
1535
1536}
1537
1538//*********
1539
1540// Select identity, colour and anticolour.
1541
1542void Sigma2ffbar2gmZgmZ::setIdColAcol() {
1543
1544 // Flavours trivial.
1545 setId( id1, id2, 23, 23);
1546
1547 // Colour flow topologies. Swap when antiquarks.
1548 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1549 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1550 if (id1 < 0) swapColAcol();
1551
1552}
1553
1554//*********
1555
1556// Evaluate correlated decay flavours of the two gamma*/Z0.
1557// Unique complication, caused by gamma*/Z0 mix different left/right.
1558
1559double Sigma2ffbar2gmZgmZ::weightDecayFlav( Event& process) {
1560
1561 // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1562 i1 = (process[3].id() < 0) ? 3 : 4;
1563 i2 = 7 - i1;
1564 i3 = (process[7].id() > 0) ? 7 : 8;
1565 i4 = 15 - i3;
1566 i5 = (process[9].id() > 0) ? 9 : 10;
1567 i6 = 19 - i5;
1568
1569 // Charge/2, left- and righthanded couplings for in- and out-fermions.
1570 int idAbs = process[i1].idAbs();
1571 double ei = 0.5 * CoupEW::ef(idAbs);
1572 double li = CoupEW::lf(idAbs);
1573 double ri = CoupEW::rf(idAbs);
1574 idAbs = process[i3].idAbs();
1575 double e3 = 0.5 * CoupEW::ef(idAbs);
1576 double l3 = CoupEW::lf(idAbs);
1577 double r3 = CoupEW::rf(idAbs);
1578 idAbs = process[i5].idAbs();
1579 double e4 = 0.5 * CoupEW::ef(idAbs);
1580 double l4 = CoupEW::lf(idAbs);
1581 double r4 = CoupEW::rf(idAbs);
1582
1583 // Left- and righthanded couplings combined with propagators.
1584 c3LL = ei * ei * gamProp3 * e3 * e3
1585 + ei * li * intProp3 * e3 * l3
1586 + li * li * resProp3 * l3 * l3;
1587 c3LR = ei * ei * gamProp3 * e3 * e3
1588 + ei * li * intProp3 * e3 * r3
1589 + li * li * resProp3 * r3 * r3;
1590 c3RL = ei * ei * gamProp3 * e3 * e3
1591 + ei * ri * intProp3 * e3 * l3
1592 + ri * ri * resProp3 * l3 * l3;
1593 c3RR = ei * ei * gamProp3 * e3 * e3
1594 + ei * ri * intProp3 * e3 * r3
1595 + ri * ri * resProp3 * r3 * r3;
1596 c4LL = ei * ei * gamProp4 * e4 * e4
1597 + ei * li * intProp4 * e4 * l4
1598 + li * li * resProp4 * l4 * l4;
1599 c4LR = ei * ei * gamProp4 * e4 * e4
1600 + ei * li * intProp4 * e4 * r4
1601 + li * li * resProp4 * r4 * r4;
1602 c4RL = ei * ei * gamProp4 * e4 * e4
1603 + ei * ri * intProp4 * e4 * l4
1604 + ri * ri * resProp4 * l4 * l4;
1605 c4RR = ei * ei * gamProp4 * e4 * e4
1606 + ei * ri * intProp4 * e4 * r4
1607 + ri * ri * resProp4 * r4 * r4;
1608
1609 // Flavour weight and maximum.
1610 flavWt = (c3LL + c3LR) * (c4LL + c4LR) + (c3RL + c3RR) * (c4RL + c4RR);
1611 double flavWtMax = (c3LL + c3LR + c3RL + c3RR) * (c4LL + c4LR + c4RL + c4RR);
1612
1613 // Done.
1614 return flavWt / flavWtMax;
1615
1616}
1617
1618//*********
1619
1620// Evaluate weight for decay angles of the two gamma*/Z0.
1621
1622double Sigma2ffbar2gmZgmZ::weightDecay( Event& process, int iResBeg,
1623 int iResEnd) {
1624
1625 // Two resonance decays, but with common weight.
1626 if (iResBeg != 5 || iResEnd != 6) return 1.;
1627
1628 // Set up four-products and internal products.
1629 setupProd( process, i1, i2, i3, i4, i5, i6);
1630
1631 // Flip tHat and uHat if first incoming is fermion.
1632 double tHres = tH;
1633 double uHres = uH;
1634 if (process[3].id() > 0) swap( tHres, uHres);
1635
1636 // Kinematics factors (norm(x) = |x|^2).
1637 double fGK135 = norm( fGK( 1, 2, 3, 4, 5, 6) / tHres
1638 + fGK( 1, 2, 5, 6, 3, 4) / uHres );
1639 double fGK145 = norm( fGK( 1, 2, 4, 3, 5, 6) / tHres
1640 + fGK( 1, 2, 5, 6, 4, 3) / uHres );
1641 double fGK136 = norm( fGK( 1, 2, 3, 4, 6, 5) / tHres
1642 + fGK( 1, 2, 6, 5, 3, 4) / uHres );
1643 double fGK146 = norm( fGK( 1, 2, 4, 3, 6, 5) / tHres
1644 + fGK( 1, 2, 6, 5, 4, 3) / uHres );
1645 double fGK253 = norm( fGK( 2, 1, 5, 6, 3, 4) / tHres
1646 + fGK( 2, 1, 3, 4, 5, 6) / uHres );
1647 double fGK263 = norm( fGK( 2, 1, 6, 5, 3, 4) / tHres
1648 + fGK( 2, 1, 3, 4, 6, 5) / uHres );
1649 double fGK254 = norm( fGK( 2, 1, 5, 6, 4, 3) / tHres
1650 + fGK( 2, 1, 4, 3, 5, 6) / uHres );
1651 double fGK264 = norm( fGK( 2, 1, 6, 5, 4, 3) / tHres
1652 + fGK( 2, 1, 4, 3, 6, 5) / uHres );
1653
1654 // Weight and maximum.
1655 double wt = c3LL * c4LL * fGK135 + c3LR * c4LL * fGK145
1656 + c3LL * c4LR * fGK136 + c3LR * c4LR * fGK146
1657 + c3RL * c4RL * fGK253 + c3RR * c4RL * fGK263
1658 + c3RL * c4RR * fGK254 + c3RR * c4RR * fGK264;
1659 double wtMax = 16. * s3 * s4 * flavWt
1660 * ( (tHres*tHres + uHres*uHres + 2. * sH * (s3 + s4)) / (tHres * uHres)
1661 - s3 * s4 * (1. / (tHres*tHres) + 1. / (uHres*uHres)) );
1662
1663 // Done.
1664 return wt / wtMax;
1665
1666}
1667
1668//**************************************************************************
1669
1670// Sigma2ffbar2ZW class.
1671// Cross section for f fbar' -> W+ W- (f is quark or lepton).
1672
1673//*********
1674
1675// Initialize process.
1676
1677void Sigma2ffbar2ZW::initProc() {
1678
1679 // Store W+- mass and width for propagator.
1680 mW = ParticleDataTable::m0(24);
1681 widW = ParticleDataTable::mWidth(24);
1682 mWS = mW*mW;
1683 mwWS = pow2(mW * widW);
1684
1685 // Left-handed couplings for up/nu- and down/e-type quarks.
1686 lun = (hasLeptonBeams) ? CoupEW::lf(12) : CoupEW::lf(2);
1687 lde = (hasLeptonBeams) ? CoupEW::lf(11) : CoupEW::lf(1);
1688
1689 // Common weak coupling factor.
1690 sin2thetaW = CoupEW::sin2thetaW();
1691 cos2thetaW = CoupEW::cos2thetaW();
1692 thetaWRat = 1. / (4. * cos2thetaW);
1693 thetaWpt = (9. - 8. * sin2thetaW) / 4.;
1694 thetaWmm = (8. * sin2thetaW - 6.) / 4.;
1695
1696 // Secondary open width fractions.
1697 openFracPos = ParticleDataTable::resOpenFrac(23, 24);
1698 openFracNeg = ParticleDataTable::resOpenFrac(23, -24);
1699
1700}
1701
1702//*********
1703
1704// Evaluate sigmaHat(sHat), part independent of incoming flavour.
1705
1706void Sigma2ffbar2ZW::sigmaKin() {
1707
1708 // Evaluate cross section.
1709 double resBW = 1. / (pow2(sH - mWS) + mwWS);
1710 sigma0 = (M_PI / sH2) * 0.5 * pow2(alpEM / sin2thetaW);
1711 sigma0 *= sH * resBW * (thetaWpt * pT2 + thetaWmm * (s3 + s4))
1712 + (sH - mWS) * resBW * sH * (pT2 - s3 - s4) * (lun / tH - lde / uH)
1713 + thetaWRat * sH * pT2 * ( lun*lun / tH2 + lde*lde / uH2 )
1714 + 2. * thetaWRat * sH * (s3 + s4) * lun * lde / (tH * uH);
1715
1716 // Protect against slightly negative cross sections,
1717 // probably caused by addition of width to the W propagator.
1718 sigma0 = max(0., sigma0);
1719
1720}
1721
1722//*********
1723
1724// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1725
1726double Sigma2ffbar2ZW::sigmaHat() {
1727
1728 // CKM and colour factors.
1729 double sigma = sigma0;
1730 if (abs(id1) < 9) sigma *= VCKM::V2id(abs(id1), abs(id2)) / 3.;
1731
1732 // Corrections for secondary widths in Z0 and W+- decays.
1733 int idUp = (abs(id1)%2 == 0) ? id1 : id2;
1734 sigma *= (idUp > 0) ? openFracPos : openFracNeg;
1735
1736 // Answer.
1737 return sigma;
1738
1739}
1740
1741//*********
1742
1743// Select identity, colour and anticolour.
1744
1745void Sigma2ffbar2ZW::setIdColAcol() {
1746
1747 // Sign of outgoing W.
1748 int sign = 1 - 2 * (abs(id1)%2);
1749 if (id1 < 0) sign = -sign;
1750 setId( id1, id2, 23, 24 * sign);
1751
1752 // tHat is defined between (f, W-) or (fbar, W+),
1753 // so OK for u/ubar on side 1, but must swap tHat <-> uHat if d/dbar.
1754 if (abs(id1)%2 == 1) swapTU = true;
1755
1756 // Colour flow topologies. Swap when antiquarks.
1757 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1758 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1759 if (id1 < 0) swapColAcol();
1760
1761}
1762
1763//*********
1764
1765// Evaluate weight for Z0 and W+- decay angles.
1766
1767double Sigma2ffbar2ZW::weightDecay( Event& process, int iResBeg, int iResEnd) {
1768
1769 // Two resonance decays, but with common weight.
1770 if (iResBeg != 5 || iResEnd != 6) return 1.;
1771
1772 // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6)
1773 // with f' fbar' from W+- and f" fbar" from Z0 (note flip Z0 <-> W+-).
1774 int i1 = (process[3].id() < 0) ? 3 : 4;
1775 int i2 = 7 - i1;
1776 int i3 = (process[9].id() > 0) ? 9 : 10;
1777 int i4 = 19 - i3;
1778 int i5 = (process[7].id() > 0) ? 7 : 8;
1779 int i6 = 15 - i5;
1780
1781 // Set up four-products and internal products.
1782 setupProd( process, i1, i2, i3, i4, i5, i6);
1783
1784 // Swap tHat and uHat if incoming fermion is downtype.
1785 double tHres = tH;
1786 double uHres = uH;
1787 if (process[i2].id()%2 == 1) swap( tHres, uHres);
1788
1789 // Couplings of incoming (anti)fermions and outgoing from Z0.
1790 int idAbs = process[i1].idAbs();
1791 double ai = CoupEW::af(idAbs);
1792 double li1 = CoupEW::lf(idAbs);
1793 idAbs = process[i2].idAbs();
1794 double li2 = CoupEW::lf(idAbs);
1795 idAbs = process[i5].idAbs();
1796 double l4 = CoupEW::lf(idAbs);
1797 double r4 = CoupEW::rf(idAbs);
1798
1799 // W propagator/interference factor.
1800 double Wint = cos2thetaW * (sH - mWS) / (pow2(sH - mWS) + mwWS);
1801
1802 // Combinations of couplings and kinematics (norm(x) = |x|^2).
1803 double aWZ = li2 / tHres - 2. * Wint * ai;
1804 double bWZ = li1 / uHres + 2. * Wint * ai;
1805 double fGK135 = norm( aWZ * fGK( 1, 2, 3, 4, 5, 6)
1806 + bWZ * fGK( 1, 2, 5, 6, 3, 4) );
1807 double fGK136 = norm( aWZ * fGK( 1, 2, 3, 4, 6, 5)
1808 + bWZ * fGK( 1, 2, 6, 5, 3, 4) );
1809 double xiT = xiGK( tHres, uHres);
1810 double xiU = xiGK( uHres, tHres);
1811 double xjTU = xjGK( tHres, uHres);
1812
1813 // Weight and maximum weight.
1814 double wt = l4*l4 * fGK135 + r4*r4 * fGK136;
1815 double wtMax = 4. * s3 * s4 * (l4*l4 + r4*r4)
1816 * (aWZ * aWZ * xiT + bWZ * bWZ * xiU + aWZ * bWZ * xjTU);
1817
1818 // Done.
1819 return wt / wtMax;
1820
1821}
1822
1823//**************************************************************************
1824
1825// Sigma2ffbar2WW class.
1826// Cross section for f fbar -> W- W+ (f is quark or lepton).
1827
1828//*********
1829
1830// Initialize process.
1831
1832void Sigma2ffbar2WW::initProc() {
1833
1834 // Store Z0 mass and width for propagator. Common coupling factor.
1835 mZ = ParticleDataTable::m0(23);
1836 widZ = ParticleDataTable::mWidth(23);
1837 mZS = mZ*mZ;
1838 mwZS = pow2(mZ * widZ);
1839 thetaWRat = 1. / (4. * CoupEW::sin2thetaW());
1840
1841 // Secondary open width fraction.
1842 openFracPair = ParticleDataTable::resOpenFrac(24, -24);
1843
1844}
1845
1846//*********
1847
1848// Evaluate sigmaHat(sHat), part independent of incoming flavour.
1849
1850void Sigma2ffbar2WW::sigmaKin() {
1851
1852 // Cross section part common for all incoming flavours.
1853 sigma0 = (M_PI / sH2) * pow2(alpEM);
1854
1855 // Z0 propagator and gamma*/Z0 interference.
1856 double Zprop = sH2 / (pow2(sH - mZS) + mwZS);
1857 double Zint = Zprop * (1. - mZS / sH);
1858
1859 // Common coupling factors (g = gamma*, Z = Z0, f = t-channel fermion).
1860 cgg = 0.5;
1861 cgZ = thetaWRat * Zint;
1862 cZZ = 0.5 * pow2(thetaWRat) * Zprop;
1863 cfg = thetaWRat;
1864 cfZ = pow2(thetaWRat) * Zint;
1865 cff = pow2(thetaWRat);
1866
1867 // Kinematical functions.
1868 double rat34 = sH * (2. * (s3 + s4) + pT2) / (s3 * s4);
1869 double lambdaS = pow2(sH - s3 - s4) - 4. * s3 * s4;
1870 double intA = (sH - s3 - s4) * rat34 / sH;
1871 double intB = 4. * (s3 + s4 - pT2);
1872 gSS = (lambdaS * rat34 + 12. * sH * pT2) / sH2;
1873 gTT = rat34 + 4. * sH * pT2 / tH2;
1874 gST = intA + intB / tH;
1875 gUU = rat34 + 4. * sH * pT2 / uH2;
1876 gSU = intA + intB / uH;
1877
1878}
1879
1880//*********
1881
1882// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1883
1884double Sigma2ffbar2WW::sigmaHat() {
1885
1886 // Flavour-specific couplings.
1887 int idAbs = abs(id1);
1888 double ei = CoupEW::ef(idAbs);
1889 double vi = CoupEW::vf(idAbs);
1890 double ai = CoupEW::af(idAbs);
1891
1892 // Combine, with different cases for up- and down-type in-flavours.
1893 double sigma = sigma0;
1894 sigma *= (idAbs%2 == 1)
1895 ? (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
1896 + (cfg * ei + cfZ * (vi + ai)) * gST + cff * gTT
1897 : (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
1898 - (cfg * ei + cfZ * (vi + ai)) * gSU + cff * gUU;
1899
1900 // Initial-state colour factor. Correction for secondary widths. Answer.
1901 if (idAbs < 9) sigma /= 3.;
1902 sigma *= openFracPair;
1903 return sigma;
1904
1905}
1906
1907//*********
1908
1909// Select identity, colour and anticolour.
1910
1911void Sigma2ffbar2WW::setIdColAcol() {
1912
1913 // Always order W- W+, i.e. W- first.
1914 setId( id1, id2, -24, 24);
1915
1916 // tHat is defined between (f, W-) or (fbar, W+),
1917 if (id1 < 0) swapTU = true;
1918
1919 // Colour flow topologies. Swap when antiquarks.
1920 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1921 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1922 if (id1 < 0) swapColAcol();
1923
1924}
1925
1926//*********
1927
1928// Evaluate weight for W+ and W- decay angles.
1929
1930double Sigma2ffbar2WW::weightDecay( Event& process, int iResBeg, int iResEnd) {
1931
1932 // Two resonance decays, but with common weight.
1933 if (iResBeg != 5 || iResEnd != 6) return 1.;
1934
1935 // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1936 // with f' fbar' from W- and f" fbar" from W+.
1937 int i1 = (process[3].id() < 0) ? 3 : 4;
1938 int i2 = 7 - i1;
1939 int i3 = (process[7].id() > 0) ? 7 : 8;
1940 int i4 = 15 - i3;
1941 int i5 = (process[9].id() > 0) ? 9 : 10;
1942 int i6 = 19 - i5;
1943
1944 // Set up four-products and internal products.
1945 setupProd( process, i1, i2, i3, i4, i5, i6);
1946
1947 // tHat and uHat of fbar f -> W- W+ opposite to previous convention.
1948 double tHres = uH;
1949 double uHres = tH;
1950
1951 // Couplings of incoming (anti)fermion.
1952 int idAbs = process[i1].idAbs();
1953 double ai = CoupEW::af(idAbs);
1954 double li = CoupEW::lf(idAbs);
1955 double ri = CoupEW::rf(idAbs);
1956
1957 // gamma*/Z0 propagator/interference factor.
1958 double Zint = mZS * (sH - mZS) / (pow2(sH - mZS) + mwZS);
1959
1960 // Combinations of couplings and kinematics (norm(x) = |x|^2).
1961 double dWW = (li * Zint + ai) / sH;
1962 double aWW = dWW + 0.5 * (ai + 1.) / tHres;
1963 double bWW = dWW + 0.5 * (ai - 1.) / uHres;
1964 double cWW = ri * Zint / sH;
1965 double fGK135 = norm( aWW * fGK( 1, 2, 3, 4, 5, 6)
1966 - bWW * fGK( 1, 2, 5, 6, 3, 4) );
1967 double fGK253 = norm( cWW * ( fGK( 2, 1, 5, 6, 3, 4)
1968 - fGK( 2, 1, 3, 4, 5, 6) ) );
1969 double xiT = xiGK( tHres, uHres);
1970 double xiU = xiGK( uHres, tHres);
1971 double xjTU = xjGK( tHres, uHres);
1972
1973 // Weight and maximum weight.
1974 double wt = fGK135 + fGK253;
1975 double wtMax = 4. * s3 * s4
1976 * ( aWW * aWW * xiT + bWW * bWW * xiU - aWW * bWW * xjTU
1977 + cWW * cWW * (xiT + xiU - xjTU) );
1978
1979 // Done.
1980 return wt / wtMax;
1981}
1982
1983//**************************************************************************
1984
1985// Sigma2ffbargmZggm class.
1986// Collects common methods for f fbar -> gamma*/Z0 g/gamma and permutations.
1987
1988//*********
1989
1990// Initialize process.
1991
1992void Sigma2ffbargmZggm::initProc() {
1993
1994 // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1995 gmZmode = Settings::mode("WeakZ0:gmZmode");
1996
1997 // Store Z0 mass and width for propagator.
1998 mRes = ParticleDataTable::m0(23);
1999 GammaRes = ParticleDataTable::mWidth(23);
2000 m2Res = mRes*mRes;
2001 GamMRat = GammaRes / mRes;
2002 thetaWRat = 1. / (16. * CoupEW::sin2thetaW() * CoupEW::cos2thetaW());
2003
2004 // Set pointer to particle properties and decay table.
2005 particlePtr = ParticleDataTable::particleDataPtr(23);
2006
2007}
2008
2009//*********
2010
2011// Evaluate sum of flavour couplings times phase space.
2012
2013void Sigma2ffbargmZggm::flavSum() {
2014
2015 // Coupling factors for Z0 subsystem.
2016 double alpSZ = alphaSPtr->alphaS(s3);
2017 double colQZ = 3. * (1. + alpSZ / M_PI);
2018
2019 // Reset quantities to sum. Declare variables in loop.
2020 gamSum = 0.;
2021 intSum = 0.;
2022 resSum = 0.;
2023 int onMode;
2024 double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
2025
2026 // Loop over all Z0 decay channels.
2027 for (int i = 0; i < particlePtr->decay.size(); ++i) {
2028 int idAbs = abs( particlePtr->decay[i].product(0) );
2029
2030 // Only contributions from three fermion generations, except top.
2031 if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
2032 mf = ParticleDataTable::m0(idAbs);
2033
2034 // Check that above threshold. Phase space.
2035 if (m3 > 2. * mf + MASSMARGIN) {
2036 mr = pow2(mf / m3);
2037 betaf = sqrtpos(1. - 4. * mr);
2038 psvec = betaf * (1. + 2. * mr);
2039 psaxi = pow3(betaf);
2040
2041 // Combine phase space with couplings.
2042 ef2 = CoupEW::ef2(idAbs) * psvec;
2043 efvf = CoupEW::efvf(idAbs) * psvec;
2044 vf2af2 = CoupEW::vf2(idAbs) * psvec + CoupEW::af2(idAbs) * psaxi;
2045 colf = (idAbs < 6) ? colQZ : 1.;
2046
2047 // Store sum of combinations. For outstate only open channels.
2048 onMode = particlePtr->decay[i].onMode();
2049 if (onMode == 1 || onMode == 2) {
2050 gamSum += colf * ef2;
2051 intSum += colf * efvf;
2052 resSum += colf * vf2af2;
2053 }
2054
2055 // End loop over fermions.
2056 }
2057 }
2058 }
2059
2060 // Done. Return values in gamSum, intSum and resSum.
2061
2062}
2063
2064//*********
2065
2066// Calculate common parts of gamma/interference/Z0 propagator terms.
2067
2068void Sigma2ffbargmZggm::propTerm() {
2069
2070 // Calculate prefactors for gamma/interference/Z0 cross section terms.
2071 gamProp = 4. * alpEM / (3. * M_PI * s3);
2072 intProp = gamProp * 2. * thetaWRat * s3 * (s3 - m2Res)
2073 / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2074 resProp = gamProp * pow2(thetaWRat * s3)
2075 / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2076
2077 // Optionally only keep gamma* or Z0 term.
2078 if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
2079 if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
2080
2081}
2082
2083//*********
2084
2085// Evaluate weight for gamma*/Z0 decay angle.
2086
2087double Sigma2ffbargmZggm::weightDecay( Event& process, int iResBeg, int iResEnd) {
2088
2089 // Z should sit in entry 5 and one more parton in entry 6.
2090 if (iResBeg != 5 || iResEnd != 6) return 1.;
2091
2092 // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2093 // where f' fbar' come from gamma*/Z0 decay.
2094 int i1, i2;
2095 int i3 = (process[7].id() > 0) ? 7 : 8;
2096 int i4 = 15 - i3;
2097
2098 // Order so that fbar(1) f(2) -> gamma*/Z0 g/gamma.
2099 if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2100 i1 = (process[3].id() < 0) ? 3 : 4;
2101 i2 = 7 - i1;
2102
2103 // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) gamma*/Z0.
2104 } else if (process[3].idAbs() < 20) {
2105 i1 = (process[3].id() < 0) ? 3 : 6;
2106 i2 = 9 - i1;
2107 } else {
2108 i1 = (process[4].id() < 0) ? 4 : 6;
2109 i2 = 10 - i1;
2110 }
2111
2112 // Charge/2, left- and righthanded couplings for in- and out-fermion.
2113 int id1Abs = process[i1].idAbs();
2114 double ei = 0.5 * CoupEW::ef(id1Abs);
2115 double li = CoupEW::lf(id1Abs);
2116 double ri = CoupEW::rf(id1Abs);
2117 int id3Abs = process[i3].idAbs();
2118 double ef = 0.5 * CoupEW::ef(id3Abs);
2119 double lf = CoupEW::lf(id3Abs);
2120 double rf = CoupEW::rf(id3Abs);
2121
2122 // Combinations of left/right for in/out, gamma*/interference/Z0.
2123 double clilf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*lf
2124 + li*li * resProp * lf*lf;
2125 double clirf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*rf
2126 + li*li * resProp * rf*rf;
2127 double crilf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*lf
2128 + ri*ri * resProp * lf*lf;
2129 double crirf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*rf
2130 + ri*ri * resProp * rf*rf;
2131
2132 // Evaluate four-vector products.
2133 double p13 = process[i1].p() * process[i3].p();
2134 double p14 = process[i1].p() * process[i4].p();
2135 double p23 = process[i2].p() * process[i3].p();
2136 double p24 = process[i2].p() * process[i4].p();
2137
2138 // Calculate weight and its maximum.
2139 double wt = (clilf + crirf) * (p13*p13 + p24*p24)
2140 + (clirf + crilf) * (p14*p14 + p23*p23) ;
2141 double wtMax = (clilf + clirf + crilf + crirf)
2142 * (pow2(p13 + p14) + pow2(p23 + p24));
2143
2144 // Done.
2145 return (wt / wtMax);
2146
2147}
2148
2149//**************************************************************************
2150
2151// Sigma2qqbar2gmZg class.
2152// Cross section for q qbar -> gamma*/Z0 g.
2153
2154//*********
2155
2156// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2157
2158void Sigma2qqbar2gmZg::sigmaKin() {
2159
2160 // Cross section part common for all incoming flavours.
2161 sigma0 = (M_PI / sH2) * (alpEM * alpS)
2162 * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2163
2164 // Calculate flavour sums for final state.
2165 flavSum();
2166
2167 // Calculate prefactors for gamma/interference/Z0 cross section terms.
2168 propTerm();
2169
2170}
2171
2172//*********
2173
2174// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2175
2176double Sigma2qqbar2gmZg::sigmaHat() {
2177
2178 // Combine gamma, interference and Z0 parts.
2179 int idAbs = abs(id1);
2180 double sigma = sigma0
2181 * ( CoupEW::ef2(idAbs) * gamProp * gamSum
2182 + CoupEW::efvf(idAbs) * intProp * intSum
2183 + CoupEW::vf2af2(idAbs) * resProp * resSum);
2184
2185 // Correct for the running-width Z0 propagater weight in PhaseSpace.
2186 sigma /= runBW3;
2187
2188 // Answer.
2189 return sigma;
2190
2191}
2192
2193//*********
2194
2195// Select identity, colour and anticolour.
2196
2197void Sigma2qqbar2gmZg::setIdColAcol() {
2198
2199 // Flavours trivial.
2200 setId( id1, id2, 23, 21);
2201
2202 // Colour flow topologies. Swap when antiquarks.
2203 setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2204 if (id1 < 0) swapColAcol();
2205
2206}
2207
2208//**************************************************************************
2209
2210// Sigma2qg2gmZq class.
2211// Cross section for q g -> gamma*/Z0 q.
2212
2213//*********
2214
2215// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2216
2217void Sigma2qg2gmZq::sigmaKin() {
2218
2219 // Cross section part common for all incoming flavours.
2220 sigma0 = (M_PI / sH2) * (alpEM * alpS)
2221 * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2222
2223 // Calculate flavour sums for final state.
2224 flavSum();
2225
2226 // Calculate prefactors for gamma/interference/Z0 cross section terms.
2227 propTerm();
2228
2229}
2230
2231//*********
2232
2233// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2234
2235double Sigma2qg2gmZq::sigmaHat() {
2236
2237 // Combine gamma, interference and Z0 parts.
2238 int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2239 double sigma = sigma0
2240 * ( CoupEW::ef2(idAbs) * gamProp * gamSum
2241 + CoupEW::efvf(idAbs) * intProp * intSum
2242 + CoupEW::vf2af2(idAbs) * resProp * resSum);
2243
2244 // Correct for the running-width Z0 propagater weight in PhaseSpace.
2245 sigma /= runBW3;
2246
2247 // Answer.
2248 return sigma;
2249
2250}
2251
2252//*********
2253
2254// Select identity, colour and anticolour.
2255
2256void Sigma2qg2gmZq::setIdColAcol() {
2257
2258 // Flavour set up for q g -> gamma*/Z0 q.
2259 int idq = (id2 == 21) ? id1 : id2;
2260 setId( id1, id2, 23, idq);
2261
2262 // tH defined between f and f': must swap tHat <-> uHat if q g in.
2263 swapTU = (id2 == 21);
2264
2265 // Colour flow topologies. Swap when antiquarks.
2266 if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2267 else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2268 if (idq < 0) swapColAcol();
2269
2270}
2271
2272//**************************************************************************
2273
2274// Sigma2ffbar2gmZgm class.
2275// Cross section for f fbar -> gamma*/Z0 gamma.
2276
2277//*********
2278
2279// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2280
2281void Sigma2ffbar2gmZgm::sigmaKin() {
2282
2283 // Cross section part common for all incoming flavours.
2284 sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2285 * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2286
2287 // Calculate flavour sums for final state.
2288 flavSum();
2289
2290 // Calculate prefactors for gamma/interference/Z0 cross section terms.
2291 propTerm();
2292
2293
2294}
2295
2296//*********
2297
2298// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2299
2300double Sigma2ffbar2gmZgm::sigmaHat() {
2301
2302 // Combine gamma, interference and Z0 parts.
2303 int idAbs = abs(id1);
2304 double sigma = sigma0 * CoupEW::ef2(idAbs)
2305 * ( CoupEW::ef2(idAbs) * gamProp * gamSum
2306 + CoupEW::efvf(idAbs) * intProp * intSum
2307 + CoupEW::vf2af2(idAbs) * resProp * resSum);
2308
2309 // Correct for the running-width Z0 propagater weight in PhaseSpace.
2310 sigma /= runBW3;
2311
2312 // Colour factor. Answer.
2313 if (idAbs < 9) sigma /= 3.;
2314 return sigma;
2315
2316}
2317
2318//*********
2319
2320// Select identity, colour and anticolour.
2321
2322void Sigma2ffbar2gmZgm::setIdColAcol() {
2323
2324 // Flavours trivial.
2325 setId( id1, id2, 23, 22);
2326
2327 // Colour flow topologies. Swap when antiquarks.
2328 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2329 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2330 if (id1 < 0) swapColAcol();
2331
2332}
2333
2334//**************************************************************************
2335
2336// Sigma2fgm2gmZf class.
2337// Cross section for f gamma -> gamma*/Z0 f'.
2338
2339//*********
2340
2341// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2342
2343void Sigma2fgm2gmZf::sigmaKin() {
2344
2345 // Cross section part common for all incoming flavours.
2346 sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2347 * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (- sH * uH);
2348
2349 // Calculate flavour sums for final state.
2350 flavSum();
2351
2352 // Calculate prefactors for gamma/interference/Z0 cross section terms.
2353 propTerm();
2354
2355}
2356
2357//*********
2358
2359// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2360
2361double Sigma2fgm2gmZf::sigmaHat() {
2362
2363 // Combine gamma, interference and Z0 parts.
2364 int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
2365 double sigma = sigma0 * CoupEW::ef2(idAbs)
2366 * ( CoupEW::ef2(idAbs) * gamProp * gamSum
2367 + CoupEW::efvf(idAbs) * intProp * intSum
2368 + CoupEW::vf2af2(idAbs) * resProp * resSum);
2369
2370 // Correct for the running-width Z0 propagater weight in PhaseSpace.
2371 sigma /= runBW3;
2372
2373 // Answer.
2374 return sigma;
2375
2376}
2377
2378//*********
2379
2380// Select identity, colour and anticolour.
2381
2382void Sigma2fgm2gmZf::setIdColAcol() {
2383
2384 // Flavour set up for q gamma -> gamma*/Z0 q.
2385 int idq = (id2 == 22) ? id1 : id2;
2386 setId( id1, id2, 23, idq);
2387
2388 // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2389 swapTU = (id2 == 22);
2390
2391 // Colour flow topologies. Swap when antiquarks.
2392 if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2393 else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2394 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2395 if (idq < 0) swapColAcol();
2396
2397}
2398
2399//**************************************************************************
2400
2401// Sigma2ffbarWggm class.
2402// Collects common methods for f fbar -> W+- g/gamma and permutations.
2403
2404//*********
2405
2406// Evaluate weight for W+- decay angle.
2407
2408double Sigma2ffbarWggm::weightDecay( Event& process, int iResBeg, int iResEnd) {
2409
2410 // W should sit in entry 5 and one more parton in entry 6.
2411 if (iResBeg != 5 || iResEnd != 6) return 1.;
2412
2413 // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2414 // where f' fbar' come from W+- decay.
2415 int i1, i2;
2416 int i3 = (process[7].id() > 0) ? 7 : 8;
2417 int i4 = 15 - i3;
2418
2419 // Order so that fbar(1) f(2) -> W+- g/gamma.
2420 if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2421 i1 = (process[3].id() < 0) ? 3 : 4;
2422 i2 = 7 - i1;
2423
2424 // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) W+-.
2425 } else if (process[3].idAbs() < 20) {
2426 i1 = (process[3].id() < 0) ? 3 : 6;
2427 i2 = 9 - i1;
2428 } else {
2429 i1 = (process[4].id() < 0) ? 4 : 6;
2430 i2 = 10 - i1;
2431 }
2432
2433 // Evaluate four-vector products.
2434 double p13 = process[i1].p() * process[i3].p();
2435 double p14 = process[i1].p() * process[i4].p();
2436 double p23 = process[i2].p() * process[i3].p();
2437 double p24 = process[i2].p() * process[i4].p();
2438
2439 // Calculate weight and its maximum.
2440 double wt = pow2(p13) + pow2(p24);
2441 double wtMax = pow2(p13 + p14) + pow2(p23 + p24);
2442
2443 // Done.
2444 return (wt / wtMax);
2445
2446}
2447
2448//**************************************************************************
2449
2450// Sigma2qqbar2Wg class.
2451// Cross section for q qbar' -> W+- g.
2452
2453//*********
2454
2455// Initialize process.
2456
2457void Sigma2qqbar2Wg::initProc() {
2458
2459 // Secondary open width fractions, relevant for top (or heavier).
2460 openFracPos = ParticleDataTable::resOpenFrac(24);
2461 openFracNeg = ParticleDataTable::resOpenFrac(-24);
2462
2463}
2464
2465//*********
2466
2467// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2468
2469void Sigma2qqbar2Wg::sigmaKin() {
2470
2471 // Cross section part common for all incoming flavours.
2472 sigma0 = (M_PI / sH2) * (alpEM * alpS / CoupEW::sin2thetaW())
2473 * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2474
2475}
2476
2477//*********
2478
2479// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2480
2481double Sigma2qqbar2Wg::sigmaHat() {
2482
2483 // CKM factor. Secondary width for W+ or W-.
2484 double sigma = sigma0 * VCKM::V2id(abs(id1), abs(id2));
2485 int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2486 sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2487
2488 // Answer.
2489 return sigma;
2490
2491}
2492
2493//*********
2494
2495// Select identity, colour and anticolour.
2496
2497void Sigma2qqbar2Wg::setIdColAcol() {
2498
2499 // Sign of outgoing W.
2500 int sign = 1 - 2 * (abs(id1)%2);
2501 if (id1 < 0) sign = -sign;
2502 setId( id1, id2, 24 * sign, 21);
2503
2504 // Colour flow topologies. Swap when antiquarks.
2505 setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2506 if (id1 < 0) swapColAcol();
2507
2508}
2509
2510//**************************************************************************
2511
2512// Sigma2qg2Wq class.
2513// Cross section for q g -> W+- q'.
2514
2515//*********
2516
2517// Initialize process.
2518
2519void Sigma2qg2Wq::initProc() {
2520
2521 // Secondary open width fractions, relevant for top (or heavier).
2522 openFracPos = ParticleDataTable::resOpenFrac(24);
2523 openFracNeg = ParticleDataTable::resOpenFrac(-24);
2524
2525}
2526
2527//*********
2528
2529// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2530
2531void Sigma2qg2Wq::sigmaKin() {
2532
2533 // Cross section part common for all incoming flavours.
2534 sigma0 = (M_PI / sH2) * (alpEM * alpS / CoupEW::sin2thetaW())
2535 * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2536
2537}
2538
2539//*********
2540
2541// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2542
2543double Sigma2qg2Wq::sigmaHat() {
2544
2545 // CKM factor. Secondary width for W+ or W-.
2546 int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2547 double sigma = sigma0 * VCKM::V2sum(idAbs);
2548 int idUp = (id2 == 21) ? id1 : id2;
2549 if (idAbs%2 == 1) idUp = -idUp;
2550 sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2551
2552 // Answer.
2553 return sigma;
2554
2555}
2556
2557//*********
2558
2559// Select identity, colour and anticolour.
2560
2561void Sigma2qg2Wq::setIdColAcol() {
2562
2563 // Sign of outgoing W. Flavour of outgoing quark.
2564 int idq = (id2 == 21) ? id1 : id2;
2565 int sign = 1 - 2 * (abs(idq)%2);
2566 if (idq < 0) sign = -sign;
2567 id4 = VCKM::V2pick(idq);
2568
2569 // Flavour set up for q g -> W q.
2570 setId( id1, id2, 24 * sign, id4);
2571
2572 // tH defined between f and f': must swap tHat <-> uHat if q g in.
2573 swapTU = (id2 == 21);
2574
2575 // Colour flow topologies. Swap when antiquarks.
2576 if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2577 else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2578 if (idq < 0) swapColAcol();
2579
2580}
2581
2582//**************************************************************************
2583
2584// Sigma2ffbar2Wgm class.
2585// Cross section for f fbar' -> W+- gamma.
2586
2587//*********
2588
2589// Initialize process.
2590
2591void Sigma2ffbar2Wgm::initProc() {
2592
2593 // Secondary open width fractions, relevant for top (or heavier).
2594 openFracPos = ParticleDataTable::resOpenFrac(24);
2595 openFracNeg = ParticleDataTable::resOpenFrac(-24);
2596
2597}
2598
2599//*********
2600
2601// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2602
2603void Sigma2ffbar2Wgm::sigmaKin() {
2604
2605 // Cross section part common for all incoming flavours.
2606 sigma0 = (M_PI / sH2) * (alpEM*alpEM / CoupEW::sin2thetaW())
2607 * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2608}
2609
2610//*********
2611
2612// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2613
2614double Sigma2ffbar2Wgm::sigmaHat() {
2615
2616 // Extrafactor different for e nu and q qbar' instate.
2617 int id1Abs = abs(id1);
2618 int id2Abs = abs(id2);
2619 double chgUp = (id1Abs > 10) ? 0. : 2./3.;
2620 double sigma = sigma0 * pow2( chgUp - tH / (tH + uH) );
2621
2622 // CKM and colour factors. Secondary width for W+ or W-.
2623 if (id1Abs < 9) sigma *= VCKM::V2id(id1Abs, id2Abs) / 3.;
2624 int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2625 sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2626
2627 // Answer.
2628 return sigma;
2629
2630}
2631
2632//*********
2633
2634// Select identity, colour and anticolour.
2635
2636void Sigma2ffbar2Wgm::setIdColAcol() {
2637
2638 // Sign of outgoing W.
2639 int sign = 1 - 2 * (abs(id1)%2);
2640 if (id1 < 0) sign = -sign;
2641 setId( id1, id2, 24 * sign, 22);
2642
2643 // tH defined between (f,W-) or (fbar',W+).
2644 swapTU = (sign * id1 > 0);
2645
2646 // Colour flow topologies. Swap when antiquarks.
2647 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2648 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2649 if (id1 < 0) swapColAcol();
2650
2651}
2652
2653//**************************************************************************
2654
2655// Sigma2fgm2Wf class.
2656// Cross section for f gamma -> W+- f'.
2657
2658//*********
2659
2660// Initialize process.
2661
2662void Sigma2fgm2Wf::initProc() {
2663
2664 // Secondary open width fractions, relevant for top (or heavier).
2665 openFracPos = ParticleDataTable::resOpenFrac(24);
2666 openFracNeg = ParticleDataTable::resOpenFrac(-24);
2667
2668}
2669
2670//*********
2671
2672// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2673
2674void Sigma2fgm2Wf::sigmaKin() {
2675
2676 // Cross section part common for all incoming flavours.
2677 sigma0 = (M_PI / sH2) * (alpEM*alpEM / CoupEW::sin2thetaW())
2678 * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (pT2 * s3 - sH * uH);
2679
2680}
2681
2682//*********
2683
2684// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2685
2686double Sigma2fgm2Wf::sigmaHat() {
2687
2688 // Extrafactor dependent on charge of incoming fermion.
2689 int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
2690 double charge = (idAbs > 10) ? 1. : ( (idAbs%2 == 1) ? 1./3. : 2./3. );
2691 double sigma = sigma0 * pow2( charge - sH / (sH + uH) );
2692
2693 // CKM factor. Secondary width for W+ or W-.
2694 sigma *= VCKM::V2sum(idAbs);
2695 int idUp = (id2 == 22) ? id1 : id2;
2696 if (idAbs%2 == 1) idUp = -idUp;
2697 sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2698
2699 // Answer.
2700 return sigma;
2701
2702}
2703
2704//*********
2705
2706// Select identity, colour and anticolour.
2707
2708void Sigma2fgm2Wf::setIdColAcol() {
2709
2710 // Sign of outgoing W. Flavour of outgoing fermion.
2711 int idq = (id2 == 22) ? id1 : id2;
2712 int sign = 1 - 2 * (abs(idq)%2);
2713 if (idq < 0) sign = -sign;
2714 id4 = VCKM::V2pick(idq);
2715
2716 // Flavour set up for q gamma -> W q.
2717 setId( id1, id2, 24 * sign, id4);
2718
2719 // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2720 swapTU = (id2 == 22);
2721
2722 // Colour flow topologies. Swap when antiquarks.
2723 if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2724 else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2725 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2726 if (idq < 0) swapColAcol();
2727
2728}
2729
2730//**************************************************************************
2731
2732} // end namespace Pythia8