« back to all changes in this revision
Viewing changes to src/SigmaEW.cc
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Package Import Robot
- Author(s): Lifeng Sun
- Date: 2012-05-22 11:43:00 UTC
- Revision ID: package-import@ubuntu.com-20120522114300-0jvsv2vl4o2bo435
Tags: upstream-8.1.65
Import upstream version 8.1.65
- files added:
- examples
- examples/outref
- hepmcinterface
- htmldoc
- include
- lhapdfdummy
- phpdoc
- phpdoc/files
- rootexamples
- src
- xmldoc
added
removed
src/SigmaEW.cc
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
11
namespace 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
22
void 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
36
double 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
49
void 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
72
void 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
86
double 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
98
void 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
119
void 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
136
void 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
173
void 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
191
void 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
205
double 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
218
void 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
240
void 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
257
void 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
294
void 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
312
void 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
327
void 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
345
double 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
380
void 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
408
void 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
421
void 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
433
double 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
461
void 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
492
void 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
517
void 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
528
double 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
573
void 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
617
double 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
636
void 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
658
void 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
722
double 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
740
void 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
756
double 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
809
void 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
827
void 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
841
double 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
857
void 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
875
double 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
910
void 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
951
double 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
967
void 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
991
void 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
1030
void 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
1067
double 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
1103
void 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
1122
double 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
1141
void 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
1184
void 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
1220
double 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
1242
void 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
1275
double 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
1294
void 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
1341
complex 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
1353
double 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
1366
double 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
1384
void 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
1406
void 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
1516
double 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
1554
void 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
1571
double 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
1634
double 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
1689
void 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
1719
void 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
1761
double 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
1780
void 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
1802
double 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
1867
void 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
1885
void 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
1919
double 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
1946
void 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
1965
double 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
2027
void 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
2049
void 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
2105
void 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
2124
double 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
2196
void 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
2214
double 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
2235
void 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
2255
void 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
2273
double 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
2294
void 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
2319
void 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
2338
double 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
2360
void 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
2381
void 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
2399
double 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
2420
void 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
2446
double 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
2496
void 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
2508
void 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
2520
double 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
2536
void 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
2558
void 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
2570
void 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
2582
double 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
2600
void 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
2630
void 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
2642
void 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
2653
double 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
2675
void 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
2701
void 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
2713
void 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
2725
double 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
2747
void 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
2778
void 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
2809
void 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
2843
void 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
Loggerhead is a web-based interface for Breezy
Version: 2.0.1