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