StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaNewGaugeBosons.cc
1 // SigmaNewGaugeBosons.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 "Pythia8/SigmaNewGaugeBosons.h"
10 
11 namespace 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 
24 void 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. * 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;
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 
71 complex 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 
83 double 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 
97 double 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 
116 void Sigma1ffbar2gmZZprime::initProc() {
117 
118  // Allow to pick only parts of full gamma*/Z0/Z'0 expression.
119  gmZmode = mode("Zprime:gmZmode");
120 
121  // Store Z'0 mass and width for propagator.
122  mRes = particleDataPtr->m0(32);
123  GammaRes = particleDataPtr->mWidth(32);
124  m2Res = mRes*mRes;
125  GamMRat = GammaRes / mRes;
126  sin2tW = coupSMPtr->sin2thetaW();
127  cos2tW = 1. - sin2tW;
128  thetaWRat = 1. / (16. * sin2tW * cos2tW);
129 
130  // Properties of Z0 resonance also needed.
131  mZ = particleDataPtr->m0(23);
132  GammaZ = particleDataPtr->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] = parm("Zprime:ad");
142  afZp[2] = parm("Zprime:au");
143  afZp[11] = parm("Zprime:ae");
144  afZp[12] = parm("Zprime:anue");
145  vfZp[1] = parm("Zprime:vd");
146  vfZp[2] = parm("Zprime:vu");
147  vfZp[11] = parm("Zprime:ve");
148  vfZp[12] = parm("Zprime:vnue");
149 
150  // Determine if the 4th generation should be included
151  bool coupZp2gen4 = flag("Zprime:coup2gen4");
152 
153  maxZpGen = (coupZp2gen4) ? 8 : 6;
154 
155  // Second and third (and possibly 4th) generation could be carbon copy
156  // of this...
157  if (flag("Zprime:universality")) {
158  for (int i = 3; i <= maxZpGen; ++i) {
159  afZp[i] = afZp[i-2];
160  vfZp[i] = vfZp[i-2];
161  afZp[i+10] = afZp[i+8];
162  vfZp[i+10] = vfZp[i+8];
163  }
164 
165  // ... or could have different couplings.
166  } else {
167  afZp[3] = parm("Zprime:as");
168  afZp[4] = parm("Zprime:ac");
169  afZp[5] = parm("Zprime:ab");
170  afZp[6] = parm("Zprime:at");
171  afZp[13] = parm("Zprime:amu");
172  afZp[14] = parm("Zprime:anumu");
173  afZp[15] = parm("Zprime:atau");
174  afZp[16] = parm("Zprime:anutau");
175  vfZp[3] = parm("Zprime:vs");
176  vfZp[4] = parm("Zprime:vc");
177  vfZp[5] = parm("Zprime:vb");
178  vfZp[6] = parm("Zprime:vt");
179  vfZp[13] = parm("Zprime:vmu");
180  vfZp[14] = parm("Zprime:vnumu");
181  vfZp[15] = parm("Zprime:vtau");
182  vfZp[16] = parm("Zprime:vnutau");
183  if( coupZp2gen4 ) {
184  afZp[7] = parm("Zprime:abPrime");
185  afZp[8] = parm("Zprime:atPrime");
186  vfZp[7] = parm("Zprime:vbPrime");
187  vfZp[8] = parm("Zprime:vtPrime");
188  afZp[17] = parm("Zprime:atauPrime");
189  afZp[18] = parm("Zprime:anutauPrime");
190  vfZp[17] = parm("Zprime:vtauPrime");
191  vfZp[18] = parm("Zprime:vnutauPrime");
192  }
193  }
194 
195  // Coupling for Z' -> W+ W- and decay angular admixture.
196  coupZpWW = parm("Zprime:coup2WW");
197  anglesZpWW = parm("Zprime:anglesWW");
198 
199  // Set pointer to particle properties and decay table.
200  particlePtr = particleDataPtr->particleDataEntryPtr(32);
201 
202 }
203 
204 //--------------------------------------------------------------------------
205 
206 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
207 
208 void Sigma1ffbar2gmZZprime::sigmaKin() {
209 
210  // Common coupling factors.
211  double colQ = 3. * (1. + alpS / M_PI);
212 
213  // Reset quantities to sum. Declare variables in loop.
214  gamSum = 0.;
215  gamZSum = 0.;
216  ZSum = 0.;
217  gamZpSum = 0.;
218  ZZpSum = 0.;
219  ZpSum = 0.;
220  int idAbs, onMode;
221  double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf,
222  ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf;
223 
224  // Loop over all open Z'0 decay channels.
225  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
226  onMode = particlePtr->channel(i).onMode();
227  if (onMode != 1 && onMode != 2) continue;
228  idAbs = abs( particlePtr->channel(i).product(0) );
229 
230  // Contributions from three/four fermion generations,
231  // and optionally also from excited fermions.
232  if ( (idAbs > 0 && idAbs <= maxZpGen)
233  || (idAbs > 10 && idAbs <= maxZpGen+10)
234  || (idAbs > 4000000 && idAbs <= 4000006)
235  || (idAbs > 4000010 && idAbs <= 4000016) ) {
236  int idAbs4 = (idAbs < 4000000) ? idAbs : idAbs - 4000000;
237  mf = particleDataPtr->m0(idAbs);
238 
239  // Check that above threshold.
240  if (mH > 2. * mf + MASSMARGIN) {
241  mr = pow2(mf / mH);
242  ps = sqrtpos(1. - 4. * mr);
243 
244  // Couplings of gamma^*/Z^0/Z'^0 to final flavour
245  ef = coupSMPtr->ef(idAbs4);
246  af = coupSMPtr->af(idAbs4);
247  vf = coupSMPtr->vf(idAbs4);
248  apf = afZp[idAbs4];
249  vpf = vfZp[idAbs4];
250 
251  // Combine couplings with kinematical factors.
252  kinFacA = pow3(ps);
253  kinFacV = ps * (1. + 2. * mr);
254  ef2 = ef * ef * kinFacV;
255  efvf = ef * vf * kinFacV;
256  vaf2 = vf * vf * kinFacV + af * af * kinFacA;
257  efvpf = ef * vpf * kinFacV;
258  vafvapf = vf * vpf * kinFacV + af * apf * kinFacA;
259  vapf2 = vpf * vpf * kinFacV + apf * apf * kinFacA;
260 
261  // Colour factor. Additionally secondary width for heavy particles.
262  colf = (idAbs4 < 9) ? colQ : 1.;
263  if ( (idAbs > 5 && idAbs < 9) || (idAbs > 17 && idAbs < 19)
264  || idAbs > 4000000)
265  colf *= particleDataPtr->resOpenFrac(idAbs, -idAbs);
266 
267  // Store sum of combinations.
268  gamSum += colf * ef2;
269  gamZSum += colf * efvf;
270  ZSum += colf * vaf2;
271  gamZpSum += colf * efvpf;
272  ZZpSum += colf * vafvapf;
273  ZpSum += colf * vapf2;
274  }
275 
276  // Optional contribution from W+ W-.
277  } else if (idAbs == 24) {
278  mf = particleDataPtr->m0(idAbs);
279  if (mH > 2. * mf + MASSMARGIN) {
280  mr = pow2(mf / mH);
281  ps = sqrtpos(1. - 4. * mr);
282  ZpSum += pow2(coupZpWW * cos2tW) * pow3(ps)
283  * (1. + 20. * mr + 12. * mr*mr)
284  * particleDataPtr->resOpenFrac(24, -24);
285  }
286  }
287  }
288 
289  // Calculate prefactors for gamma/Z0/Z'0 cross section terms.
290  double propZ = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
291  double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
292  gamNorm = 4. * M_PI * pow2(alpEM) / (3. * sH);
293  gamZNorm = gamNorm * 2. * thetaWRat * (sH - m2Z) * propZ;
294  ZNorm = gamNorm * pow2(thetaWRat) * sH * propZ;
295  gamZpNorm = gamNorm * 2. * thetaWRat * (sH - m2Res) * propZp;
296  ZZpNorm = gamNorm * 2. * pow2(thetaWRat) * ((sH - m2Z) * (sH - m2Res)
297  + sH * GamMRatZ * sH * GamMRat) * propZ * propZp;
298  ZpNorm = gamNorm * pow2(thetaWRat) * sH * propZp;
299 
300  // Optionally only keep some of gamma*, Z0 and Z' terms.
301  if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.;
302  ZZpNorm = 0.; ZpNorm = 0.;}
303  if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;
304  ZZpNorm = 0.; ZpNorm = 0.;}
305  if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
306  gamZpNorm = 0.; ZZpNorm = 0.;}
307  if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
308  if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
309  if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
310 
311 }
312 
313 //--------------------------------------------------------------------------
314 
315 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
316 
317 double Sigma1ffbar2gmZZprime::sigmaHat() {
318 
319  // Couplings to an incoming flavour.
320  int idAbs = abs(id1);
321  double ei = coupSMPtr->ef(idAbs);
322  double ai = coupSMPtr->af(idAbs);
323  double vi = coupSMPtr->vf(idAbs);
324  double api = afZp[idAbs];
325  double vpi = vfZp[idAbs];
326  double ei2 = ei * ei;
327  double eivi = ei * vi;
328  double vai2 = vi * vi + ai * ai;
329  double eivpi = ei * vpi;
330  double vaivapi = vi * vpi + ai * api;;
331  double vapi2 = vpi * vpi + api * api;
332 
333  // Combine gamma, interference and Z0 parts.
334  double sigma = ei2 * gamNorm * gamSum + eivi * gamZNorm * gamZSum
335  + vai2 * ZNorm * ZSum + eivpi * gamZpNorm * gamZpSum
336  + vaivapi * ZZpNorm * ZZpSum + vapi2 * ZpNorm * ZpSum;
337 
338  // Colour factor. Answer.
339  if (idAbs < 9) sigma /= 3.;
340  return sigma;
341 
342 }
343 
344 //--------------------------------------------------------------------------
345 
346 // Select identity, colour and anticolour.
347 
348 void Sigma1ffbar2gmZZprime::setIdColAcol() {
349 
350  // Flavours trivial.
351  setId( id1, id2, 32);
352 
353  // Colour flow topologies. Swap when antiquarks.
354  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
355  else setColAcol( 0, 0, 0, 0, 0, 0);
356  if (id1 < 0) swapColAcol();
357 
358 }
359 
360 //--------------------------------------------------------------------------
361 
362 // Evaluate weight for gamma*/Z0/Z'0 decay angle.
363 
364 double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg,
365  int iResEnd) {
366 
367  // Default values, in- and out-flavours in process.
368  double wt = 1.;
369  double wtMax = 1.;
370  int idInAbs = process[3].idAbs();
371  int idOutAbs = process[6].idAbs();
372 
373  // Angular weight for outgoing fermion pair.
374  if (iResBeg == 5 && iResEnd == 5 && (idOutAbs <= maxZpGen
375  || (idOutAbs > 10 && idOutAbs <= maxZpGen+10) || idOutAbs > 4000000) ) {
376 
377  // Couplings for in- and out-flavours.
378  double ei = coupSMPtr->ef(idInAbs);
379  double vi = coupSMPtr->vf(idInAbs);
380  double ai = coupSMPtr->af(idInAbs);
381  double vpi = vfZp[idInAbs];
382  double api = afZp[idInAbs];
383  int idOutAbs4 = (idOutAbs < 4000000) ? idOutAbs : idOutAbs - 4000000;
384  double ef = coupSMPtr->ef(idOutAbs4);
385  double vf = coupSMPtr->vf(idOutAbs4);
386  double af = coupSMPtr->af(idOutAbs4);
387  double vpf = vfZp[idOutAbs4];
388  double apf = afZp[idOutAbs4];
389 
390  // Phase space factors. (One power of beta left out in formulae.)
391  double mr1 = pow2(process[6].m()) / sH;
392  double mr2 = pow2(process[7].m()) / sH;
393  double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
394  double mrAvg = 0.5 * (mr1 + mr2) - 0.25 * pow2(mr1 - mr2);
395 
396  // Coefficients of angular expression.
397  double coefTran = ei*ei * gamNorm * ef*ef + ei * vi * gamZNorm * ef * vf
398  + (vi*vi + ai*ai) * ZNorm * (vf*vf + ps*ps * af*af)
399  + ei * vpi * gamZpNorm * ef * vpf
400  + (vi * vpi + ai * api) * ZZpNorm * (vf * vpf + ps*ps * af * apf)
401  + (vpi*vpi + api*api) * ZpNorm * (vpf*vpf + ps*ps * apf*apf);
402  double coefLong = 4. * mrAvg * ( ei*ei * gamNorm * ef*ef
403  + ei * vi * gamZNorm * ef * vf + (vi*vi + ai*ai) * ZNorm * vf*vf
404  + ei * vpi * gamZpNorm * ef * vpf
405  + (vi * vpi + ai * api) * ZZpNorm * vf * vpf
406  + (vpi*vpi + api*api) * ZpNorm * vpf*vpf );
407  double coefAsym = ps * ( ei * ai * gamZNorm * ef * af
408  + 4. * vi * ai * ZNorm * vf * af + ei * api * gamZpNorm * ef * apf
409  + (vi * api + vpi * ai) * ZZpNorm * (vf * apf + vpf * af)
410  + 4. * vpi * api * ZpNorm * vpf * apf );
411 
412  // Flip asymmetry for in-fermion + out-antifermion.
413  if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
414 
415  // Reconstruct decay angle and weight for it.
416  double cosThe = (process[3].p() - process[4].p())
417  * (process[7].p() - process[6].p()) / (sH * ps);
418  wt = coefTran * (1. + pow2(cosThe))
419  + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
420  wtMax = 2. * (coefTran + abs(coefAsym));
421  }
422 
423  // Angular weight for Z' -> W+ W-.
424  else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
425  double mr1 = pow2(process[6].m()) / sH;
426  double mr2 = pow2(process[7].m()) / sH;
427  double ps = sqrtpos(pow2(1. - mr1 -mr2) - 4. * mr1 * mr2);
428  double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2
429  + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
430  double cFlat = -cCos2 + 0.5 * (mr1 + mr2)
431  * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
432 
433  // Reconstruct decay angle and weight for it.
434  double cosThe = (process[3].p() - process[4].p())
435  * (process[7].p() - process[6].p()) / (sH * ps);
436  wt = cFlat + cCos2 * cosThe*cosThe;
437  wtMax = cFlat + max(0., cCos2);
438  }
439 
440  // Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions.
441  else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) {
442 
443  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
444  // with f' fbar' from W- and f" fbar" from W+.
445  int i1 = (process[3].id() < 0) ? 3 : 4;
446  int i2 = 7 - i1;
447  int i3 = (process[8].id() > 0) ? 8 : 9;
448  int i4 = 17 - i3;
449  int i5 = (process[10].id() > 0) ? 10 : 11;
450  int i6 = 21 - i5;
451  if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);}
452 
453  // Decay distribution like in f fbar -> Z^* -> W+ W-.
454  if (rndmPtr->flat() > anglesZpWW) {
455 
456  // Set up four-products and internal products.
457  setupProd( process, i1, i2, i3, i4, i5, i6);
458 
459  // tHat and uHat of fbar f -> W- W+, and their squared masses.
460  int iNeg = (process[6].id() < 0) ? 6 : 7;
461  int iPos = 13 - iNeg;
462  double tHres = (process[i1].p() - process[iNeg].p()).m2Calc();
463  double uHres = (process[i1].p() - process[iPos].p()).m2Calc();
464  double s3now = process[iNeg].m2();
465  double s4now = process[iPos].m2();
466 
467  // Kinematics combinations (norm(x) = |x|^2).
468  double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
469  double fGK253 = norm(fGK( 2, 1, 5, 6, 3, 4) - fGK( 2, 1, 3, 4, 5, 6) );
470  double xiT = xiGK( tHres, uHres, s3now, s4now);
471  double xiU = xiGK( uHres, tHres, s3now, s4now);
472  double xjTU = xjGK( tHres, uHres, s3now, s4now);
473 
474  // Couplings of incoming (anti)fermion. Combine with kinematics.
475  int idAbs = process[i1].idAbs();
476  double li = 0.5 * (vfZp[idAbs] + afZp[idAbs]);
477  double ri = 0.5 * (vfZp[idAbs] - afZp[idAbs]);
478  wt = li*li * fGK135 + ri*ri * fGK253;
479  wtMax = 4. * s3now * s4now * (li*li + ri*ri)
480  * (xiT + xiU - xjTU);
481 
482  // Decay distribution like in f fbar -> h^0 -> W+ W-.
483  } else {
484  double p35 = 2. * process[i3].p() * process[i5].p();
485  double p46 = 2. * process[i4].p() * process[i6].p();
486  wt = 16. * p35 * p46;
487  wtMax = sH2;
488  }
489  }
490 
491  // Angular weight in top decay by standard routine.
492  else if (process[process[iResBeg].mother1()].idAbs() == 6)
493  return weightTopDecay( process, iResBeg, iResEnd);
494 
495  // Angular weight for fourth generation or excited fermions not implemented.
496 
497  // Done.
498  return (wt / wtMax);
499 
500 }
501 
502 //==========================================================================
503 
504 // Sigma1ffbar2Wprime class.
505 // Cross section for f fbar' -> W'+- (f is quark or lepton).
506 
507 //--------------------------------------------------------------------------
508 
509 // Initialize process.
510 
511 void Sigma1ffbar2Wprime::initProc() {
512 
513  // Store W+- mass and width for propagator.
514  mRes = particleDataPtr->m0(34);
515  GammaRes = particleDataPtr->mWidth(34);
516  m2Res = mRes*mRes;
517  GamMRat = GammaRes / mRes;
518  thetaWRat = 1. / (12. * coupSMPtr->sin2thetaW());
519 
520  // Axial and vector couplings of fermions.
521  aqWp = parm("Wprime:aq");
522  vqWp = parm("Wprime:vq");
523  alWp = parm("Wprime:al");
524  vlWp = parm("Wprime:vl");
525 
526  // Coupling for W' -> W Z and decay angular admixture.
527  coupWpWZ = parm("Wprime:coup2WZ");
528  anglesWpWZ = parm("Wprime:anglesWZ");
529 
530  // Set pointer to particle properties and decay table.
531  particlePtr = particleDataPtr->particleDataEntryPtr(34);
532 
533 }
534 
535 //--------------------------------------------------------------------------
536 
537 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
538 
539 void Sigma1ffbar2Wprime::sigmaKin() {
540 
541  // Set up Breit-Wigner. Cross section for W+ and W- separately.
542  double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
543  double preFac = alpEM * thetaWRat * mH;
544  sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(34, mH);
545  sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-34, mH);
546 
547 }
548 
549 //--------------------------------------------------------------------------
550 
551 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
552 
553 double Sigma1ffbar2Wprime::sigmaHat() {
554 
555  // Secondary width for W+ or W-. CKM and colour factors.
556  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
557  double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
558  if (abs(id1) < 7) sigma *= coupSMPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
559 
560  // Couplings.
561  if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp);
562  else sigma *= 0.5 * (alWp * alWp + vlWp * vlWp);
563 
564  // Answer.
565  return sigma;
566 
567 }
568 
569 //--------------------------------------------------------------------------
570 
571 // Select identity, colour and anticolour.
572 
573 void Sigma1ffbar2Wprime::setIdColAcol() {
574 
575  // Sign of outgoing W.
576  int sign = 1 - 2 * (abs(id1)%2);
577  if (id1 < 0) sign = -sign;
578  setId( id1, id2, 34 * sign);
579 
580  // Colour flow topologies. Swap when antiquarks.
581  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
582  else setColAcol( 0, 0, 0, 0, 0, 0);
583  if (id1 < 0) swapColAcol();
584 
585 }
586 
587 //--------------------------------------------------------------------------
588 
589 // Evaluate weight for W decay angle.
590 
591 double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg,
592  int iResEnd) {
593 
594  // Default values, in- and out-flavours in process.
595  double wt = 1.;
596  double wtMax = 1.;
597  int idInAbs = process[3].idAbs();
598  int idOutAbs = process[6].idAbs();
599 
600  // Angular weight for outgoing fermion pair.
601  if (iResBeg == 5 && iResEnd == 5 &&
602  (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
603 
604  // Couplings for in- and out-flavours.
605  double ai = (idInAbs < 9) ? aqWp : alWp;
606  double vi = (idInAbs < 9) ? vqWp : vlWp;
607  double af = (idOutAbs < 9) ? aqWp : alWp;
608  double vf = (idOutAbs < 9) ? vqWp : vlWp;
609 
610  // Asymmetry expression.
611  double coefAsym = 8. * vi * ai * vf * af
612  / ((vi*vi + ai*ai) * (vf*vf + af*af));
613 
614  // Flip asymmetry for in-fermion + out-antifermion.
615  if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
616 
617  // Phase space factors.
618  double mr1 = pow2(process[6].m()) / sH;
619  double mr2 = pow2(process[7].m()) / sH;
620  double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
621 
622  // Reconstruct decay angle and weight for it.
623  double cosThe = (process[3].p() - process[4].p())
624  * (process[7].p() - process[6].p()) / (sH * ps);
625  wt = 1. + coefAsym * cosThe + cosThe * cosThe;
626  wtMax = 2. + abs(coefAsym);
627  }
628 
629  // Angular weight for W' -> W Z.
630  else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
631  double mr1 = pow2(process[6].m()) / sH;
632  double mr2 = pow2(process[7].m()) / sH;
633  double ps = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
634  double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2
635  + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
636  double cFlat = -cCos2 + 0.5 * (mr1 + mr2)
637  * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
638 
639  // Reconstruct decay angle and weight for it.
640  double cosThe = (process[3].p() - process[4].p())
641  * (process[7].p() - process[6].p()) / (sH * ps);
642  wt = cFlat + cCos2 * cosThe*cosThe;
643  wtMax = cFlat + max(0., cCos2);
644  }
645 
646  // Angular weight for f + fbar -> W' -> W + Z -> 4 fermions.
647  else if (iResBeg == 6 && iResEnd == 7
648  && (idOutAbs == 24 || idOutAbs == 23)) {
649 
650  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
651  // with f' fbar' from W and f" fbar" from Z.
652  int i1 = (process[3].id() < 0) ? 3 : 4;
653  int i2 = 7 - i1;
654  int i3 = (process[8].id() > 0) ? 8 : 9;
655  int i4 = 17 - i3;
656  int i5 = (process[10].id() > 0) ? 10 : 11;
657  int i6 = 21 - i5;
658  if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);}
659 
660  // Decay distribution like in f fbar -> Z^* -> W+ W-.
661  if (rndmPtr->flat() > anglesWpWZ) {
662 
663  // Set up four-products and internal products.
664  setupProd( process, i1, i2, i3, i4, i5, i6);
665 
666  // tHat and uHat of fbar f -> W Z, and their squared masses.
667  int iW = (process[6].id() == 23) ? 7 : 6;
668  int iZ = 13 - iW;
669  double tHres = (process[i1].p() - process[iW].p()).m2Calc();
670  double uHres = (process[i1].p() - process[iZ].p()).m2Calc();
671  double s3now = process[iW].m2();
672  double s4now = process[iZ].m2();
673 
674  // Kinematics combinations (norm(x) = |x|^2).
675  double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
676  double fGK136 = norm(fGK( 1, 2, 3, 4, 6, 5) - fGK( 1, 2, 6, 5, 3, 4) );
677  double xiT = xiGK( tHres, uHres, s3now, s4now);
678  double xiU = xiGK( uHres, tHres, s3now, s4now);
679  double xjTU = xjGK( tHres, uHres, s3now, s4now);
680 
681  // Couplings of outgoing fermion from Z. Combine with kinematics.
682  int idAbs = process[i5].idAbs();
683  double lfZ = coupSMPtr->lf(idAbs);
684  double rfZ = coupSMPtr->rf(idAbs);
685  wt = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136;
686  wtMax = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ)
687  * (xiT + xiU - xjTU);
688 
689  // Decay distribution like in f fbar -> H^+- -> W+- Z0.
690  } else {
691  double p35 = 2. * process[i3].p() * process[i5].p();
692  double p46 = 2. * process[i4].p() * process[i6].p();
693  wt = 16. * p35 * p46;
694  wtMax = sH2;
695  }
696  }
697 
698  // Angular weight in top decay by standard routine.
699  else if (process[process[iResBeg].mother1()].idAbs() == 6)
700  return weightTopDecay( process, iResBeg, iResEnd);
701 
702  // Done.
703  return (wt / wtMax);
704 
705 }
706 
707 
708 //==========================================================================
709 
710 // Sigma1ffbar2Rhorizontal class.
711 // Cross section for f fbar' -> R^0 (f is a quark or lepton).
712 
713 //--------------------------------------------------------------------------
714 
715 // Initialize process.
716 
717 void Sigma1ffbar2Rhorizontal::initProc() {
718 
719  // Store R^0 mass and width for propagator.
720  mRes = particleDataPtr->m0(41);
721  GammaRes = particleDataPtr->mWidth(41);
722  m2Res = mRes*mRes;
723  GamMRat = GammaRes / mRes;
724  thetaWRat = 1. / (12. * coupSMPtr->sin2thetaW());
725 
726  // Set pointer to particle properties and decay table.
727  particlePtr = particleDataPtr->particleDataEntryPtr(41);
728 
729 }
730 
731 //--------------------------------------------------------------------------
732 
733 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
734 
735 void Sigma1ffbar2Rhorizontal::sigmaKin() {
736 
737  // Set up Breit-Wigner. Cross section for W+ and W- separately.
738  double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
739  double preFac = alpEM * thetaWRat * mH;
740  sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(41, mH);
741  sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-41, mH);
742 
743 }
744 
745 //--------------------------------------------------------------------------
746 
747 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
748 
749 double Sigma1ffbar2Rhorizontal::sigmaHat() {
750 
751  // Check for allowed flavour combinations, one generation apart.
752  if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.;
753 
754  // Find whether R0 or R0bar. Colour factors.
755  double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg;
756  if (abs(id1) < 7) sigma /= 3.;
757 
758  // Answer.
759  return sigma;
760 
761 }
762 
763 //--------------------------------------------------------------------------
764 
765 // Select identity, colour and anticolour.
766 
767 void Sigma1ffbar2Rhorizontal::setIdColAcol() {
768 
769  // Outgoing R0 or R0bar.
770  id3 = (id1 +id2 > 0) ? 41 : -41;
771  setId( id1, id2, id3);
772 
773  // Colour flow topologies. Swap when antiquarks.
774  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
775  else setColAcol( 0, 0, 0, 0, 0, 0);
776  if (id1 < 0) swapColAcol();
777 
778 }
779 
780 //==========================================================================
781 
782 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5