1
// SigmaNewGaugeBosons.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.
6
// Function definitions (not found in the header) for the
7
// leptoquark simulation classes.
9
#include "SigmaNewGaugeBosons.h"
14
//==========================================================================
16
// Sigma1ffbarZprimeWprime class.
17
// Collects common methods for f fbar -> Z'/W' -> WW/WZ -> 4 fermions.
18
// Copied from SigmaEW for gauge-boson-pair production.
20
//--------------------------------------------------------------------------
22
// Calculate and store internal products.
24
void Sigma1ffbarZprimeWprime::setupProd( Event& process, int i1, int i2,
25
int i3, int i4, int i5, int i6) {
27
// Store incoming and outgoing momenta,
28
pRot[1] = process[i1].p();
29
pRot[2] = process[i2].p();
30
pRot[3] = process[i3].p();
31
pRot[4] = process[i4].p();
32
pRot[5] = process[i5].p();
33
pRot[6] = process[i6].p();
35
// Do random rotation to avoid accidental zeroes in HA expressions.
39
double thetaNow = acos(2. * rndmPtr->flat() - 1.);
40
double phiNow = 2. * M_PI * rndmPtr->flat();
41
for (int i = 1; i <= 6; ++i) {
42
pRot[i].rot( thetaNow, phiNow);
43
if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
47
// Calculate internal products.
48
for (int i = 1; i < 6; ++i) {
49
for (int j = i + 1; j <= 6; ++j) {
51
sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
52
/ pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
53
- sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
54
/ pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
55
hC[i][j] = conj( hA[i][j] );
57
hA[i][j] *= complex( 0., 1.);
58
hC[i][j] *= complex( 0., 1.);
60
hA[j][i] = - hA[i][j];
61
hC[j][i] = - hC[i][j];
67
//--------------------------------------------------------------------------
69
// Evaluate the F function of Gunion and Kunszt.
71
complex Sigma1ffbarZprimeWprime::fGK(int j1, int j2, int j3, int j4,
74
return 4. * hA[j1][j3] * hC[j2][j6]
75
* ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
79
//--------------------------------------------------------------------------
81
// Evaluate the Xi function of Gunion and Kunszt.
83
double Sigma1ffbarZprimeWprime::xiGK( double tHnow, double uHnow,
84
double s3now, double s4now) {
86
return - 4. * s3now * s4now + tHnow * (3. * tHnow + 4. * uHnow)
87
+ tHnow * tHnow * ( tHnow * uHnow / (s3now * s4now)
88
- 2. * (1. / s3now + 1./s4now) * (tHnow + uHnow)
89
+ 2. * (s3now / s4now + s4now / s3now) );
93
//--------------------------------------------------------------------------
95
// Evaluate the Xj function of Gunion and Kunszt.
97
double Sigma1ffbarZprimeWprime::xjGK( double tHnow, double uHnow,
98
double s3now, double s4now) {
100
return 8. * pow2(s3now + s4now) - 8. * (s3now + s4now) * (tHnow + uHnow)
101
- 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
102
/ (s3now * s4now) - 2. * (1. / s3now + 1. / s4now) * (tHnow + uHnow)
103
+ 2. * (s3now / s4now + s4now / s3now) );
107
//==========================================================================
109
// Sigma1ffbar2gmZZprime class.
110
// Cross section for f fbar -> gamma*/Z0/Z'0 (f is quark or lepton).
112
//--------------------------------------------------------------------------
114
// Initialize process.
116
void Sigma1ffbar2gmZZprime::initProc() {
118
// Allow to pick only parts of full gamma*/Z0/Z'0 expression.
119
gmZmode = settingsPtr->mode("Zprime:gmZmode");
121
// Store Z'0 mass and width for propagator.
122
mRes = particleDataPtr->m0(32);
123
GammaRes = particleDataPtr->mWidth(32);
125
GamMRat = GammaRes / mRes;
126
sin2tW = couplingsPtr->sin2thetaW();
127
cos2tW = 1. - sin2tW;
128
thetaWRat = 1. / (16. * sin2tW * cos2tW);
130
// Properties of Z0 resonance also needed.
131
mZ = particleDataPtr->m0(23);
132
GammaZ = particleDataPtr->mWidth(23);
134
GamMRatZ = GammaZ / mZ;
136
// Ensure that arrays initially are empty.
137
for (int i = 0; i < 20; ++i) afZp[i] = 0.;
138
for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
140
// Store first-generation axial and vector couplings.
141
afZp[1] = settingsPtr->parm("Zprime:ad");
142
afZp[2] = settingsPtr->parm("Zprime:au");
143
afZp[11] = settingsPtr->parm("Zprime:ae");
144
afZp[12] = settingsPtr->parm("Zprime:anue");
145
vfZp[1] = settingsPtr->parm("Zprime:vd");
146
vfZp[2] = settingsPtr->parm("Zprime:vu");
147
vfZp[11] = settingsPtr->parm("Zprime:ve");
148
vfZp[12] = settingsPtr->parm("Zprime:vnue");
150
// Second and third generation could be carbon copy of this...
151
if (settingsPtr->flag("Zprime:universality")) {
152
for (int i = 3; i <= 6; ++i) {
155
afZp[i+10] = afZp[i+8];
156
vfZp[i+10] = vfZp[i+8];
159
// ... or could have different couplings.
161
afZp[3] = settingsPtr->parm("Zprime:as");
162
afZp[4] = settingsPtr->parm("Zprime:ac");
163
afZp[5] = settingsPtr->parm("Zprime:ab");
164
afZp[6] = settingsPtr->parm("Zprime:at");
165
afZp[13] = settingsPtr->parm("Zprime:amu");
166
afZp[14] = settingsPtr->parm("Zprime:anumu");
167
afZp[15] = settingsPtr->parm("Zprime:atau");
168
afZp[16] = settingsPtr->parm("Zprime:anutau");
169
vfZp[3] = settingsPtr->parm("Zprime:vs");
170
vfZp[4] = settingsPtr->parm("Zprime:vc");
171
vfZp[5] = settingsPtr->parm("Zprime:vb");
172
vfZp[6] = settingsPtr->parm("Zprime:vt");
173
vfZp[13] = settingsPtr->parm("Zprime:vmu");
174
vfZp[14] = settingsPtr->parm("Zprime:vnumu");
175
vfZp[15] = settingsPtr->parm("Zprime:vtau");
176
vfZp[16] = settingsPtr->parm("Zprime:vnutau");
179
// Coupling for Z' -> W+ W- and decay angular admixture.
180
coupZpWW = settingsPtr->parm("Zprime:coup2WW");
181
anglesZpWW = settingsPtr->parm("Zprime:anglesWW");
183
// Set pointer to particle properties and decay table.
184
particlePtr = particleDataPtr->particleDataEntryPtr(32);
188
//--------------------------------------------------------------------------
190
// Evaluate sigmaHat(sHat), part independent of incoming flavour.
192
void Sigma1ffbar2gmZZprime::sigmaKin() {
194
// Common coupling factors.
195
double colQ = 3. * (1. + alpS / M_PI);
197
// Reset quantities to sum. Declare variables in loop.
205
double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf,
206
ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf;
208
// Loop over all open Z'0 decay channels.
209
for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
210
onMode = particlePtr->channel(i).onMode();
211
if (onMode != 1 && onMode != 2) continue;
212
idAbs = abs( particlePtr->channel(i).product(0) );
214
// Contributions from three fermion generations.
215
if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) {
216
mf = particleDataPtr->m0(idAbs);
218
// Check that above threshold.
219
if (mH > 2. * mf + MASSMARGIN) {
221
ps = sqrtpos(1. - 4. * mr);
223
// Couplings of gamma^*/Z^0/Z'^0 to final flavour
224
ef = couplingsPtr->ef(idAbs);
225
af = couplingsPtr->af(idAbs);
226
vf = couplingsPtr->vf(idAbs);
230
// Combine couplings with kinematical factors.
232
kinFacV = ps * (1. + 2. * mr);
233
ef2 = ef * ef * kinFacV;
234
efvf = ef * vf * kinFacV;
235
vaf2 = vf * vf * kinFacV + af * af * kinFacA;
236
efvpf = ef * vpf * kinFacV;
237
vafvapf = vf * vpf * kinFacV + af * apf * kinFacA;
238
vapf2 = vpf * vpf * kinFacV + apf * apf * kinFacA;
240
// Colour factor. Additionally secondary width for top.
241
colf = (idAbs < 9) ? colQ : 1.;
242
if (idAbs == 6) colf *= particleDataPtr->resOpenFrac(6, -6);
244
// Store sum of combinations.
245
gamSum += colf * ef2;
246
gamZSum += colf * efvf;
248
gamZpSum += colf * efvpf;
249
ZZpSum += colf * vafvapf;
250
ZpSum += colf * vapf2;
253
// Optional contribution from W+ W-.
254
} else if (idAbs == 24) {
255
mf = particleDataPtr->m0(idAbs);
256
if (mH > 2. * mf + MASSMARGIN) {
258
ps = sqrtpos(1. - 4. * mr);
259
ZpSum += pow2(coupZpWW * cos2tW) * pow3(ps)
260
* (1. + 20. * mr + 12. * mr*mr)
261
* particleDataPtr->resOpenFrac(24, -24);
266
// Calculate prefactors for gamma/Z0/Z'0 cross section terms.
267
double propZ = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
268
double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
269
gamNorm = 4. * M_PI * pow2(alpEM) / (3. * sH);
270
gamZNorm = gamNorm * 2. * thetaWRat * (sH - m2Z) * propZ;
271
ZNorm = gamNorm * pow2(thetaWRat) * sH * propZ;
272
gamZpNorm = gamNorm * 2. * thetaWRat * (sH - m2Res) * propZp;
273
ZZpNorm = gamNorm * 2. * pow2(thetaWRat) * ((sH - m2Z) * (sH - m2Res)
274
+ sH * GamMRatZ * sH * GamMRat) * propZ * propZp;
275
ZpNorm = gamNorm * pow2(thetaWRat) * sH * propZp;
277
// Optionally only keep some of gamma*, Z0 and Z' terms.
278
if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.;
279
ZZpNorm = 0.; ZpNorm = 0.;}
280
if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;
281
ZZpNorm = 0.; ZpNorm = 0.;}
282
if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
283
gamZpNorm = 0.; ZZpNorm = 0.;}
284
if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
285
if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
286
if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
290
//--------------------------------------------------------------------------
292
// Evaluate sigmaHat(sHat), including incoming flavour dependence.
294
double Sigma1ffbar2gmZZprime::sigmaHat() {
296
// Couplings to an incoming flavour.
297
int idAbs = abs(id1);
298
double ei = couplingsPtr->ef(idAbs);
299
double ai = couplingsPtr->af(idAbs);
300
double vi = couplingsPtr->vf(idAbs);
301
double api = afZp[idAbs];
302
double vpi = vfZp[idAbs];
303
double ei2 = ei * ei;
304
double eivi = ei * vi;
305
double vai2 = vi * vi + ai * ai;
306
double eivpi = ei * vpi;
307
double vaivapi = vi * vpi + ai * api;;
308
double vapi2 = vpi * vpi + api * api;
310
// Combine gamma, interference and Z0 parts.
311
double sigma = ei2 * gamNorm * gamSum + eivi * gamZNorm * gamZSum
312
+ vai2 * ZNorm * ZSum + eivpi * gamZpNorm * gamZpSum
313
+ vaivapi * ZZpNorm * ZZpSum + vapi2 * ZpNorm * ZpSum;
315
// Colour factor. Answer.
316
if (idAbs < 9) sigma /= 3.;
321
//--------------------------------------------------------------------------
323
// Select identity, colour and anticolour.
325
void Sigma1ffbar2gmZZprime::setIdColAcol() {
328
setId( id1, id2, 32);
330
// Colour flow topologies. Swap when antiquarks.
331
if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
332
else setColAcol( 0, 0, 0, 0, 0, 0);
333
if (id1 < 0) swapColAcol();
337
//--------------------------------------------------------------------------
339
// Evaluate weight for gamma*/Z0/Z'0 decay angle.
341
double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg,
344
// Default values, in- and out-flavours in process.
347
int idInAbs = process[3].idAbs();
348
int idOutAbs = process[6].idAbs();
350
// Angular weight for outgoing fermion pair.
351
if (iResBeg == 5 && iResEnd == 5 &&
352
(idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
354
// Couplings for in- and out-flavours.
355
double ei = couplingsPtr->ef(idInAbs);
356
double vi = couplingsPtr->vf(idInAbs);
357
double ai = couplingsPtr->af(idInAbs);
358
double vpi = vfZp[idInAbs];
359
double api = afZp[idInAbs];
360
double ef = couplingsPtr->ef(idOutAbs);
361
double vf = couplingsPtr->vf(idOutAbs);
362
double af = couplingsPtr->af(idOutAbs);
363
double vpf = vfZp[idOutAbs];
364
double apf = afZp[idOutAbs];
366
// Phase space factors. (One power of beta left out in formulae.)
367
double mr1 = pow2(process[6].m()) / sH;
368
double mr2 = pow2(process[7].m()) / sH;
369
double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
370
double mrAvg = 0.5 * (mr1 + mr2) - 0.25 * pow2(mr1 - mr2);
372
// Coefficients of angular expression.
373
double coefTran = ei*ei * gamNorm * ef*ef + ei * vi * gamZNorm * ef * vf
374
+ (vi*vi + ai*ai) * ZNorm * (vf*vf + ps*ps * af*af)
375
+ ei * vpi * gamZpNorm * ef * vpf
376
+ (vi * vpi + ai * api) * ZZpNorm * (vf * vpf + ps*ps * af * apf)
377
+ (vpi*vpi + api*api) * ZpNorm * (vpf*vpf + ps*ps * apf*apf);
378
double coefLong = 4. * mrAvg * ( ei*ei * gamNorm * ef*ef
379
+ ei * vi * gamZNorm * ef * vf + (vi*vi + ai*ai) * ZNorm * vf*vf
380
+ ei * vpi * gamZpNorm * ef * vpf
381
+ (vi * vpi + ai * api) * ZZpNorm * vf * vpf
382
+ (vpi*vpi + api*api) * ZpNorm * vpf*vpf );
383
double coefAsym = ps * ( ei * ai * gamZNorm * ef * af
384
+ 4. * vi * ai * ZNorm * vf * af + ei * api * gamZpNorm * ef * apf
385
+ (vi * api + vpi * ai) * ZZpNorm * (vf * apf + vpf * af)
386
+ 4. * vpi * api * ZpNorm * vpf * apf );
388
// Flip asymmetry for in-fermion + out-antifermion.
389
if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
391
// Reconstruct decay angle and weight for it.
392
double cosThe = (process[3].p() - process[4].p())
393
* (process[7].p() - process[6].p()) / (sH * ps);
394
wt = coefTran * (1. + pow2(cosThe))
395
+ coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
396
wtMax = 2. * (coefTran + abs(coefAsym));
399
// Angular weight for Z' -> W+ W-.
400
else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
401
double mr1 = pow2(process[6].m()) / sH;
402
double mr2 = pow2(process[7].m()) / sH;
403
double ps = sqrtpos(pow2(1. - mr1 -mr2) - 4. * mr1 * mr2);
404
double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2
405
+ mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
406
double cFlat = -cCos2 + 0.5 * (mr1 + mr2)
407
* (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
409
// Reconstruct decay angle and weight for it.
410
double cosThe = (process[3].p() - process[4].p())
411
* (process[7].p() - process[6].p()) / (sH * ps);
412
wt = cFlat + cCos2 * cosThe*cosThe;
413
wtMax = cFlat + max(0., cCos2);
416
// Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions.
417
else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) {
419
// Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
420
// with f' fbar' from W- and f" fbar" from W+.
421
int i1 = (process[3].id() < 0) ? 3 : 4;
423
int i3 = (process[8].id() > 0) ? 8 : 9;
425
int i5 = (process[10].id() > 0) ? 10 : 11;
427
if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);}
429
// Decay distribution like in f fbar -> Z^* -> W+ W-.
430
if (rndmPtr->flat() > anglesZpWW) {
432
// Set up four-products and internal products.
433
setupProd( process, i1, i2, i3, i4, i5, i6);
435
// tHat and uHat of fbar f -> W- W+, and their squared masses.
436
int iNeg = (process[6].id() < 0) ? 6 : 7;
437
int iPos = 13 - iNeg;
438
double tHres = (process[i1].p() - process[iNeg].p()).m2Calc();
439
double uHres = (process[i1].p() - process[iPos].p()).m2Calc();
440
double s3now = process[iNeg].m2();
441
double s4now = process[iPos].m2();
443
// Kinematics combinations (norm(x) = |x|^2).
444
double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
445
double fGK253 = norm(fGK( 2, 1, 5, 6, 3, 4) - fGK( 2, 1, 3, 4, 5, 6) );
446
double xiT = xiGK( tHres, uHres, s3now, s4now);
447
double xiU = xiGK( uHres, tHres, s3now, s4now);
448
double xjTU = xjGK( tHres, uHres, s3now, s4now);
450
// Couplings of incoming (anti)fermion. Combine with kinematics.
451
int idAbs = process[i1].idAbs();
452
double li = 0.5 * (vfZp[idAbs] + afZp[idAbs]);
453
double ri = 0.5 * (vfZp[idAbs] - afZp[idAbs]);
454
wt = li*li * fGK135 + ri*ri * fGK253;
455
wtMax = 4. * s3now * s4now * (li*li + ri*ri)
456
* (xiT + xiU - xjTU);
458
// Decay distribution like in f fbar -> h^0 -> W+ W-.
460
double p35 = 2. * process[i3].p() * process[i5].p();
461
double p46 = 2. * process[i4].p() * process[i6].p();
462
wt = 16. * p35 * p46;
467
// Angular weight in top decay by standard routine.
468
else if (process[process[iResBeg].mother1()].idAbs() == 6)
469
return weightTopDecay( process, iResBeg, iResEnd);
477
//==========================================================================
479
// Sigma1ffbar2Wprime class.
480
// Cross section for f fbar' -> W'+- (f is quark or lepton).
482
//--------------------------------------------------------------------------
484
// Initialize process.
486
void Sigma1ffbar2Wprime::initProc() {
488
// Store W+- mass and width for propagator.
489
mRes = particleDataPtr->m0(34);
490
GammaRes = particleDataPtr->mWidth(34);
492
GamMRat = GammaRes / mRes;
493
thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
495
// Axial and vector couplings of fermions.
496
aqWp = settingsPtr->parm("Wprime:aq");
497
vqWp = settingsPtr->parm("Wprime:vq");
498
alWp = settingsPtr->parm("Wprime:al");
499
vlWp = settingsPtr->parm("Wprime:vl");
501
// Coupling for W' -> W Z and decay angular admixture.
502
coupWpWZ = settingsPtr->parm("Wprime:coup2WZ");
503
anglesWpWZ = settingsPtr->parm("Wprime:anglesWZ");
505
// Set pointer to particle properties and decay table.
506
particlePtr = particleDataPtr->particleDataEntryPtr(34);
510
//--------------------------------------------------------------------------
512
// Evaluate sigmaHat(sHat), part independent of incoming flavour.
514
void Sigma1ffbar2Wprime::sigmaKin() {
516
// Set up Breit-Wigner. Cross section for W+ and W- separately.
517
double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
518
double preFac = alpEM * thetaWRat * mH;
519
sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(34, mH);
520
sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-34, mH);
524
//--------------------------------------------------------------------------
526
// Evaluate sigmaHat(sHat), including incoming flavour dependence.
528
double Sigma1ffbar2Wprime::sigmaHat() {
530
// Secondary width for W+ or W-. CKM and colour factors.
531
int idUp = (abs(id1)%2 == 0) ? id1 : id2;
532
double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
533
if (abs(id1) < 7) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
536
if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp);
537
else sigma *= 0.5 * (alWp * alWp + vlWp * vlWp);
544
//--------------------------------------------------------------------------
546
// Select identity, colour and anticolour.
548
void Sigma1ffbar2Wprime::setIdColAcol() {
550
// Sign of outgoing W.
551
int sign = 1 - 2 * (abs(id1)%2);
552
if (id1 < 0) sign = -sign;
553
setId( id1, id2, 34 * sign);
555
// Colour flow topologies. Swap when antiquarks.
556
if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
557
else setColAcol( 0, 0, 0, 0, 0, 0);
558
if (id1 < 0) swapColAcol();
562
//--------------------------------------------------------------------------
564
// Evaluate weight for W decay angle.
566
double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg,
569
// Default values, in- and out-flavours in process.
572
int idInAbs = process[3].idAbs();
573
int idOutAbs = process[6].idAbs();
575
// Angular weight for outgoing fermion pair.
576
if (iResBeg == 5 && iResEnd == 5 &&
577
(idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
579
// Couplings for in- and out-flavours.
580
double ai = (idInAbs < 9) ? aqWp : alWp;
581
double vi = (idInAbs < 9) ? vqWp : vlWp;
582
double af = (idOutAbs < 9) ? aqWp : alWp;
583
double vf = (idOutAbs < 9) ? vqWp : vlWp;
585
// Asymmetry expression.
586
double coefAsym = 8. * vi * ai * vf * af
587
/ ((vi*vi + ai*ai) * (vf*vf + af*af));
589
// Flip asymmetry for in-fermion + out-antifermion.
590
if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
592
// Phase space factors.
593
double mr1 = pow2(process[6].m()) / sH;
594
double mr2 = pow2(process[7].m()) / sH;
595
double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
597
// Reconstruct decay angle and weight for it.
598
double cosThe = (process[3].p() - process[4].p())
599
* (process[7].p() - process[6].p()) / (sH * ps);
600
wt = 1. + coefAsym * cosThe + cosThe * cosThe;
601
wtMax = 2. + abs(coefAsym);
604
// Angular weight for W' -> W Z.
605
else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
606
double mr1 = pow2(process[6].m()) / sH;
607
double mr2 = pow2(process[7].m()) / sH;
608
double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
609
double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2
610
+ mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
611
double cFlat = -cCos2 + 0.5 * (mr1 + mr2)
612
* (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
614
// Reconstruct decay angle and weight for it.
615
double cosThe = (process[3].p() - process[4].p())
616
* (process[7].p() - process[6].p()) / (sH * ps);
617
wt = cFlat + cCos2 * cosThe*cosThe;
618
wtMax = cFlat + max(0., cCos2);
621
// Angular weight for f + fbar -> W' -> W + Z -> 4 fermions.
622
else if (iResBeg == 6 && iResEnd == 7
623
&& (idOutAbs == 24 || idOutAbs == 23)) {
625
// Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
626
// with f' fbar' from W and f" fbar" from Z.
627
int i1 = (process[3].id() < 0) ? 3 : 4;
629
int i3 = (process[8].id() > 0) ? 8 : 9;
631
int i5 = (process[10].id() > 0) ? 10 : 11;
633
if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);}
635
// Decay distribution like in f fbar -> Z^* -> W+ W-.
636
if (rndmPtr->flat() > anglesWpWZ) {
638
// Set up four-products and internal products.
639
setupProd( process, i1, i2, i3, i4, i5, i6);
641
// tHat and uHat of fbar f -> W Z, and their squared masses.
642
int iW = (process[6].id() == 23) ? 7 : 6;
644
double tHres = (process[i1].p() - process[iW].p()).m2Calc();
645
double uHres = (process[i1].p() - process[iZ].p()).m2Calc();
646
double s3now = process[iW].m2();
647
double s4now = process[iZ].m2();
649
// Kinematics combinations (norm(x) = |x|^2).
650
double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
651
double fGK136 = norm(fGK( 1, 2, 3, 4, 6, 5) - fGK( 1, 2, 6, 5, 3, 4) );
652
double xiT = xiGK( tHres, uHres, s3now, s4now);
653
double xiU = xiGK( uHres, tHres, s3now, s4now);
654
double xjTU = xjGK( tHres, uHres, s3now, s4now);
656
// Couplings of outgoing fermion from Z. Combine with kinematics.
657
int idAbs = process[i5].idAbs();
658
double lfZ = couplingsPtr->lf(idAbs);
659
double rfZ = couplingsPtr->rf(idAbs);
660
wt = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136;
661
wtMax = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ)
662
* (xiT + xiU - xjTU);
664
// Decay distribution like in f fbar -> H^+- -> W+- Z0.
666
double p35 = 2. * process[i3].p() * process[i5].p();
667
double p46 = 2. * process[i4].p() * process[i6].p();
668
wt = 16. * p35 * p46;
673
// Angular weight in top decay by standard routine.
674
else if (process[process[iResBeg].mother1()].idAbs() == 6)
675
return weightTopDecay( process, iResBeg, iResEnd);
683
//==========================================================================
685
// Sigma1ffbar2Rhorizontal class.
686
// Cross section for f fbar' -> R^0 (f is a quark or lepton).
688
//--------------------------------------------------------------------------
690
// Initialize process.
692
void Sigma1ffbar2Rhorizontal::initProc() {
694
// Store R^0 mass and width for propagator.
695
mRes = particleDataPtr->m0(41);
696
GammaRes = particleDataPtr->mWidth(41);
698
GamMRat = GammaRes / mRes;
699
thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
701
// Set pointer to particle properties and decay table.
702
particlePtr = particleDataPtr->particleDataEntryPtr(41);
706
//--------------------------------------------------------------------------
708
// Evaluate sigmaHat(sHat), part independent of incoming flavour.
710
void Sigma1ffbar2Rhorizontal::sigmaKin() {
712
// Set up Breit-Wigner. Cross section for W+ and W- separately.
713
double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
714
double preFac = alpEM * thetaWRat * mH;
715
sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(41, mH);
716
sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-41, mH);
720
//--------------------------------------------------------------------------
722
// Evaluate sigmaHat(sHat), including incoming flavour dependence.
724
double Sigma1ffbar2Rhorizontal::sigmaHat() {
726
// Check for allowed flavour combinations, one generation apart.
727
if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.;
729
// Find whether R0 or R0bar. Colour factors.
730
double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg;
731
if (abs(id1) < 7) sigma /= 3.;
738
//--------------------------------------------------------------------------
740
// Select identity, colour and anticolour.
742
void Sigma1ffbar2Rhorizontal::setIdColAcol() {
744
// Outgoing R0 or R0bar.
745
id3 = (id1 +id2 > 0) ? 41 : -41;
746
setId( id1, id2, id3);
748
// Colour flow topologies. Swap when antiquarks.
749
if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
750
else setColAcol( 0, 0, 0, 0, 0, 0);
751
if (id1 < 0) swapColAcol();
755
//==========================================================================
757
} // end namespace Pythia8
added
removed
src/SigmaNewGaugeBosons.cc
