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) 2014 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 "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 = settingsPtr->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 = couplingsPtr->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] = 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");
149 
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) {
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] = 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");
177  }
178 
179  // Coupling for Z' -> W+ W- and decay angular admixture.
180  coupZpWW = settingsPtr->parm("Zprime:coup2WW");
181  anglesZpWW = settingsPtr->parm("Zprime:anglesWW");
182 
183  // Set pointer to particle properties and decay table.
184  particlePtr = particleDataPtr->particleDataEntryPtr(32);
185 
186 }
187 
188 //--------------------------------------------------------------------------
189 
190 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
191 
192 void 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->sizeChannels(); ++i) {
210  onMode = particlePtr->channel(i).onMode();
211  if (onMode != 1 && onMode != 2) continue;
212  idAbs = abs( particlePtr->channel(i).product(0) );
213 
214  // Contributions from three fermion generations.
215  if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) {
216  mf = particleDataPtr->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 = couplingsPtr->ef(idAbs);
225  af = couplingsPtr->af(idAbs);
226  vf = couplingsPtr->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 *= particleDataPtr->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 = particleDataPtr->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  * particleDataPtr->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 
294 double Sigma1ffbar2gmZZprime::sigmaHat() {
295 
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;
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 
325 void 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 
341 double 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 = 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];
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 (rndmPtr->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 
486 void Sigma1ffbar2Wprime::initProc() {
487 
488  // Store W+- mass and width for propagator.
489  mRes = particleDataPtr->m0(34);
490  GammaRes = particleDataPtr->mWidth(34);
491  m2Res = mRes*mRes;
492  GamMRat = GammaRes / mRes;
493  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
494 
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");
500 
501  // Coupling for W' -> W Z and decay angular admixture.
502  coupWpWZ = settingsPtr->parm("Wprime:coup2WZ");
503  anglesWpWZ = settingsPtr->parm("Wprime:anglesWZ");
504 
505  // Set pointer to particle properties and decay table.
506  particlePtr = particleDataPtr->particleDataEntryPtr(34);
507 
508 }
509 
510 //--------------------------------------------------------------------------
511 
512 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
513 
514 void 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 
528 double 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 *= couplingsPtr->V2CKMid(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 
548 void 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 
566 double 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 (rndmPtr->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 = 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);
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 
692 void Sigma1ffbar2Rhorizontal::initProc() {
693 
694  // Store R^0 mass and width for propagator.
695  mRes = particleDataPtr->m0(41);
696  GammaRes = particleDataPtr->mWidth(41);
697  m2Res = mRes*mRes;
698  GamMRat = GammaRes / mRes;
699  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
700 
701  // Set pointer to particle properties and decay table.
702  particlePtr = particleDataPtr->particleDataEntryPtr(41);
703 
704 }
705 
706 //--------------------------------------------------------------------------
707 
708 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
709 
710 void 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 
724 double 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 
742 void 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
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5