]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PYTHIA8/pythia8130/src/SigmaNewGaugeBosons.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaNewGaugeBosons.cxx
... / ...
CommitLineData
1// SigmaNewGaugeBosons.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for the
7// leptoquark simulation classes.
8
9#include "SigmaNewGaugeBosons.h"
10
11namespace Pythia8 {
12
13
14//**************************************************************************
15
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.
19
20//*********
21
22// Calculate and store internal products.
23
24void Sigma1ffbarZprimeWprime::setupProd( Event& process, int i1, int i2,
25 int i3, int i4, int i5, int i6) {
26
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();
34
35 // Do random rotation to avoid accidental zeroes in HA expressions.
36 bool smallPT = false;
37 do {
38 smallPT = false;
39 double thetaNow = acos(2. * Rndm::flat() - 1.);
40 double phiNow = 2. * M_PI * Rndm::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;
44 }
45 } while (smallPT);
46
47 // Calculate internal products.
48 for (int i = 1; i < 6; ++i) {
49 for (int j = i + 1; j <= 6; ++j) {
50 hA[i][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] );
56 if (i <= 2) {
57 hA[i][j] *= complex( 0., 1.);
58 hC[i][j] *= complex( 0., 1.);
59 }
60 hA[j][i] = - hA[i][j];
61 hC[j][i] = - hC[i][j];
62 }
63 }
64
65}
66
67//*********
68
69// Evaluate the F function of Gunion and Kunszt.
70
71complex Sigma1ffbarZprimeWprime::fGK(int j1, int j2, int j3, int j4,
72 int j5, int j6) {
73
74 return 4. * hA[j1][j3] * hC[j2][j6]
75 * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
76
77}
78
79//*********
80
81// Evaluate the Xi function of Gunion and Kunszt.
82
83double Sigma1ffbarZprimeWprime::xiGK( double tHnow, double uHnow,
84 double s3now, double s4now) {
85
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) );
90
91}
92
93//*********
94
95// Evaluate the Xj function of Gunion and Kunszt.
96
97double Sigma1ffbarZprimeWprime::xjGK( double tHnow, double uHnow,
98 double s3now, double s4now) {
99
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) );
104
105}
106
107//*************************************************************************
108
109// Sigma1ffbar2gmZZprime class.
110// Cross section for f fbar -> gamma*/Z0/Z'0 (f is quark or lepton).
111
112//*********
113
114// Initialize process.
115
116void Sigma1ffbar2gmZZprime::initProc() {
117
118 // Allow to pick only parts of full gamma*/Z0/Z'0 expression.
119 gmZmode = Settings::mode("Zprime:gmZmode");
120
121 // Store Z'0 mass and width for propagator.
122 mRes = ParticleDataTable::m0(32);
123 GammaRes = ParticleDataTable::mWidth(32);
124 m2Res = mRes*mRes;
125 GamMRat = GammaRes / mRes;
126 sin2tW = CoupEW::sin2thetaW();
127 cos2tW = 1. - sin2tW;
128 thetaWRat = 1. / (16. * sin2tW * cos2tW);
129
130 // Properties of Z0 resonance also needed.
131 mZ = ParticleDataTable::m0(23);
132 GammaZ = ParticleDataTable::mWidth(23);
133 m2Z = mZ*mZ;
134 GamMRatZ = GammaZ / mZ;
135
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.;
139
140 // Store first-generation axial and vector couplings.
141 afZp[1] = Settings::parm("Zprime:ad");
142 afZp[2] = Settings::parm("Zprime:au");
143 afZp[11] = Settings::parm("Zprime:ae");
144 afZp[12] = Settings::parm("Zprime:anue");
145 vfZp[1] = Settings::parm("Zprime:vd");
146 vfZp[2] = Settings::parm("Zprime:vu");
147 vfZp[11] = Settings::parm("Zprime:ve");
148 vfZp[12] = Settings::parm("Zprime:vnue");
149
150 // Second and third generation could be carbon copy of this...
151 if (Settings::flag("Zprime:universality")) {
152 for (int i = 3; i <= 6; ++i) {
153 afZp[i] = afZp[i-2];
154 vfZp[i] = vfZp[i-2];
155 afZp[i+10] = afZp[i+8];
156 vfZp[i+10] = vfZp[i+8];
157 }
158
159 // ... or could have different couplings.
160 } else {
161 afZp[3] = Settings::parm("Zprime:as");
162 afZp[4] = Settings::parm("Zprime:ac");
163 afZp[5] = Settings::parm("Zprime:ab");
164 afZp[6] = Settings::parm("Zprime:at");
165 afZp[13] = Settings::parm("Zprime:amu");
166 afZp[14] = Settings::parm("Zprime:anumu");
167 afZp[15] = Settings::parm("Zprime:atau");
168 afZp[16] = Settings::parm("Zprime:anutau");
169 vfZp[3] = Settings::parm("Zprime:vs");
170 vfZp[4] = Settings::parm("Zprime:vc");
171 vfZp[5] = Settings::parm("Zprime:vb");
172 vfZp[6] = Settings::parm("Zprime:vt");
173 vfZp[13] = Settings::parm("Zprime:vmu");
174 vfZp[14] = Settings::parm("Zprime:vnumu");
175 vfZp[15] = Settings::parm("Zprime:vtau");
176 vfZp[16] = Settings::parm("Zprime:vnutau");
177 }
178
179 // Coupling for Z' -> W+ W- and decay angular admixture.
180 coupZpWW = Settings::parm("Zprime:coup2WW");
181 anglesZpWW = Settings::parm("Zprime:anglesWW");
182
183 // Set pointer to particle properties and decay table.
184 particlePtr = ParticleDataTable::particleDataPtr(32);
185
186}
187
188//*********
189
190// Evaluate sigmaHat(sHat), part independent of incoming flavour.
191
192void Sigma1ffbar2gmZZprime::sigmaKin() {
193
194 // Common coupling factors.
195 double colQ = 3. * (1. + alpS / M_PI);
196
197 // Reset quantities to sum. Declare variables in loop.
198 gamSum = 0.;
199 gamZSum = 0.;
200 ZSum = 0.;
201 gamZpSum = 0.;
202 ZZpSum = 0.;
203 ZpSum = 0.;
204 int idAbs, onMode;
205 double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf,
206 ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf;
207
208 // Loop over all open Z'0 decay channels.
209 for (int i = 0; i < particlePtr->decay.size(); ++i) {
210 onMode = particlePtr->decay[i].onMode();
211 if (onMode != 1 && onMode != 2) continue;
212 idAbs = abs( particlePtr->decay[i].product(0) );
213
214 // Contributions from three fermion generations.
215 if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) {
216 mf = ParticleDataTable::m0(idAbs);
217
218 // Check that above threshold.
219 if (mH > 2. * mf + MASSMARGIN) {
220 mr = pow2(mf / mH);
221 ps = sqrtpos(1. - 4. * mr);
222
223 // Couplings of gamma^*/Z^0/Z'^0 to final flavour
224 ef = CoupEW::ef(idAbs);
225 af = CoupEW::af(idAbs);
226 vf = CoupEW::vf(idAbs);
227 apf = afZp[idAbs];
228 vpf = vfZp[idAbs];
229
230 // Combine couplings with kinematical factors.
231 kinFacA = pow3(ps);
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;
239
240 // Colour factor. Additionally secondary width for top.
241 colf = (idAbs < 9) ? colQ : 1.;
242 if (idAbs == 6) colf *= ParticleDataTable::resOpenFrac(6, -6);
243
244 // Store sum of combinations.
245 gamSum += colf * ef2;
246 gamZSum += colf * efvf;
247 ZSum += colf * vaf2;
248 gamZpSum += colf * efvpf;
249 ZZpSum += colf * vafvapf;
250 ZpSum += colf * vapf2;
251 }
252
253 // Optional contribution from W+ W-.
254 } else if (idAbs == 24) {
255 mf = ParticleDataTable::m0(idAbs);
256 if (mH > 2. * mf + MASSMARGIN) {
257 mr = pow2(mf / mH);
258 ps = sqrtpos(1. - 4. * mr);
259 ZpSum += pow2(coupZpWW * cos2tW) * pow3(ps)
260 * (1. + 20. * mr + 12. * mr*mr)
261 * ParticleDataTable::resOpenFrac(24, -24);
262 }
263 }
264 }
265
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;
276
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.;}
287
288}
289
290//*********
291
292// Evaluate sigmaHat(sHat), including incoming flavour dependence.
293
294double Sigma1ffbar2gmZZprime::sigmaHat() {
295
296 // Couplings to an incoming flavour.
297 int idAbs = abs(id1);
298 double ei = CoupEW::ef(idAbs);
299 double ai = CoupEW::af(idAbs);
300 double vi = CoupEW::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;
309
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;
314
315 // Colour factor. Answer.
316 if (idAbs < 9) sigma /= 3.;
317 return sigma;
318
319}
320
321//*********
322
323// Select identity, colour and anticolour.
324
325void Sigma1ffbar2gmZZprime::setIdColAcol() {
326
327 // Flavours trivial.
328 setId( id1, id2, 32);
329
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();
334
335}
336
337//*********
338
339// Evaluate weight for gamma*/Z0/Z'0 decay angle.
340
341double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg,
342 int iResEnd) {
343
344 // Default values, in- and out-flavours in process.
345 double wt = 1.;
346 double wtMax = 1.;
347 int idInAbs = process[3].idAbs();
348 int idOutAbs = process[6].idAbs();
349
350 // Angular weight for outgoing fermion pair.
351 if (iResBeg == 5 && iResEnd == 5 &&
352 (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
353
354 // Couplings for in- and out-flavours.
355 double ei = CoupEW::ef(idInAbs);
356 double vi = CoupEW::vf(idInAbs);
357 double ai = CoupEW::af(idInAbs);
358 double vpi = vfZp[idInAbs];
359 double api = afZp[idInAbs];
360 double ef = CoupEW::ef(idOutAbs);
361 double vf = CoupEW::vf(idOutAbs);
362 double af = CoupEW::af(idOutAbs);
363 double vpf = vfZp[idOutAbs];
364 double apf = afZp[idOutAbs];
365
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);
371
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 );
387
388 // Flip asymmetry for in-fermion + out-antifermion.
389 if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
390
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));
397 }
398
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));
408
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);
414 }
415
416 // Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions..
417 else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) {
418
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;
422 int i2 = 7 - i1;
423 int i3 = (process[8].id() > 0) ? 8 : 9;
424 int i4 = 17 - i3;
425 int i5 = (process[10].id() > 0) ? 10 : 11;
426 int i6 = 21 - i5;
427 if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);}
428
429 // Decay distribution like in f fbar -> Z^* -> W+ W-.
430 if (Rndm::flat() > anglesZpWW) {
431
432 // Set up four-products and internal products.
433 setupProd( process, i1, i2, i3, i4, i5, i6);
434
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();
442
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);
449
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);
457
458 // Decay distribution like in f fbar -> h^0 -> W+ W-.
459 } else {
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;
463 wtMax = sH2;
464 }
465 }
466
467 // Angular weight in top decay by standard routine.
468 else if (process[process[iResBeg].mother1()].idAbs() == 6)
469 return weightTopDecay( process, iResBeg, iResEnd);
470
471
472 // Done.
473 return (wt / wtMax);
474
475}
476
477//**************************************************************************
478
479// Sigma1ffbar2Wprime class.
480// Cross section for f fbar' -> W'+- (f is quark or lepton).
481
482//*********
483
484// Initialize process.
485
486void Sigma1ffbar2Wprime::initProc() {
487
488 // Store W+- mass and width for propagator.
489 mRes = ParticleDataTable::m0(34);
490 GammaRes = ParticleDataTable::mWidth(34);
491 m2Res = mRes*mRes;
492 GamMRat = GammaRes / mRes;
493 thetaWRat = 1. / (12. * CoupEW::sin2thetaW());
494
495 // Axial and vector couplings of fermions.
496 aqWp = Settings::parm("Wprime:aq");
497 vqWp = Settings::parm("Wprime:vq");
498 alWp = Settings::parm("Wprime:al");
499 vlWp = Settings::parm("Wprime:vl");
500
501 // Coupling for W' -> W Z and decay angular admixture.
502 coupWpWZ = Settings::parm("Wprime:coup2WZ");
503 anglesWpWZ = Settings::parm("Wprime:anglesWZ");
504
505 // Set pointer to particle properties and decay table.
506 particlePtr = ParticleDataTable::particleDataPtr(34);
507
508}
509
510//*********
511
512// Evaluate sigmaHat(sHat), part independent of incoming flavour.
513
514void Sigma1ffbar2Wprime::sigmaKin() {
515
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);
521
522}
523
524//*********
525
526// Evaluate sigmaHat(sHat), including incoming flavour dependence.
527
528double Sigma1ffbar2Wprime::sigmaHat() {
529
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 *= VCKM::V2id(abs(id1), abs(id2)) / 3.;
534
535 // Couplings.
536 if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp);
537 else sigma *= 0.5 * (alWp * alWp + vlWp * vlWp);
538
539 // Answer.
540 return sigma;
541
542}
543
544//*********
545
546// Select identity, colour and anticolour.
547
548void Sigma1ffbar2Wprime::setIdColAcol() {
549
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);
554
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();
559
560}
561
562//*********
563
564// Evaluate weight for W decay angle.
565
566double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg,
567 int iResEnd) {
568
569 // Default values, in- and out-flavours in process.
570 double wt = 1.;
571 double wtMax = 1.;
572 int idInAbs = process[3].idAbs();
573 int idOutAbs = process[6].idAbs();
574
575 // Angular weight for outgoing fermion pair.
576 if (iResBeg == 5 && iResEnd == 5 &&
577 (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
578
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;
584
585 // Asymmetry expression.
586 double coefAsym = 8. * vi * ai * vf * af
587 / ((vi*vi + ai*ai) * (vf*vf + af*af));
588
589 // Flip asymmetry for in-fermion + out-antifermion.
590 if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
591
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);
596
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);
602 }
603
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));
613
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);
619 }
620
621 // Angular weight for f + fbar -> W' -> W + Z -> 4 fermions.
622 else if (iResBeg == 6 && iResEnd == 7
623 && (idOutAbs == 24 || idOutAbs == 23)) {
624
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;
628 int i2 = 7 - i1;
629 int i3 = (process[8].id() > 0) ? 8 : 9;
630 int i4 = 17 - i3;
631 int i5 = (process[10].id() > 0) ? 10 : 11;
632 int i6 = 21 - i5;
633 if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);}
634
635 // Decay distribution like in f fbar -> Z^* -> W+ W-.
636 if (Rndm::flat() > anglesWpWZ) {
637
638 // Set up four-products and internal products.
639 setupProd( process, i1, i2, i3, i4, i5, i6);
640
641 // tHat and uHat of fbar f -> W Z, and their squared masses.
642 int iW = (process[6].id() == 23) ? 7 : 6;
643 int iZ = 13 - iW;
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();
648
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);
655
656 // Couplings of outgoing fermion from Z. Combine with kinematics.
657 int idAbs = process[i5].idAbs();
658 double lfZ = CoupEW::lf(idAbs);
659 double rfZ = CoupEW::rf(idAbs);
660 wt = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136;
661 wtMax = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ)
662 * (xiT + xiU - xjTU);
663
664 // Decay distribution like in f fbar -> H^+- -> W+- Z0.
665 } else {
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;
669 wtMax = sH2;
670 }
671 }
672
673 // Angular weight in top decay by standard routine.
674 else if (process[process[iResBeg].mother1()].idAbs() == 6)
675 return weightTopDecay( process, iResBeg, iResEnd);
676
677 // Done.
678 return (wt / wtMax);
679
680}
681
682
683//**************************************************************************
684
685// Sigma1ffbar2Rhorizontal class.
686// Cross section for f fbar' -> R^0 (f is a quark or lepton).
687
688//*********
689
690// Initialize process.
691
692void Sigma1ffbar2Rhorizontal::initProc() {
693
694 // Store R^0 mass and width for propagator.
695 mRes = ParticleDataTable::m0(41);
696 GammaRes = ParticleDataTable::mWidth(41);
697 m2Res = mRes*mRes;
698 GamMRat = GammaRes / mRes;
699 thetaWRat = 1. / (12. * CoupEW::sin2thetaW());
700
701 // Set pointer to particle properties and decay table.
702 particlePtr = ParticleDataTable::particleDataPtr(41);
703
704}
705
706//*********
707
708// Evaluate sigmaHat(sHat), part independent of incoming flavour.
709
710void Sigma1ffbar2Rhorizontal::sigmaKin() {
711
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);
717
718}
719
720//*********
721
722// Evaluate sigmaHat(sHat), including incoming flavour dependence.
723
724double Sigma1ffbar2Rhorizontal::sigmaHat() {
725
726 // Check for allowed flavour combinations, one generation apart.
727 if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.;
728
729 // Find whether R0 or R0bar. Colour factors.
730 double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg;
731 if (abs(id1) < 7) sigma /= 3.;
732
733 // Answer.
734 return sigma;
735
736}
737
738//*********
739
740// Select identity, colour and anticolour.
741
742void Sigma1ffbar2Rhorizontal::setIdColAcol() {
743
744 // Outgoing R0 or R0bar.
745 id3 = (id1 +id2 > 0) ? 41 : -41;
746 setId( id1, id2, id3);
747
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();
752
753}
754
755//**************************************************************************
756
757} // end namespace Pythia8