StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaEW.cc
1 // SigmaEW.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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 // electroweak simulation classes (including photon processes).
8 
9 #include "Pythia8/SigmaEW.h"
10 
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // Sigma2qg2qgamma class.
16 // Cross section for q g -> q gamma.
17 
18 //--------------------------------------------------------------------------
19 
20 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
21 
22 void Sigma2qg2qgamma::sigmaKin() {
23 
24  // Calculate kinematics dependence.
25  sigUS = (1./3.) * (sH2 + uH2) / (-sH * uH);
26 
27  // Answer.
28  sigma0 = (M_PI/sH2) * alpS * alpEM * sigUS;
29 
30  }
31 
32 //--------------------------------------------------------------------------
33 
34 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
35 
36 double Sigma2qg2qgamma::sigmaHat() {
37 
38  // Incoming flavour gives charge factor.
39  int idNow = (id2 == 21) ? id1 : id2;
40  double eNow = couplingsPtr->ef( abs(idNow) );
41  return sigma0 * pow2(eNow);
42 
43 }
44 
45 //--------------------------------------------------------------------------
46 
47 // Select identity, colour and anticolour.
48 
49 void Sigma2qg2qgamma::setIdColAcol() {
50 
51  // Construct outgoing flavours.
52  id3 = (id1 == 21) ? 22 : id1;
53  id4 = (id2 == 21) ? 22 : id2;
54  setId( id1, id2, id3, id4);
55 
56  // Colour flow topology. Swap if first is gluon, or when antiquark.
57  setColAcol( 1, 0, 2, 1, 2, 0, 0, 0);
58  if (id1 == 21) swapCol1234();
59  if (id1 < 0 || id2 < 0) swapColAcol();
60 
61 }
62 
63 //==========================================================================
64 
65 // Sigma2qqbar2ggamma class.
66 // Cross section for q qbar -> g gamma.
67 
68 //--------------------------------------------------------------------------
69 
70 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
71 
72 void Sigma2qqbar2ggamma::sigmaKin() {
73 
74  // Calculate kinematics dependence.
75  double sigTU = (8./9.) * (tH2 + uH2) / (tH * uH);
76 
77  // Answer.
78  sigma0 = (M_PI/sH2) * alpS * alpEM * sigTU;
79 
80 }
81 
82 //--------------------------------------------------------------------------
83 
84 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
85 
86 double Sigma2qqbar2ggamma::sigmaHat() {
87 
88  // Incoming flavour gives charge factor.
89  double eNow = couplingsPtr->ef( abs(id1) );
90  return sigma0 * pow2(eNow);
91 
92 }
93 
94 //--------------------------------------------------------------------------
95 
96 // Select identity, colour and anticolour.
97 
98 void Sigma2qqbar2ggamma::setIdColAcol() {
99 
100  // Outgoing flavours trivial.
101  setId( id1, id2, 21, 22);
102 
103  // One colour flow topology. Swap if first is antiquark.
104  setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
105  if (id1 < 0) swapColAcol();
106 
107 }
108 
109 //==========================================================================
110 
111 // Sigma2gg2ggamma class.
112 // Cross section for g g -> g gamma.
113 // Proceeds through a quark box, by default using 5 massless quarks.
114 
115 //--------------------------------------------------------------------------
116 
117 // Initialize process, especially parton-flux object.
118 
119 void Sigma2gg2ggamma::initProc() {
120 
121  // Maximum quark flavour in loop.
122  int nQuarkLoop = settingsPtr->mode("PromptPhoton:nQuarkLoop");
123 
124  // Calculate charge factor from the allowed quarks in the box.
125  chargeSum = - 1./3. + 2./3. - 1./3.;
126  if (nQuarkLoop >= 4) chargeSum += 2./3.;
127  if (nQuarkLoop >= 5) chargeSum -= 1./3.;
128  if (nQuarkLoop >= 6) chargeSum += 2./3.;
129 
130 }
131 
132 //--------------------------------------------------------------------------
133 
134 // Evaluate d(sigmaHat)/d(tHat).
135 
136 void Sigma2gg2ggamma::sigmaKin() {
137 
138  // Logarithms of Mandelstam variable ratios.
139  double logST = log( -sH / tH );
140  double logSU = log( -sH / uH );
141  double logTU = log( tH / uH );
142 
143  // Real and imaginary parts of separate amplitudes.
144  double b0stuRe = 1. + (tH - uH) / sH * logTU
145  + 0.5 * (tH2 + uH2) / sH2 * (pow2(logTU) + pow2(M_PI));
146  double b0stuIm = 0.;
147  double b0tsuRe = 1. + (sH - uH) / tH * logSU
148  + 0.5 * (sH2 + uH2) / tH2 * pow2(logSU);
149  double b0tsuIm = -M_PI * ( (sH - uH) / tH + (sH2 + uH2) / tH2 * logSU);
150  double b0utsRe = 1. + (sH - tH) / uH * logST
151  + 0.5 * (sH2 + tH2) / uH2 * pow2(logST);
152  double b0utsIm = -M_PI * ( (sH - tH) / uH + (sH2 + tH2) / uH2 * logST);
153  double b1stuRe = -1.;
154  double b1stuIm = 0.;
155  double b2stuRe = -1.;
156  double b2stuIm = 0.;
157 
158  // Calculate kinematics dependence.
159  double sigBox = pow2(b0stuRe) + pow2(b0stuIm) + pow2(b0tsuRe)
160  + pow2(b0tsuIm) + pow2(b0utsRe) + pow2(b0utsIm) + 4. * pow2(b1stuRe)
161  + 4. * pow2(b1stuIm) + pow2(b2stuRe) + pow2(b2stuIm);
162 
163  // Answer.
164  sigma = (5. / (192. * M_PI * sH2)) * pow2(chargeSum)
165  * pow3(alpS) * alpEM * sigBox;
166 
167 }
168 
169 //--------------------------------------------------------------------------
170 
171 // Select identity, colour and anticolour.
172 
173 void Sigma2gg2ggamma::setIdColAcol() {
174 
175  // Flavours and colours are trivial.
176  setId( id1, id2, 21, 22);
177  setColAcol( 1, 2, 2, 3, 1, 3, 0, 0);
178  if (rndmPtr->flat() > 0.5) swapColAcol();
179 
180 }
181 
182 //==========================================================================
183 
184 // Sigma2ffbar2gammagamma class.
185 // Cross section for q qbar -> gamma gamma.
186 
187 //--------------------------------------------------------------------------
188 
189 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
190 
191 void Sigma2ffbar2gammagamma::sigmaKin() {
192 
193  // Calculate kinematics dependence.
194  sigTU = 2. * (tH2 + uH2) / (tH * uH);
195 
196  // Answer contains factor 1/2 from identical photons.
197  sigma0 = (M_PI/sH2) * pow2(alpEM) * 0.5 * sigTU;
198 
199 }
200 
201 //--------------------------------------------------------------------------
202 
203 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
204 
205 double Sigma2ffbar2gammagamma::sigmaHat() {
206 
207  // Incoming flavour gives charge and colour factors.
208  double eNow = couplingsPtr->ef( abs(id1) );
209  double colFac = (abs(id1) < 9) ? 1. / 3. : 1.;
210  return sigma0 * pow4(eNow) * colFac;
211 
212 }
213 
214 //--------------------------------------------------------------------------
215 
216 // Select identity, colour and anticolour.
217 
218 void Sigma2ffbar2gammagamma::setIdColAcol() {
219 
220  // Outgoing flavours trivial.
221  setId( id1, id2, 22, 22);
222 
223  // No colours at all or one flow topology. Swap if first is antiquark.
224  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
225  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
226  if (id1 < 0) swapColAcol();
227 
228 }
229 
230 //==========================================================================
231 
232 // Sigma2gg2gammagamma class.
233 // Cross section for g g -> gamma gamma.
234 // Proceeds through a quark box, by default using 5 massless quarks.
235 
236 //--------------------------------------------------------------------------
237 
238 // Initialize process.
239 
240 void Sigma2gg2gammagamma::initProc() {
241 
242  // Maximum quark flavour in loop.
243  int nQuarkLoop = settingsPtr->mode("PromptPhoton:nQuarkLoop");
244 
245  // Calculate charge factor from the allowed quarks in the box.
246  charge2Sum = 1./9. + 4./9. + 1./9.;
247  if (nQuarkLoop >= 4) charge2Sum += 4./9.;
248  if (nQuarkLoop >= 5) charge2Sum += 1./9.;
249  if (nQuarkLoop >= 6) charge2Sum += 4./9.;
250 
251 }
252 
253 //--------------------------------------------------------------------------
254 
255 // Evaluate d(sigmaHat)/d(tHat).
256 
257 void Sigma2gg2gammagamma::sigmaKin() {
258 
259  // Logarithms of Mandelstam variable ratios.
260  double logST = log( -sH / tH );
261  double logSU = log( -sH / uH );
262  double logTU = log( tH / uH );
263 
264  // Real and imaginary parts of separate amplitudes.
265  double b0stuRe = 1. + (tH - uH) / sH * logTU
266  + 0.5 * (tH2 + uH2) / sH2 * (pow2(logTU) + pow2(M_PI));
267  double b0stuIm = 0.;
268  double b0tsuRe = 1. + (sH - uH) / tH * logSU
269  + 0.5 * (sH2 + uH2) / tH2 * pow2(logSU);
270  double b0tsuIm = -M_PI * ( (sH - uH) / tH + (sH2 + uH2) / tH2 * logSU);
271  double b0utsRe = 1. + (sH - tH) / uH * logST
272  + 0.5 * (sH2 + tH2) / uH2 * pow2(logST);
273  double b0utsIm = -M_PI * ( (sH - tH) / uH + (sH2 + tH2) / uH2 * logST);
274  double b1stuRe = -1.;
275  double b1stuIm = 0.;
276  double b2stuRe = -1.;
277  double b2stuIm = 0.;
278 
279  // Calculate kinematics dependence.
280  double sigBox = pow2(b0stuRe) + pow2(b0stuIm) + pow2(b0tsuRe)
281  + pow2(b0tsuIm) + pow2(b0utsRe) + pow2(b0utsIm) + 4. * pow2(b1stuRe)
282  + 4. * pow2(b1stuIm) + pow2(b2stuRe) + pow2(b2stuIm);
283 
284  // Answer contains factor 1/2 from identical photons.
285  sigma = (0.5 / (16. * M_PI * sH2)) * pow2(charge2Sum)
286  * pow2(alpS) * pow2(alpEM) * sigBox;
287 
288 }
289 
290 //--------------------------------------------------------------------------
291 
292 // Select identity, colour and anticolour.
293 
294 void Sigma2gg2gammagamma::setIdColAcol() {
295 
296  // Flavours and colours are trivial.
297  setId( id1, id2, 22, 22);
298  setColAcol( 1, 2, 2, 1, 0, 0, 0, 0);
299 
300 }
301 
302 //==========================================================================
303 
304 // Sigma2ff2fftgmZ class.
305 // Cross section for f f' -> f f' via t-channel gamma*/Z0 exchange
306 // (f is quark or lepton).
307 
308 //--------------------------------------------------------------------------
309 
310 // Initialize process.
311 
312 void Sigma2ff2fftgmZ::initProc() {
313 
314  // Store Z0 mass for propagator. Common coupling factor.
315  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
316  mZ = particleDataPtr->m0(23);
317  mZS = mZ*mZ;
318  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
319  * couplingsPtr->cos2thetaW());
320 
321 }
322 
323 //--------------------------------------------------------------------------
324 
325 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
326 
327 void Sigma2ff2fftgmZ::sigmaKin() {
328 
329  // Cross section part common for all incoming flavours.
330  double sigma0 = (M_PI / sH2) * pow2(alpEM);
331 
332  // Kinematical functions for gamma-gamma, gamma-Z and Z-Z parts.
333  sigmagmgm = sigma0 * 2. * (sH2 + uH2) / tH2;
334  sigmagmZ = sigma0 * 4. * thetaWRat * sH2 / (tH * (tH - mZS));
335  sigmaZZ = sigma0 * 2. * pow2(thetaWRat) * sH2 / pow2(tH - mZS);
336  if (gmZmode == 1) {sigmagmZ = 0.; sigmaZZ = 0.;}
337  if (gmZmode == 2) {sigmagmgm = 0.; sigmagmZ = 0.;}
338 
339 }
340 
341 //--------------------------------------------------------------------------
342 
343 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
344 
345 double Sigma2ff2fftgmZ::sigmaHat() {
346 
347  // Couplings for current flavour combination.
348  int id1Abs = abs(id1);
349  double e1 = couplingsPtr->ef(id1Abs);
350  double v1 = couplingsPtr->vf(id1Abs);
351  double a1 = couplingsPtr->af(id1Abs);
352  int id2Abs = abs(id2);
353  double e2 = couplingsPtr->ef(id2Abs);
354  double v2 = couplingsPtr->vf(id2Abs);
355  double a2 = couplingsPtr->af(id2Abs);
356 
357  // Distinguish same-sign and opposite-sign fermions.
358  double epsi = (id1 * id2 > 0) ? 1. : -1.;
359 
360  // Flavour-dependent cross section.
361  double sigma = sigmagmgm * pow2(e1 * e2)
362  + sigmagmZ * e1 * e2 * (v1 * v2 * (1. + uH2 / sH2)
363  + a1 * a2 * epsi * (1. - uH2 / sH2))
364  + sigmaZZ * ((v1*v1 + a1*a1) * (v2*v2 + a2*a2) * (1. + uH2 / sH2)
365  + 4. * v1 * a1 * v2 * a2 * epsi * (1. - uH2 / sH2));
366 
367  // Spin-state extra factor 2 per incoming neutrino.
368  if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
369  if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
370 
371  // Answer.
372  return sigma;
373 
374 }
375 
376 //--------------------------------------------------------------------------
377 
378 // Select identity, colour and anticolour.
379 
380 void Sigma2ff2fftgmZ::setIdColAcol() {
381 
382  // Trivial flavours: out = in.
383  setId( id1, id2, id1, id2);
384 
385  // Colour flow topologies. Swap when antiquarks.
386  if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
387  setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
388  else if (abs(id1) < 9 && abs(id2) < 9)
389  setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
390  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
391  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
392  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
393  if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
394  swapColAcol();
395 
396 }
397 
398 //==========================================================================
399 
400 // Sigma2ff2fftW class.
401 // Cross section for f_1 f_2 -> f_3 f_4 via t-channel W+- exchange
402 // (f is quark or lepton).
403 
404 //--------------------------------------------------------------------------
405 
406 // Initialize process.
407 
408 void Sigma2ff2fftW::initProc() {
409 
410  // Store W+- mass for propagator. Common coupling factor.
411  mW = particleDataPtr->m0(24);
412  mWS = mW*mW;
413  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
414 
415 }
416 
417 //--------------------------------------------------------------------------
418 
419 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
420 
421 void Sigma2ff2fftW::sigmaKin() {
422 
423  // Cross section part common for all incoming flavours.
424  sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat)
425  * 4. * sH2 / pow2(tH - mWS);
426 
427 }
428 
429 //--------------------------------------------------------------------------
430 
431 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
432 
433 double Sigma2ff2fftW::sigmaHat() {
434 
435  // Some flavour combinations not possible.
436  int id1Abs = abs(id1);
437  int id2Abs = abs(id2);
438  if ( (id1Abs%2 == id2Abs%2 && id1 * id2 > 0)
439  || (id1Abs%2 != id2Abs%2 && id1 * id2 < 0) ) return 0.;
440 
441  // Basic cross section.
442  double sigma = sigma0;
443  if (id1 * id2 < 0) sigma *= uH2 / sH2;
444 
445  // CKM factors for final states.
446  sigma *= couplingsPtr->V2CKMsum(id1Abs) * couplingsPtr->V2CKMsum(id2Abs);
447 
448  // Spin-state extra factor 2 per incoming neutrino.
449  if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
450  if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
451 
452  // Answer.
453  return sigma;
454 
455 }
456 
457 //--------------------------------------------------------------------------
458 
459 // Select identity, colour and anticolour.
460 
461 void Sigma2ff2fftW::setIdColAcol() {
462 
463  // Pick out-flavours by relative CKM weights.
464  id3 = couplingsPtr->V2CKMpick(id1);
465  id4 = couplingsPtr->V2CKMpick(id2);
466  setId( id1, id2, id3, id4);
467 
468  // Colour flow topologies. Swap when antiquarks.
469  if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
470  setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
471  else if (abs(id1) < 9 && abs(id2) < 9)
472  setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
473  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
474  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
475  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
476  if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
477  swapColAcol();
478 
479 }
480 
481 
482 //==========================================================================
483 
484 // Sigma2qq2QqtW class.
485 // Cross section for q q' -> Q q" via t-channel W+- exchange.
486 // Related to Sigma2ff2ffViaW class, but with massive matrix elements.
487 
488 //--------------------------------------------------------------------------
489 
490 // Initialize process.
491 
492 void Sigma2qq2QqtW::initProc() {
493 
494  // Process name.
495  nameSave = "q q -> Q q (t-channel W+-)";
496  if (idNew == 4) nameSave = "q q -> c q (t-channel W+-)";
497  if (idNew == 5) nameSave = "q q -> b q (t-channel W+-)";
498  if (idNew == 6) nameSave = "q q -> t q (t-channel W+-)";
499  if (idNew == 7) nameSave = "q q -> b' q (t-channel W+-)";
500  if (idNew == 8) nameSave = "q q -> t' q (t-channel W+-)";
501 
502  // Store W+- mass for propagator. Common coupling factor.
503  mW = particleDataPtr->m0(24);
504  mWS = mW*mW;
505  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
506 
507  // Secondary open width fractions, relevant for top (or heavier).
508  openFracPos = particleDataPtr->resOpenFrac(idNew);
509  openFracNeg = particleDataPtr->resOpenFrac(-idNew);
510 
511 }
512 
513 //--------------------------------------------------------------------------
514 
515 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
516 
517 void Sigma2qq2QqtW::sigmaKin() {
518 
519  // Cross section part common for all incoming flavours.
520  sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat) * 4. / pow2(tH - mWS);
521 
522 }
523 
524 //--------------------------------------------------------------------------
525 
526 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
527 
528 double Sigma2qq2QqtW::sigmaHat() {
529 
530  // Some flavour combinations not possible.
531  int id1Abs = abs(id1);
532  int id2Abs = abs(id2);
533  bool diff12 = (id1Abs%2 != id2Abs%2);
534  if ( (!diff12 && id1 * id2 > 0)
535  || ( diff12 && id1 * id2 < 0) ) return 0.;
536 
537  // Basic cross section.
538  double sigma = sigma0;
539  sigma *= (id1 * id2 > 0) ? sH * (sH - s3) : uH * (uH - s3);
540 
541  // Secondary width if t or tbar produced on either side.
542  double openFrac1 = (id1 > 0) ? openFracPos : openFracNeg;
543  double openFrac2 = (id2 > 0) ? openFracPos : openFracNeg;
544 
545  // CKM factors for final states; further impossible case.
546  bool diff1N = (id1Abs%2 != idNew%2);
547  bool diff2N = (id2Abs%2 != idNew%2);
548  if (diff1N && diff2N)
549  sigma *= ( couplingsPtr->V2CKMid(id1Abs, idNew) * openFrac1
550  * couplingsPtr->V2CKMsum(id2Abs) + couplingsPtr->V2CKMsum(id1Abs)
551  * couplingsPtr->V2CKMid(id2Abs, idNew) * openFrac2 );
552  else if (diff1N)
553  sigma *= couplingsPtr->V2CKMid(id1Abs, idNew) * openFrac1
554  * couplingsPtr->V2CKMsum(id2Abs);
555  else if (diff2N)
556  sigma *= couplingsPtr->V2CKMsum(id1Abs)
557  * couplingsPtr->V2CKMid(id2Abs, idNew) * openFrac2;
558  else sigma = 0.;
559 
560  // Spin-state extra factor 2 per incoming neutrino.
561  if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
562  if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
563 
564  // Answer.
565  return sigma;
566 
567 }
568 
569 //--------------------------------------------------------------------------
570 
571 // Select identity, colour and anticolour.
572 
573 void Sigma2qq2QqtW::setIdColAcol() {
574 
575  // For topologies like d dbar -> (t/c/u) (t/c/u)bar pick side.
576  int id1Abs = abs(id1);
577  int id2Abs = abs(id2);
578  int side = 1;
579  if ( (id1Abs + idNew)%2 == 1 && (id2Abs + idNew)%2 == 1 ) {
580  double prob1 = couplingsPtr->V2CKMid(id1Abs, idNew)
581  * couplingsPtr->V2CKMsum(id2Abs);
582  prob1 *= (id1 > 0) ? openFracPos : openFracNeg;
583  double prob2 = couplingsPtr->V2CKMid(id2Abs, idNew)
584  * couplingsPtr->V2CKMsum(id1Abs);
585  prob2 *= (id2 > 0) ? openFracPos : openFracNeg;
586  if (prob2 > rndmPtr->flat() * (prob1 + prob2)) side = 2;
587  }
588  else if ((id2Abs + idNew)%2 == 1) side = 2;
589 
590  // Pick out-flavours by relative CKM weights.
591  if (side == 1) {
592  // q q' -> t q" : correct order from start.
593  id3 = (id1 > 0) ? idNew : -idNew;
594  id4 = couplingsPtr->V2CKMpick(id2);
595  setId( id1, id2, id3, id4);
596  } else {
597  // q q' -> q" t : stored as t q" so swap tHat <-> uHat.
598  swapTU = true;
599  id3 = couplingsPtr->V2CKMpick(id1);
600  id4 = (id2 > 0) ? idNew : -idNew;
601  setId( id1, id2, id4, id3);
602  }
603 
604  // Colour flow topologies. Swap when antiquarks on side 1.
605  if (side == 1 && id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
606  else if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
607  else if (side == 1) setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
608  else setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
609  if (id1 < 0) swapColAcol();
610 
611 }
612 
613 //--------------------------------------------------------------------------
614 
615 // Evaluate weight for decay angles of W in top decay.
616 
617 double Sigma2qq2QqtW::weightDecay( Event& process, int iResBeg,
618  int iResEnd) {
619 
620  // For top decay hand over to standard routine, else done.
621  if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
622  return weightTopDecay( process, iResBeg, iResEnd);
623  else return 1.;
624 
625 }
626 
627 //==========================================================================
628 
629 // Sigma1ffbar2gmZ class.
630 // Cross section for f fbar -> gamma*/Z0 (f is quark or lepton).
631 
632 //--------------------------------------------------------------------------
633 
634 // Initialize process.
635 
636 void Sigma1ffbar2gmZ::initProc() {
637 
638  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
639  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
640 
641  // Store Z0 mass and width for propagator.
642  mRes = particleDataPtr->m0(23);
643  GammaRes = particleDataPtr->mWidth(23);
644  m2Res = mRes*mRes;
645  GamMRat = GammaRes / mRes;
646  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
647  * couplingsPtr->cos2thetaW());
648 
649  // Set pointer to particle properties and decay table.
650  particlePtr = particleDataPtr->particleDataEntryPtr(23);
651 
652 }
653 
654 //--------------------------------------------------------------------------
655 
656 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
657 
658 void Sigma1ffbar2gmZ::sigmaKin() {
659 
660  // Common coupling factors.
661  double colQ = 3. * (1. + alpS / M_PI);
662 
663  // Reset quantities to sum. Declare variables in loop.
664  gamSum = 0.;
665  intSum = 0.;
666  resSum = 0.;
667  int idAbs, onMode;
668  double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
669 
670  // Loop over all Z0 decay channels.
671  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
672  idAbs = abs( particlePtr->channel(i).product(0) );
673 
674  // Only contributions from three fermion generations, except top.
675  if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
676  mf = particleDataPtr->m0(idAbs);
677 
678  // Check that above threshold. Phase space.
679  if (mH > 2. * mf + MASSMARGIN) {
680  mr = pow2(mf / mH);
681  betaf = sqrtpos(1. - 4. * mr);
682  psvec = betaf * (1. + 2. * mr);
683  psaxi = pow3(betaf);
684 
685  // Combine phase space with couplings.
686  ef2 = couplingsPtr->ef2(idAbs) * psvec;
687  efvf = couplingsPtr->efvf(idAbs) * psvec;
688  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
689  + couplingsPtr->af2(idAbs) * psaxi;
690  colf = (idAbs < 6) ? colQ : 1.;
691 
692  // Store sum of combinations. For outstate only open channels.
693  onMode = particlePtr->channel(i).onMode();
694  if (onMode == 1 || onMode == 2) {
695  gamSum += colf * ef2;
696  intSum += colf * efvf;
697  resSum += colf * vf2af2;
698  }
699 
700  // End loop over fermions.
701  }
702  }
703  }
704 
705  // Calculate prefactors for gamma/interference/Z0 cross section terms.
706  gamProp = 4. * M_PI * pow2(alpEM) / (3. * sH);
707  intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
708  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
709  resProp = gamProp * pow2(thetaWRat * sH)
710  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
711 
712  // Optionally only keep gamma* or Z0 term.
713  if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
714  if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
715 
716 }
717 
718 //--------------------------------------------------------------------------
719 
720 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
721 
722 double Sigma1ffbar2gmZ::sigmaHat() {
723 
724  // Combine gamma, interference and Z0 parts.
725  int idAbs = abs(id1);
726  double sigma = couplingsPtr->ef2(idAbs) * gamProp * gamSum
727  + couplingsPtr->efvf(idAbs) * intProp * intSum
728  + couplingsPtr->vf2af2(idAbs) * resProp * resSum;
729 
730  // Colour factor. Answer.
731  if (idAbs < 9) sigma /= 3.;
732  return sigma;
733 
734 }
735 
736 //--------------------------------------------------------------------------
737 
738 // Select identity, colour and anticolour.
739 
740 void Sigma1ffbar2gmZ::setIdColAcol() {
741 
742  // Flavours trivial.
743  setId( id1, id2, 23);
744 
745  // Colour flow topologies. Swap when antiquarks.
746  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
747  else setColAcol( 0, 0, 0, 0, 0, 0);
748  if (id1 < 0) swapColAcol();
749 
750 }
751 
752 //--------------------------------------------------------------------------
753 
754 // Evaluate weight for gamma*/Z0 decay angle.
755 
756 double Sigma1ffbar2gmZ::weightDecay( Event& process, int iResBeg,
757  int iResEnd) {
758 
759  // Z should sit in entry 5.
760  if (iResBeg != 5 || iResEnd != 5) return 1.;
761 
762  // Couplings for in- and out-flavours.
763  int idInAbs = process[3].idAbs();
764  double ei = couplingsPtr->ef(idInAbs);
765  double vi = couplingsPtr->vf(idInAbs);
766  double ai = couplingsPtr->af(idInAbs);
767  int idOutAbs = process[6].idAbs();
768  double ef = couplingsPtr->ef(idOutAbs);
769  double vf = couplingsPtr->vf(idOutAbs);
770  double af = couplingsPtr->af(idOutAbs);
771 
772  // Phase space factors. (One power of beta left out in formulae.)
773  double mf = process[6].m();
774  double mr = mf*mf / sH;
775  double betaf = sqrtpos(1. - 4. * mr);
776 
777  // Coefficients of angular expression.
778  double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
779  + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
780  double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
781  + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
782  double coefAsym = betaf * ( ei * ai * intProp * ef * af
783  + 4. * vi * ai * resProp * vf * af );
784 
785  // Flip asymmetry for in-fermion + out-antifermion.
786  if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
787 
788  // Reconstruct decay angle and weight for it.
789  double cosThe = (process[3].p() - process[4].p())
790  * (process[7].p() - process[6].p()) / (sH * betaf);
791  double wtMax = 2. * (coefTran + abs(coefAsym));
792  double wt = coefTran * (1. + pow2(cosThe))
793  + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
794 
795  // Done.
796  return (wt / wtMax);
797 
798 }
799 
800 //==========================================================================
801 
802 // Sigma1ffbar2W class.
803 // Cross section for f fbar' -> W+- (f is quark or lepton).
804 
805 //--------------------------------------------------------------------------
806 
807 // Initialize process.
808 
809 void Sigma1ffbar2W::initProc() {
810 
811  // Store W+- mass and width for propagator.
812  mRes = particleDataPtr->m0(24);
813  GammaRes = particleDataPtr->mWidth(24);
814  m2Res = mRes*mRes;
815  GamMRat = GammaRes / mRes;
816  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
817 
818  // Set pointer to particle properties and decay table.
819  particlePtr = particleDataPtr->particleDataEntryPtr(24);
820 
821 }
822 
823 //--------------------------------------------------------------------------
824 
825 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
826 
827 void Sigma1ffbar2W::sigmaKin() {
828 
829  // Set up Breit-Wigner. Cross section for W+ and W- separately.
830  double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
831  double preFac = alpEM * thetaWRat * mH;
832  sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(24, mH);
833  sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-24, mH);
834 
835 }
836 
837 //--------------------------------------------------------------------------
838 
839 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
840 
841 double Sigma1ffbar2W::sigmaHat() {
842 
843  // Secondary width for W+ or W-. CKM and colour factors.
844  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
845  double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
846  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
847 
848  // Answer.
849  return sigma;
850 
851 }
852 
853 //--------------------------------------------------------------------------
854 
855 // Select identity, colour and anticolour.
856 
857 void Sigma1ffbar2W::setIdColAcol() {
858 
859  // Sign of outgoing W.
860  int sign = 1 - 2 * (abs(id1)%2);
861  if (id1 < 0) sign = -sign;
862  setId( id1, id2, 24 * sign);
863 
864  // Colour flow topologies. Swap when antiquarks.
865  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
866  else setColAcol( 0, 0, 0, 0, 0, 0);
867  if (id1 < 0) swapColAcol();
868 
869 }
870 
871 //--------------------------------------------------------------------------
872 
873 // Evaluate weight for W decay angle.
874 
875 double Sigma1ffbar2W::weightDecay( Event& process, int iResBeg,
876  int iResEnd) {
877 
878  // W should sit in entry 5.
879  if (iResBeg != 5 || iResEnd != 5) return 1.;
880 
881  // Phase space factors.
882  double mr1 = pow2(process[6].m()) / sH;
883  double mr2 = pow2(process[7].m()) / sH;
884  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
885 
886  // Sign of asymmetry.
887  double eps = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
888 
889  // Reconstruct decay angle and weight for it.
890  double cosThe = (process[3].p() - process[4].p())
891  * (process[7].p() - process[6].p()) / (sH * betaf);
892  double wtMax = 4.;
893  double wt = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
894 
895  // Done.
896  return (wt / wtMax);
897 
898 }
899 
900 //==========================================================================
901 
902 // Sigma2ffbar2ffbarsgm class.
903 // Cross section f fbar -> gamma* -> f' fbar', for multiparton interactions.
904 
905 //--------------------------------------------------------------------------
906 
907 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
908 
909 void Sigma2ffbar2ffbarsgm::sigmaKin() {
910 
911  // Pick new flavour. Allow three leptons and five quarks.
912  double colQ = 1. + (alpS / M_PI);
913  double flavWt = 3. + colQ * 11. / 3.;
914  double flavRndm = rndmPtr->flat() * flavWt;
915  if (flavRndm < 3.) {
916  if (flavRndm < 1.) idNew = 11;
917  else if (flavRndm < 2.) idNew = 13;
918  else idNew = 15;
919  } else {
920  flavRndm = 3. * (flavRndm - 3.) / colQ;
921  if (flavRndm < 4.) idNew = 2;
922  else if (flavRndm < 8.) idNew = 4;
923  else if (flavRndm < 9.) idNew = 1;
924  else if (flavRndm < 10.) idNew = 3;
925  else idNew = 5;
926  }
927  double mNew = particleDataPtr->m0(idNew);
928  double m2New = mNew*mNew;
929 
930  // Calculate kinematics dependence. Give correct mass factors for
931  // tHat, uHat defined as if massless kinematics, d(sigma)/d(Omega)
932  // = beta (1 + cos^2(theta) + (1 - beta^2) sin^2(theta)).
933  // Special case related to phase space form in multiparton interactions.
934  double sigS = 0.;
935  if (sH > 4. * m2New) {
936  double beta = sqrt(1. - 4. * m2New / sH);
937  sigS = beta * (2.* (tH2 + uH2) + 4. * (1. - beta * beta) * tH * uH)
938  / sH2;
939  }
940 
941  // Answer is proportional to number of outgoing flavours.
942  sigma0 = (M_PI/sH2) * pow2(alpEM) * sigS * flavWt;
943 
944 }
945 
946 //--------------------------------------------------------------------------
947 
948 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
949 
950 double Sigma2ffbar2ffbarsgm::sigmaHat() {
951 
952  // Charge and colour factors.
953  double eNow = couplingsPtr->ef( abs(id1) );
954  double sigma = sigma0 * pow2(eNow);
955  if (abs(id1) < 9) sigma /= 3.;
956 
957  // Answer.
958  return sigma;
959 
960 }
961 
962 //--------------------------------------------------------------------------
963 
964 // Select identity, colour and anticolour.
965 
966 void Sigma2ffbar2ffbarsgm::setIdColAcol() {
967 
968  // Set outgoing flavours.
969  id3 = (id1 > 0) ? idNew : -idNew;
970  setId( id1, id2, id3, -id3);
971 
972  // Colour flow topologies. Swap when antiquarks.
973  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
974  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
975  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
976  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
977  if (id1 < 0) swapColAcol();
978 
979 }
980 
981 //==========================================================================
982 
983 // Sigma2ffbar2ffbarsgmZ class.
984 // Cross section f fbar -> gamma*/Z0 -> f' fbar',
985 // i.e. gamma*/Z0 decay as part of the hard process.
986 
987 //--------------------------------------------------------------------------
988 
989 // Initialize process.
990 
991 void Sigma2ffbar2ffbarsgmZ::initProc() {
992 
993  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
994  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
995 
996  // Store Z0 mass and width for propagator.
997  mRes = particleDataPtr->m0(23);
998  GammaRes = particleDataPtr->mWidth(23);
999  m2Res = mRes*mRes;
1000  GamMRat = GammaRes / mRes;
1001  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1002  * couplingsPtr->cos2thetaW());
1003 
1004  // Set pointer to particle properties and decay table.
1005  particlePtr = particleDataPtr->particleDataEntryPtr(23);
1006 
1007 }
1008 
1009 //--------------------------------------------------------------------------
1010 
1011 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1012 
1013 void Sigma2ffbar2ffbarsgmZ::sigmaKin() {
1014 
1015  // Common coupling factor.
1016  colQ = 3. * (1. + alpS / M_PI);
1017 
1018  // Reset vectors and sums. Declare variables in loop.
1019  idVec.resize(0);
1020  gamT.resize(0);
1021  gamL.resize(0);
1022  intT.resize(0);
1023  intL.resize(0);
1024  intA.resize(0);
1025  resT.resize(0);
1026  resL.resize(0);
1027  resA.resize(0);
1028  gamSumT = 0.;
1029  gamSumL = 0.;
1030  intSumT = 0.;
1031  intSumL = 0.;
1032  intSumA = 0.;
1033  resSumT = 0.;
1034  resSumL = 0.;
1035  resSumA = 0.;
1036  int onMode, idAbs;
1037  double mf, mr, betaf, ef, vf, af, colf, gamTf, gamLf, intTf, intLf,
1038  intAf, resTf, resLf, resAf;
1039 
1040  // Loop over all Z0 decay channels.
1041  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1042  onMode = particlePtr->channel(i).onMode();
1043  idAbs = abs( particlePtr->channel(i).product(0) );
1044 
1045  // Only contributions from three fermion generations, except top.
1046  if ( (onMode == 1 || onMode == 2) && ((idAbs > 0 && idAbs < 6)
1047  || ( idAbs > 10 && idAbs < 17)) ) {
1048  mf = particleDataPtr->m0(idAbs);
1049 
1050  // Check that channel above threshold. Phase space.
1051  if (mH > 2. * mf + MASSMARGIN) {
1052  mr = pow2(mf / mH);
1053  betaf = sqrtpos(1. - 4. * mr);
1054 
1055  // Combine couplings (including colour) with phase space.
1056  ef = couplingsPtr->ef(idAbs);
1057  vf = couplingsPtr->vf(idAbs);
1058  af = couplingsPtr->af(idAbs);
1059  colf = (idAbs < 6) ? colQ : 1.;
1060  gamTf = colf * ef * ef * betaf;
1061  gamLf = gamTf * 4. * mr;
1062  intTf = colf * ef * vf * betaf;
1063  intLf = intTf * 4. * mr;
1064  intAf = colf * ef * af * betaf;
1065  resTf = colf * (vf * vf * betaf + af * af * pow3(betaf));
1066  resLf = colf * vf * vf * betaf * 4. * mr;
1067  resAf = colf * vf * af * betaf * 4.;
1068 
1069  // Store individual coplings and their sums.
1070  idVec.push_back(idAbs);
1071  gamT.push_back(gamTf);
1072  gamL.push_back(gamLf);
1073  intT.push_back(intTf);
1074  intL.push_back(intLf);
1075  intA.push_back(intAf);
1076  resT.push_back(resTf);
1077  resL.push_back(resLf);
1078  resA.push_back(resAf);
1079  gamSumT += gamTf;
1080  gamSumL += gamLf;
1081  intSumT += intTf;
1082  intSumL += intLf;
1083  intSumA += intAf;
1084  resSumT += resTf;
1085  resSumL += resLf;
1086  resSumA += resAf;
1087 
1088  // End loop over Z0 decay channels.
1089  }
1090  }
1091  }
1092 
1093  // Calculate prefactors for gamma/interference/Z0 cross section terms.
1094  gamProp = M_PI * pow2(alpEM) / sH2;
1095  intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1096  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1097  resProp = gamProp * pow2(thetaWRat * sH)
1098  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1099 
1100  // Optionally only keep gamma* or Z0 term.
1101  if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1102  if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1103 
1104  // Scattering angle in subsystem rest frame.
1105  cThe = (tH - uH) / sH;
1106 
1107 }
1108 
1109 //--------------------------------------------------------------------------
1110 
1111 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1112 
1113 double Sigma2ffbar2ffbarsgmZ::sigmaHat() {
1114 
1115  // Couplings for current in-flavour.
1116  int id1Abs = abs(id1);
1117  double ei = couplingsPtr->ef(id1Abs);
1118  double vi = couplingsPtr->vf(id1Abs);
1119  double ai = couplingsPtr->af(id1Abs);
1120 
1121  // Coefficients of angular expression.
1122  double coefT = ei*ei * gamProp * gamSumT + ei*vi * intProp * intSumT
1123  + (vi*vi + ai*ai) * resProp * resSumT;
1124  double coefL = ei*ei * gamProp * gamSumL + ei*vi * intProp * intSumL
1125  + (vi*vi + ai*ai) * resProp * resSumL;
1126  double coefA = ei*ai * intProp * intSumA + vi*ai * resProp * resSumA;
1127 
1128  // Colour factor. Answer.
1129  double sigma = coefT * (1. + pow2(cThe)) + coefL * (1. - pow2(cThe))
1130  + 2. * coefA * cThe;
1131  if (id1Abs < 9) sigma /= 3.;
1132  return sigma;
1133 
1134 }
1135 
1136 //--------------------------------------------------------------------------
1137 
1138 // Select identity, colour and anticolour.
1139 
1140 void Sigma2ffbar2ffbarsgmZ::setIdColAcol() {
1141 
1142  // Couplings for chosen in-flavour.
1143  int id1Abs = abs(id1);
1144  double ei = couplingsPtr->ef(id1Abs);
1145  double vi = couplingsPtr->vf(id1Abs);
1146  double ai = couplingsPtr->af(id1Abs);
1147 
1148  // Contribution from each allowed out-flavour.
1149  sigTLA.resize(0);
1150  for (int i = 0; i < int(idVec.size()); ++i) {
1151  double coefT = ei*ei * gamProp * gamT[i] + ei*vi * intProp * intT[i]
1152  + (vi*vi + ai*ai) * resProp * resT[i];
1153  double coefL = ei*ei * gamProp * gamL[i] + ei*vi * intProp * intL[i]
1154  + (vi*vi + ai*ai) * resProp * resL[i];
1155  double coefA = ei*ai * intProp * intA[i] + vi*ai * resProp * resA[i];
1156  double sigma = coefT * (1. + pow2(cThe)) + coefL * (1. - pow2(cThe))
1157  + 2. * coefA * cThe;
1158  sigTLA.push_back(sigma);
1159  }
1160 
1161  // Pick outgoing flavours.
1162  int idNew = idVec[ rndmPtr->pick(sigTLA) ];
1163  id3 = (id1 > 0) ? idNew : -idNew;
1164  setId( id1, id2, id3, -id3);
1165 
1166  // Colour flow topologies. Swap when antiquarks.
1167  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1168  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1169  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1170  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1171  if (id1 < 0) swapColAcol();
1172 
1173 }
1174 
1175 //==========================================================================
1176 
1177 // Sigma2ffbar2ffbarsW class.
1178 // Cross section f_1 fbar_2 -> W+- -> f_3 fbar_4,
1179 // i.e. W decay as part of the hard process.
1180 
1181 //--------------------------------------------------------------------------
1182 
1183 // Initialize process.
1184 
1185 void Sigma2ffbar2ffbarsW::initProc() {
1186 
1187  // Store W+- mass and width for propagator.
1188  mRes = particleDataPtr->m0(24);
1189  GammaRes = particleDataPtr->mWidth(24);
1190  m2Res = mRes*mRes;
1191  GamMRat = GammaRes / mRes;
1192  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1193 
1194  // Set pointer to particle properties and decay table.
1195  particlePtr = particleDataPtr->particleDataEntryPtr(24);
1196 
1197 }
1198 
1199 //--------------------------------------------------------------------------
1200 
1201 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1202 
1203 void Sigma2ffbar2ffbarsW::sigmaKin() {
1204 
1205  // Set up Breit-Wigner. Common sigmaHat cross section for W+ and W-.
1206  double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1207  double preFac = alpEM * thetaWRat * mH;
1208  sigma0 = preFac * sigBW * particlePtr->resWidthOpen(24, mH);
1209 
1210  // Convert to d(sigmaHat)/d(tHat). Fixed so integrates to sigmaHat.
1211  sigma0 *= 3. * uH2 / (sH2 * sH);
1212 
1213  // Pick a decay channel.
1214  if (!particlePtr->preparePick(24, mH)) {
1215  sigma0 = 0.;
1216  return;
1217  }
1218  DecayChannel& channel = particlePtr->pickChannel();
1219  id3New = channel.product(0);
1220  id4New = channel.product(1);
1221 
1222 }
1223 
1224 //--------------------------------------------------------------------------
1225 
1226 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1227 
1228 double Sigma2ffbar2ffbarsW::sigmaHat() {
1229 
1230  // Secondary width for W+-. CKM and colour factors.
1231  double sigma = sigma0;
1232  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1233 
1234  // Answer.
1235  return sigma;
1236 
1237 }
1238 
1239 //--------------------------------------------------------------------------
1240 
1241 // Select identity, colour and anticolour.
1242 
1243 void Sigma2ffbar2ffbarsW::setIdColAcol() {
1244 
1245  // Find sign of W and its decay products. Set outgoing flavours.
1246  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
1247  id3 = (idUp > 0) ? id3New : -id3New;
1248  id4 = (idUp > 0) ? id4New : -id4New;
1249  if (id1 * id3 < 0) swap(id3, id4);
1250  setId( id1, id2, id3, id4);
1251 
1252  // Colour flow topologies. Swap when antifermions in 1 and 3.
1253  if (abs(id1) < 9 && abs(id3) < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1254  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1255  else if (abs(id3) < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1256  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1257  if (id1 < 0) swapColAcol();
1258 
1259 }
1260 
1261 //==========================================================================
1262 
1263 // Sigma2ffbar2FFbarsgmZ class.
1264 // Cross section f fbar -> gamma*/Z0 -> F Fbar.
1265 
1266 //--------------------------------------------------------------------------
1267 
1268 // Initialize process.
1269 
1270 void Sigma2ffbar2FFbarsgmZ::initProc() {
1271 
1272  // Process name.
1273  nameSave = "f fbar -> F Fbar (s-channel gamma*/Z0)";
1274  if (idNew == 4) nameSave = "f fbar -> c cbar (s-channel gamma*/Z0)";
1275  if (idNew == 5) nameSave = "f fbar -> b bbar (s-channel gamma*/Z0)";
1276  if (idNew == 6) nameSave = "f fbar -> t tbar (s-channel gamma*/Z0)";
1277  if (idNew == 7) nameSave = "f fbar -> b' b'bar (s-channel gamma*/Z0)";
1278  if (idNew == 8) nameSave = "f fbar -> t' t'bar (s-channel gamma*/Z0)";
1279  if (idNew == 15) nameSave = "f fbar -> tau+ tau- (s-channel gamma*/Z0)";
1280  if (idNew == 17) nameSave = "f fbar -> tau'+ tau'- (s-channel gamma*/Z0)";
1281  if (idNew == 18)
1282  nameSave = "f fbar -> nu'_tau nu'bar_tau (s-channel gamma*/Z0)";
1283 
1284  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1285  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
1286 
1287  // Store Z0 mass and width for propagator.
1288  mRes = particleDataPtr->m0(23);
1289  GammaRes = particleDataPtr->mWidth(23);
1290  m2Res = mRes*mRes;
1291  GamMRat = GammaRes / mRes;
1292  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1293  * couplingsPtr->cos2thetaW());
1294 
1295  // Store couplings of F.
1296  ef = couplingsPtr->ef(idNew);
1297  vf = couplingsPtr->vf(idNew);
1298  af = couplingsPtr->af(idNew);
1299 
1300  // Secondary open width fraction, relevant for top (or heavier).
1301  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
1302 
1303 }
1304 
1305 //--------------------------------------------------------------------------
1306 
1307 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1308 
1309 void Sigma2ffbar2FFbarsgmZ::sigmaKin() {
1310 
1311  // Check that above threshold.
1312  isPhysical = true;
1313  if (mH < m3 + m4 + MASSMARGIN) {
1314  isPhysical = false;
1315  return;
1316  }
1317 
1318  // Define average F, Fbar mass so same beta. Phase space.
1319  double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1320  mr = s34Avg / sH;
1321  betaf = sqrtpos(1. - 4. * mr);
1322 
1323  // Final-state colour factor.
1324  double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) : 1.;
1325 
1326  // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1327  cosThe = (tH - uH) / (betaf * sH);
1328 
1329  // Calculate prefactors for gamma/interference/Z0 cross section terms.
1330  gamProp = colF * M_PI * pow2(alpEM) / sH2;
1331  intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1332  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1333  resProp = gamProp * pow2(thetaWRat * sH)
1334  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1335 
1336  // Optionally only keep gamma* or Z0 term.
1337  if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1338  if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1339 
1340 }
1341 
1342 //--------------------------------------------------------------------------
1343 
1344 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1345 
1346 double Sigma2ffbar2FFbarsgmZ::sigmaHat() {
1347 
1348  // Fail if below threshold.
1349  if (!isPhysical) return 0.;
1350 
1351  // Couplings for in-flavours.
1352  int idAbs = abs(id1);
1353  double ei = couplingsPtr->ef(idAbs);
1354  double vi = couplingsPtr->vf(idAbs);
1355  double ai = couplingsPtr->af(idAbs);
1356 
1357  // Coefficients of angular expression.
1358  double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
1359  + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
1360  double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
1361  + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
1362  double coefAsym = betaf * ( ei * ai * intProp * ef * af
1363  + 4. * vi * ai * resProp * vf * af );
1364 
1365  // Combine gamma, interference and Z0 parts.
1366  double sigma = coefTran * (1. + pow2(cosThe))
1367  + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
1368 
1369  // Top: corrections for closed decay channels.
1370  sigma *= openFracPair;
1371 
1372  // Initial-state colour factor. Answer.
1373  if (idAbs < 9) sigma /= 3.;
1374  return sigma;
1375 
1376 }
1377 
1378 //--------------------------------------------------------------------------
1379 
1380 // Select identity, colour and anticolour.
1381 
1382 void Sigma2ffbar2FFbarsgmZ::setIdColAcol() {
1383 
1384  // Set outgoing flavours.
1385  id3 = (id1 > 0) ? idNew : -idNew;
1386  setId( id1, id2, id3, -id3);
1387 
1388  // Colour flow topologies. Swap when antiquarks.
1389  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1390  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1391  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1392  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1393  if (id1 < 0) swapColAcol();
1394 
1395 }
1396 
1397 //--------------------------------------------------------------------------
1398 
1399 // Evaluate weight for decay angles of W in top decay.
1400 
1401 double Sigma2ffbar2FFbarsgmZ::weightDecay( Event& process, int iResBeg,
1402  int iResEnd) {
1403 
1404  // For top decay hand over to standard routine, else done.
1405  if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1406  return weightTopDecay( process, iResBeg, iResEnd);
1407  else return 1.;
1408 
1409 }
1410 
1411 //==========================================================================
1412 
1413 // Sigma2ffbar2FfbarsW class.
1414 // Cross section f fbar' -> W+- -> F fbar".
1415 
1416 //--------------------------------------------------------------------------
1417 
1418 // Initialize process.
1419 
1420 void Sigma2ffbar2FfbarsW::initProc() {
1421 
1422  // Process name.
1423  nameSave = "f fbar -> F fbar (s-channel W+-)";
1424  if (idNew == 4) nameSave = "f fbar -> c qbar (s-channel W+-)";
1425  if (idNew == 5) nameSave = "f fbar -> b qbar (s-channel W+-)";
1426  if (idNew == 6) nameSave = "f fbar -> t qbar (s-channel W+-)";
1427  if (idNew == 7) nameSave = "f fbar -> b' qbar (s-channel W+-)";
1428  if (idNew == 8) nameSave = "f fbar -> t' qbar (s-channel W+-)";
1429  if (idNew == 7 && idNew2 == 6)
1430  nameSave = "f fbar -> b' tbar (s-channel W+-)";
1431  if (idNew == 8 && idNew2 == 7)
1432  nameSave = "f fbar -> t' b'bar (s-channel W+-)";
1433  if (idNew == 15 || idNew == 16)
1434  nameSave = "f fbar -> tau nu_taubar (s-channel W+-)";
1435  if (idNew == 17 || idNew == 18)
1436  nameSave = "f fbar -> tau' nu'_taubar (s-channel W+-)";
1437 
1438  // Store W+- mass and width for propagator.
1439  mRes = particleDataPtr->m0(24);
1440  GammaRes = particleDataPtr->mWidth(24);
1441  m2Res = mRes*mRes;
1442  GamMRat = GammaRes / mRes;
1443  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1444 
1445  // For t/t' want to use at least b mass.
1446  idPartner = idNew2;
1447  if ( (idNew == 6 || idNew == 8) && idNew2 == 0 ) idPartner = 5;
1448 
1449  // Sum of CKM weights for quarks.
1450  V2New = (idNew < 9) ? couplingsPtr->V2CKMsum(idNew) : 1.;
1451  if (idNew2 != 0) V2New = couplingsPtr->V2CKMid(idNew, idNew2);
1452 
1453  // Secondary open width fractions, relevant for top or heavier.
1454  openFracPos = particleDataPtr->resOpenFrac( idNew, -idNew2);
1455  openFracNeg = particleDataPtr->resOpenFrac(-idNew, idNew2);
1456 
1457 }
1458 
1459 //--------------------------------------------------------------------------
1460 
1461 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1462 
1463 void Sigma2ffbar2FfbarsW::sigmaKin() {
1464 
1465  // Check that above threshold.
1466  isPhysical = true;
1467  if (mH < m3 + m4 + MASSMARGIN) {
1468  isPhysical = false;
1469  return;
1470  }
1471 
1472  // Phase space factors.
1473  double mr1 = s3 / sH;
1474  double mr2 = s4 / sH;
1475  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
1476 
1477  // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1478  double cosThe = (tH - uH) / (betaf * sH);
1479 
1480  // Set up Breit-Wigner and in- and out-widths.
1481  double sigBW = 9. * M_PI * pow2(alpEM * thetaWRat)
1482  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1483 
1484  // Initial-state colour factor.
1485  double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) * V2New : 1.;
1486 
1487  // Angular dependence.
1488  double wt = pow2(1. + betaf * cosThe) - pow2(mr1 - mr2);
1489 
1490  // Temporary answer.
1491  sigma0 = sigBW * colF * wt;
1492 
1493 }
1494 
1495 //--------------------------------------------------------------------------
1496 
1497 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1498 
1499 double Sigma2ffbar2FfbarsW::sigmaHat() {
1500 
1501  // Fail if below threshold.
1502  if (!isPhysical) return 0.;
1503 
1504  // CKM and colour factors.
1505  double sigma = sigma0;
1506  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1507 
1508  // Correction for secondary width in top (or heavier) decay.
1509  int idSame = ((abs(id1) + idNew)%2 == 0) ? id1 : id2;
1510  sigma *= (idSame > 0) ? openFracPos : openFracNeg;
1511 
1512  // Answer.
1513  return sigma;
1514 
1515 }
1516 
1517 //--------------------------------------------------------------------------
1518 
1519 // Select identity, colour and anticolour.
1520 
1521 void Sigma2ffbar2FfbarsW::setIdColAcol() {
1522 
1523  // Set outgoing flavours.
1524  id3 = idNew;
1525  id4 = (idNew2 != 0) ? idNew2 : couplingsPtr->V2CKMpick(idNew);
1526  if (idNew%2 == 0) {
1527  int idInUp = (abs(id1)%2 == 0) ? id1 : id2;
1528  if (idInUp > 0) id4 = -id4;
1529  else id3 = -id3;
1530  } else {
1531  int idInDn = (abs(id1)%2 == 1) ? id1 : id2;
1532  if (idInDn > 0) id4 = -id4;
1533  else id3 = -id3;
1534  }
1535  setId( id1, id2, id3, id4);
1536 
1537  // Swap tHat and uHat for fbar' f -> F f".
1538  if (id1 * id3 < 0) swapTU = true;
1539 
1540  // Colour flow topologies. Swap when antiquarks.
1541  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1542  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1543  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1544  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1545  if (id1 < 0) swapCol12();
1546  if (id3 < 0) swapCol34();
1547 
1548 }
1549 
1550 //--------------------------------------------------------------------------
1551 
1552 // Evaluate weight for decay angles of W in top decay.
1553 
1554 double Sigma2ffbar2FfbarsW::weightDecay( Event& process, int iResBeg,
1555  int iResEnd) {
1556 
1557  // For top decay hand over to standard routine, else done.
1558  if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1559  return weightTopDecay( process, iResBeg, iResEnd);
1560  else return 1.;
1561 
1562 }
1563 
1564 //==========================================================================
1565 
1566 // Sigma2ffbargmZWgmZW class.
1567 // Collects common methods for f fbar -> gamma*/Z0/W+- gamma*/Z0/W-+.
1568 
1569 //--------------------------------------------------------------------------
1570 
1571 // Calculate and store internal products.
1572 
1573 void Sigma2ffbargmZWgmZW::setupProd( Event& process, int i1, int i2,
1574  int i3, int i4, int i5, int i6) {
1575 
1576  // Store incoming and outgoing momenta,
1577  pRot[1] = process[i1].p();
1578  pRot[2] = process[i2].p();
1579  pRot[3] = process[i3].p();
1580  pRot[4] = process[i4].p();
1581  pRot[5] = process[i5].p();
1582  pRot[6] = process[i6].p();
1583 
1584  // Do random rotation to avoid accidental zeroes in HA expressions.
1585  bool smallPT = false;
1586  do {
1587  smallPT = false;
1588  double thetaNow = acos(2. * rndmPtr->flat() - 1.);
1589  double phiNow = 2. * M_PI * rndmPtr->flat();
1590  for (int i = 1; i <= 6; ++i) {
1591  pRot[i].rot( thetaNow, phiNow);
1592  if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
1593  }
1594  } while (smallPT);
1595 
1596  // Calculate internal products.
1597  for (int i = 1; i < 6; ++i) {
1598  for (int j = i + 1; j <= 6; ++j) {
1599  hA[i][j] =
1600  sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
1601  / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
1602  - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
1603  / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
1604  hC[i][j] = conj( hA[i][j] );
1605  if (i <= 2) {
1606  hA[i][j] *= complex( 0., 1.);
1607  hC[i][j] *= complex( 0., 1.);
1608  }
1609  hA[j][i] = - hA[i][j];
1610  hC[j][i] = - hC[i][j];
1611  }
1612  }
1613 
1614 }
1615 
1616 //--------------------------------------------------------------------------
1617 
1618 // Evaluate the F function of Gunion and Kunszt.
1619 
1620 complex Sigma2ffbargmZWgmZW::fGK(int j1, int j2, int j3, int j4, int j5,
1621  int j6) {
1622 
1623  return 4. * hA[j1][j3] * hC[j2][j6]
1624  * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
1625 
1626 }
1627 
1628 //--------------------------------------------------------------------------
1629 
1630 // Evaluate the Xi function of Gunion and Kunszt.
1631 
1632 double Sigma2ffbargmZWgmZW::xiGK( double tHnow, double uHnow) {
1633 
1634  return - 4. * s3 * s4 + tHnow * (3. * tHnow + 4. * uHnow)
1635  + tHnow * tHnow * ( tHnow * uHnow / (s3 * s4)
1636  - 2. * (1. / s3 + 1./s4) * (tHnow + uHnow)
1637  + 2. * (s3 / s4 + s4 / s3) );
1638 
1639 }
1640 
1641 //--------------------------------------------------------------------------
1642 
1643 // Evaluate the Xj function of Gunion and Kunszt.
1644 
1645 double Sigma2ffbargmZWgmZW::xjGK( double tHnow, double uHnow) {
1646 
1647  return 8. * pow2(s3 + s4) - 8. * (s3 + s4) * (tHnow + uHnow)
1648  - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
1649  / (s3 * s4) - 2. * (1. / s3 + 1. / s4) * (tHnow + uHnow)
1650  + 2. * (s3 / s4 + s4 / s3) );
1651 
1652 }
1653 
1654 //==========================================================================
1655 
1656 // Sigma2ffbar2gmZgmZ class.
1657 // Cross section for f fbar -> gamma*/Z0 gamma*/Z0 (f is quark or lepton).
1658 
1659 //--------------------------------------------------------------------------
1660 
1661 // Initialize process.
1662 
1663 void Sigma2ffbar2gmZgmZ::initProc() {
1664 
1665  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1666  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
1667 
1668  // Store Z0 mass and width for propagator.
1669  mRes = particleDataPtr->m0(23);
1670  GammaRes = particleDataPtr->mWidth(23);
1671  m2Res = mRes*mRes;
1672  GamMRat = GammaRes / mRes;
1673  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1674  * couplingsPtr->cos2thetaW());
1675 
1676  // Set pointer to particle properties and decay table.
1677  particlePtr = particleDataPtr->particleDataEntryPtr(23);
1678 
1679 }
1680 
1681 //--------------------------------------------------------------------------
1682 
1683 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1684 
1685 void Sigma2ffbar2gmZgmZ::sigmaKin() {
1686 
1687  // Cross section part common for all incoming flavours.
1688  sigma0 = (M_PI / sH2) * pow2(alpEM) * 0.5
1689  * ( (tH2 + uH2 + 2. * (s3 + s4) * sH) / (tH * uH)
1690  - s3 * s4 * (1./tH2 + 1./uH2) );
1691 
1692  // Common coupling factors at the resonance masses
1693  double alpEM3 = couplingsPtr->alphaEM(s3);
1694  double alpS3 = couplingsPtr->alphaS(s3);
1695  double colQ3 = 3. * (1. + alpS3 / M_PI);
1696  double alpEM4 = couplingsPtr->alphaEM(s4);
1697  double alpS4 = couplingsPtr->alphaS(s4);
1698  double colQ4 = 3. * (1. + alpS4 / M_PI);
1699 
1700  // Reset quantities to sum. Declare variables in loop.
1701  gamSum3 = 0.;
1702  intSum3 = 0.;
1703  resSum3 = 0.;
1704  gamSum4 = 0.;
1705  intSum4 = 0.;
1706  resSum4 = 0.;
1707  int onMode;
1708  double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
1709 
1710  // Loop over all Z0 decay channels.
1711  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1712  int idAbs = abs( particlePtr->channel(i).product(0) );
1713 
1714  // Only contributions from three fermion generations, except top.
1715  if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
1716  mf = particleDataPtr->m0(idAbs);
1717  onMode = particlePtr->channel(i).onMode();
1718 
1719  // First Z0: check that above threshold. Phase space.
1720  if (m3 > 2. * mf + MASSMARGIN) {
1721  mr = pow2(mf / m3);
1722  betaf = sqrtpos(1. - 4. * mr);
1723  psvec = betaf * (1. + 2. * mr);
1724  psaxi = pow3(betaf);
1725 
1726  // First Z0: combine phase space with couplings.
1727  ef2 = couplingsPtr->ef2(idAbs) * psvec;
1728  efvf = couplingsPtr->efvf(idAbs) * psvec;
1729  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1730  + couplingsPtr->af2(idAbs) * psaxi;
1731  colf = (idAbs < 6) ? colQ3 : 1.;
1732 
1733  // First Z0: store sum of combinations for open outstate channels.
1734  if (onMode == 1 || onMode == 2) {
1735  gamSum3 += colf * ef2;
1736  intSum3 += colf * efvf;
1737  resSum3 += colf * vf2af2;
1738  }
1739  }
1740 
1741  // Second Z0: check that above threshold. Phase space.
1742  if (m4 > 2. * mf + MASSMARGIN) {
1743  mr = pow2(mf / m4);
1744  betaf = sqrtpos(1. - 4. * mr);
1745  psvec = betaf * (1. + 2. * mr);
1746  psaxi = pow3(betaf);
1747 
1748  // Second Z0: combine phase space with couplings.
1749  ef2 = couplingsPtr->ef2(idAbs) * psvec;
1750  efvf = couplingsPtr->efvf(idAbs) * psvec;
1751  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1752  + couplingsPtr->af2(idAbs) * psaxi;
1753  colf = (idAbs < 6) ? colQ4 : 1.;
1754 
1755  // Second Z0: store sum of combinations for open outstate channels.
1756  if (onMode == 1 || onMode == 2) {
1757  gamSum4 += colf * ef2;
1758  intSum4 += colf * efvf;
1759  resSum4 += colf * vf2af2;
1760  }
1761  }
1762 
1763  // End loop over fermions.
1764  }
1765  }
1766 
1767  // First Z0: calculate prefactors for gamma/interference/Z0 terms.
1768  gamProp3 = 4. * alpEM3 / (3. * M_PI * s3);
1769  intProp3 = gamProp3 * 2. * thetaWRat * s3 * (s3 - m2Res)
1770  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1771  resProp3 = gamProp3 * pow2(thetaWRat * s3)
1772  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1773 
1774  // First Z0: optionally only keep gamma* or Z0 term.
1775  if (gmZmode == 1) {intProp3 = 0.; resProp3 = 0.;}
1776  if (gmZmode == 2) {gamProp3 = 0.; intProp3 = 0.;}
1777 
1778  // Second Z0: calculate prefactors for gamma/interference/Z0 terms.
1779  gamProp4 = 4. * alpEM4 / (3. * M_PI * s4);
1780  intProp4 = gamProp4 * 2. * thetaWRat * s4 * (s4 - m2Res)
1781  / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1782  resProp4 = gamProp4 * pow2(thetaWRat * s4)
1783  / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1784 
1785  // Second Z0: optionally only keep gamma* or Z0 term.
1786  if (gmZmode == 1) {intProp4 = 0.; resProp4 = 0.;}
1787  if (gmZmode == 2) {gamProp4 = 0.; intProp4 = 0.;}
1788 
1789 }
1790 
1791 //--------------------------------------------------------------------------
1792 
1793 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1794 
1795 double Sigma2ffbar2gmZgmZ::sigmaHat() {
1796 
1797  // Charge/2, left- and righthanded couplings for in-fermion.
1798  int idAbs = abs(id1);
1799  double ei = 0.5 * couplingsPtr->ef(idAbs);
1800  double li = couplingsPtr->lf(idAbs);
1801  double ri = couplingsPtr->rf(idAbs);
1802 
1803  // Combine left/right gamma, interference and Z0 parts for each Z0.
1804  double left3 = ei * ei * gamProp3 * gamSum3
1805  + ei * li * intProp3 * intSum3
1806  + li * li * resProp3 * resSum3;
1807  double right3 = ei * ei * gamProp3 * gamSum3
1808  + ei * ri * intProp3 * intSum3
1809  + ri * ri * resProp3 * resSum3;
1810  double left4 = ei * ei * gamProp4 * gamSum4
1811  + ei * li * intProp4 * intSum4
1812  + li * li * resProp4 * resSum4;
1813  double right4 = ei * ei * gamProp4 * gamSum4
1814  + ei * ri * intProp4 * intSum4
1815  + ri * ri * resProp4 * resSum4;
1816 
1817  // Combine left- and right-handed couplings for the two Z0's.
1818  double sigma = sigma0 * (left3 * left4 + right3 * right4);
1819 
1820  // Correct for the running-width Z0 propagators weight in PhaseSpace.
1821  sigma /= (runBW3 * runBW4);
1822 
1823  // Initial-state colour factor. Answer.
1824  if (idAbs < 9) sigma /= 3.;
1825  return sigma;
1826 
1827 }
1828 
1829 //--------------------------------------------------------------------------
1830 
1831 // Select identity, colour and anticolour.
1832 
1833 void Sigma2ffbar2gmZgmZ::setIdColAcol() {
1834 
1835  // Flavours trivial.
1836  setId( id1, id2, 23, 23);
1837 
1838  // Colour flow topologies. Swap when antiquarks.
1839  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1840  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1841  if (id1 < 0) swapColAcol();
1842 
1843 }
1844 
1845 //--------------------------------------------------------------------------
1846 
1847 // Evaluate correlated decay flavours of the two gamma*/Z0.
1848 // Unique complication, caused by gamma*/Z0 mix different left/right.
1849 
1850 double Sigma2ffbar2gmZgmZ::weightDecayFlav( Event& process) {
1851 
1852  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1853  i1 = (process[3].id() < 0) ? 3 : 4;
1854  i2 = 7 - i1;
1855  i3 = (process[7].id() > 0) ? 7 : 8;
1856  i4 = 15 - i3;
1857  i5 = (process[9].id() > 0) ? 9 : 10;
1858  i6 = 19 - i5;
1859 
1860  // Charge/2, left- and righthanded couplings for in- and out-fermions.
1861  int idAbs = process[i1].idAbs();
1862  double ei = 0.5 * couplingsPtr->ef(idAbs);
1863  double li = couplingsPtr->lf(idAbs);
1864  double ri = couplingsPtr->rf(idAbs);
1865  idAbs = process[i3].idAbs();
1866  double e3 = 0.5 * couplingsPtr->ef(idAbs);
1867  double l3 = couplingsPtr->lf(idAbs);
1868  double r3 = couplingsPtr->rf(idAbs);
1869  idAbs = process[i5].idAbs();
1870  double e4 = 0.5 * couplingsPtr->ef(idAbs);
1871  double l4 = couplingsPtr->lf(idAbs);
1872  double r4 = couplingsPtr->rf(idAbs);
1873 
1874  // Left- and righthanded couplings combined with propagators.
1875  c3LL = ei * ei * gamProp3 * e3 * e3
1876  + ei * li * intProp3 * e3 * l3
1877  + li * li * resProp3 * l3 * l3;
1878  c3LR = ei * ei * gamProp3 * e3 * e3
1879  + ei * li * intProp3 * e3 * r3
1880  + li * li * resProp3 * r3 * r3;
1881  c3RL = ei * ei * gamProp3 * e3 * e3
1882  + ei * ri * intProp3 * e3 * l3
1883  + ri * ri * resProp3 * l3 * l3;
1884  c3RR = ei * ei * gamProp3 * e3 * e3
1885  + ei * ri * intProp3 * e3 * r3
1886  + ri * ri * resProp3 * r3 * r3;
1887  c4LL = ei * ei * gamProp4 * e4 * e4
1888  + ei * li * intProp4 * e4 * l4
1889  + li * li * resProp4 * l4 * l4;
1890  c4LR = ei * ei * gamProp4 * e4 * e4
1891  + ei * li * intProp4 * e4 * r4
1892  + li * li * resProp4 * r4 * r4;
1893  c4RL = ei * ei * gamProp4 * e4 * e4
1894  + ei * ri * intProp4 * e4 * l4
1895  + ri * ri * resProp4 * l4 * l4;
1896  c4RR = ei * ei * gamProp4 * e4 * e4
1897  + ei * ri * intProp4 * e4 * r4
1898  + ri * ri * resProp4 * r4 * r4;
1899 
1900  // Flavour weight and maximum.
1901  flavWt = (c3LL + c3LR) * (c4LL + c4LR) + (c3RL + c3RR) * (c4RL + c4RR);
1902  double flavWtMax = (c3LL + c3LR + c3RL + c3RR) * (c4LL + c4LR + c4RL + c4RR);
1903 
1904  // Done.
1905  return flavWt / flavWtMax;
1906 
1907 }
1908 
1909 //--------------------------------------------------------------------------
1910 
1911 // Evaluate weight for decay angles of the two gamma*/Z0.
1912 
1913 double Sigma2ffbar2gmZgmZ::weightDecay( Event& process, int iResBeg,
1914  int iResEnd) {
1915 
1916  // Two resonance decays, but with common weight.
1917  if (iResBeg != 5 || iResEnd != 6) return 1.;
1918 
1919  // Set up four-products and internal products.
1920  setupProd( process, i1, i2, i3, i4, i5, i6);
1921 
1922  // Flip tHat and uHat if first incoming is fermion.
1923  double tHres = tH;
1924  double uHres = uH;
1925  if (process[3].id() > 0) swap( tHres, uHres);
1926 
1927  // Kinematics factors (norm(x) = |x|^2).
1928  double fGK135 = norm( fGK( 1, 2, 3, 4, 5, 6) / tHres
1929  + fGK( 1, 2, 5, 6, 3, 4) / uHres );
1930  double fGK145 = norm( fGK( 1, 2, 4, 3, 5, 6) / tHres
1931  + fGK( 1, 2, 5, 6, 4, 3) / uHres );
1932  double fGK136 = norm( fGK( 1, 2, 3, 4, 6, 5) / tHres
1933  + fGK( 1, 2, 6, 5, 3, 4) / uHres );
1934  double fGK146 = norm( fGK( 1, 2, 4, 3, 6, 5) / tHres
1935  + fGK( 1, 2, 6, 5, 4, 3) / uHres );
1936  double fGK253 = norm( fGK( 2, 1, 5, 6, 3, 4) / tHres
1937  + fGK( 2, 1, 3, 4, 5, 6) / uHres );
1938  double fGK263 = norm( fGK( 2, 1, 6, 5, 3, 4) / tHres
1939  + fGK( 2, 1, 3, 4, 6, 5) / uHres );
1940  double fGK254 = norm( fGK( 2, 1, 5, 6, 4, 3) / tHres
1941  + fGK( 2, 1, 4, 3, 5, 6) / uHres );
1942  double fGK264 = norm( fGK( 2, 1, 6, 5, 4, 3) / tHres
1943  + fGK( 2, 1, 4, 3, 6, 5) / uHres );
1944 
1945  // Weight and maximum.
1946  double wt = c3LL * c4LL * fGK135 + c3LR * c4LL * fGK145
1947  + c3LL * c4LR * fGK136 + c3LR * c4LR * fGK146
1948  + c3RL * c4RL * fGK253 + c3RR * c4RL * fGK263
1949  + c3RL * c4RR * fGK254 + c3RR * c4RR * fGK264;
1950  double wtMax = 16. * s3 * s4 * flavWt
1951  * ( (tHres*tHres + uHres*uHres + 2. * sH * (s3 + s4)) / (tHres * uHres)
1952  - s3 * s4 * (1. / (tHres*tHres) + 1. / (uHres*uHres)) );
1953 
1954  // Done.
1955  return wt / wtMax;
1956 
1957 }
1958 
1959 //==========================================================================
1960 
1961 // Sigma2ffbar2ZW class.
1962 // Cross section for f fbar' -> Z0 W+- (f is quark or lepton).
1963 
1964 //--------------------------------------------------------------------------
1965 
1966 // Initialize process.
1967 
1968 void Sigma2ffbar2ZW::initProc() {
1969 
1970  // Store W+- mass and width for propagator.
1971  mW = particleDataPtr->m0(24);
1972  widW = particleDataPtr->mWidth(24);
1973  mWS = mW*mW;
1974  mwWS = pow2(mW * widW);
1975 
1976  // Left-handed couplings for up/nu- and down/e-type quarks.
1977  lun = (hasLeptonBeams) ? couplingsPtr->lf(12) : couplingsPtr->lf(2);
1978  lde = (hasLeptonBeams) ? couplingsPtr->lf(11) : couplingsPtr->lf(1);
1979 
1980  // Common weak coupling factor.
1981  sin2thetaW = couplingsPtr->sin2thetaW();
1982  cos2thetaW = couplingsPtr->cos2thetaW();
1983  thetaWRat = 1. / (4. * cos2thetaW);
1984  cotT = sqrt(cos2thetaW / sin2thetaW);
1985  thetaWpt = (9. - 8. * sin2thetaW) / 4.;
1986  thetaWmm = (8. * sin2thetaW - 6.) / 4.;
1987 
1988  // Secondary open width fractions.
1989  openFracPos = particleDataPtr->resOpenFrac(23, 24);
1990  openFracNeg = particleDataPtr->resOpenFrac(23, -24);
1991 
1992 }
1993 
1994 //--------------------------------------------------------------------------
1995 
1996 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1997 
1998 void Sigma2ffbar2ZW::sigmaKin() {
1999 
2000  // Evaluate cross section, as programmed by Merlin Kole (after tidying),
2001  // based on Brown, Sahdev, Mikaelian, Phys Rev. D20 (1979) 1069.
2002  /*
2003  double resBW = 1. / (pow2(sH - mWS) + mwWS);
2004  double prefac = 12.0 * M_PI * pow2(alpEM) / (sH2 * 8. * sin2thetaW);
2005  double temp1 = tH * uH - s3 * s4;
2006  double temp2 = temp1 / (s3 * s4);
2007  double temp3 = (s3 + s4) / sH;
2008  double temp4 = s3 * s4 / sH;
2009  double partA = temp2 * (0.25 - 0.5 * temp3
2010  + (pow2(s3 + s4) + 8. * s3 * s4)/(4. * sH2) )
2011  + (s3 + s4)/(s3 * s4) * (sH/2. - s3 - s4 + pow2(s3 - s4)/(2. * sH));
2012  double partB1 = lun * (0.25 * temp2 * (1. - temp3 - 4. * temp4 / tH)
2013  + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / tH) );
2014  double partB2 = lde * ( 0.25 * temp2 * (1.- temp3 - 4. * temp4 / uH)
2015  + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / uH) );
2016  double partE = 0.25 * temp2 + 0.5 *(s3 + s4) / temp4;
2017  sigma0 = prefac * (pow2(cotT) * sH2 * resBW * partA
2018  + 2.* sH * cotT * resBW * (sH - mWS) * (partB2 - partB1)
2019  + pow2(lun - lde) * partE + pow2(lde) * temp1/uH2
2020  + pow2(lun) * temp1/tH2 + 2. * lun * lde * sH * (s3 + s4) / (uH * tH));
2021  */
2022 
2023  // Evaluate cross section. Expression from EHLQ, with bug fix,
2024  // but can still give negative cross section so suspect.
2025  double resBW = 1. / (pow2(sH - mWS) + mwWS);
2026  sigma0 = (M_PI / sH2) * 0.5 * pow2(alpEM / sin2thetaW);
2027  sigma0 *= sH * resBW * (thetaWpt * pT2 + thetaWmm * (s3 + s4))
2028  + (sH - mWS) * resBW * sH * (pT2 - s3 - s4) * (lun / tH - lde / uH)
2029  + thetaWRat * sH * pT2 * ( lun*lun / tH2 + lde*lde / uH2 )
2030  + 2. * thetaWRat * sH * (s3 + s4) * lun * lde / (tH * uH);
2031  // Need to protect against negative cross sections at times.
2032  sigma0 = max(0., sigma0);
2033 
2034 }
2035 
2036 //--------------------------------------------------------------------------
2037 
2038 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2039 
2040 double Sigma2ffbar2ZW::sigmaHat() {
2041 
2042  // CKM and colour factors.
2043  double sigma = sigma0;
2044  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
2045 
2046  // Corrections for secondary widths in Z0 and W+- decays.
2047  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2048  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2049 
2050  // Answer.
2051  return sigma;
2052 
2053 }
2054 
2055 //--------------------------------------------------------------------------
2056 
2057 // Select identity, colour and anticolour.
2058 
2059 void Sigma2ffbar2ZW::setIdColAcol() {
2060 
2061  // Sign of outgoing W.
2062  int sign = 1 - 2 * (abs(id1)%2);
2063  if (id1 < 0) sign = -sign;
2064  setId( id1, id2, 23, 24 * sign);
2065 
2066  // tHat is defined between (f, W-) or (fbar, W+),
2067  // so OK for u/ubar on side 1, but must swap tHat <-> uHat if d/dbar.
2068  if (abs(id1)%2 == 1) swapTU = true;
2069 
2070  // Colour flow topologies. Swap when antiquarks.
2071  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2072  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2073  if (id1 < 0) swapColAcol();
2074 
2075 }
2076 
2077 //--------------------------------------------------------------------------
2078 
2079 // Evaluate weight for Z0 and W+- decay angles.
2080 
2081 double Sigma2ffbar2ZW::weightDecay( Event& process, int iResBeg, int iResEnd) {
2082 
2083  // Two resonance decays, but with common weight.
2084  if (iResBeg != 5 || iResEnd != 6) return 1.;
2085 
2086  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6)
2087  // with f' fbar' from W+- and f" fbar" from Z0 (note flip Z0 <-> W+-).
2088  int i1 = (process[3].id() < 0) ? 3 : 4;
2089  int i2 = 7 - i1;
2090  int i3 = (process[9].id() > 0) ? 9 : 10;
2091  int i4 = 19 - i3;
2092  int i5 = (process[7].id() > 0) ? 7 : 8;
2093  int i6 = 15 - i5;
2094 
2095  // Set up four-products and internal products.
2096  setupProd( process, i1, i2, i3, i4, i5, i6);
2097 
2098  // Swap tHat and uHat if incoming fermion is downtype.
2099  double tHres = tH;
2100  double uHres = uH;
2101  if (process[i2].id()%2 == 1) swap( tHres, uHres);
2102 
2103  // Couplings of incoming (anti)fermions and outgoing from Z0.
2104  int idAbs = process[i1].idAbs();
2105  double ai = couplingsPtr->af(idAbs);
2106  double li1 = couplingsPtr->lf(idAbs);
2107  idAbs = process[i2].idAbs();
2108  double li2 = couplingsPtr->lf(idAbs);
2109  idAbs = process[i5].idAbs();
2110  double l4 = couplingsPtr->lf(idAbs);
2111  double r4 = couplingsPtr->rf(idAbs);
2112 
2113  // W propagator/interference factor.
2114  double Wint = cos2thetaW * (sH - mWS) / (pow2(sH - mWS) + mwWS);
2115 
2116  // Combinations of couplings and kinematics (norm(x) = |x|^2).
2117  double aWZ = li2 / tHres - 2. * Wint * ai;
2118  double bWZ = li1 / uHres + 2. * Wint * ai;
2119  double fGK135 = norm( aWZ * fGK( 1, 2, 3, 4, 5, 6)
2120  + bWZ * fGK( 1, 2, 5, 6, 3, 4) );
2121  double fGK136 = norm( aWZ * fGK( 1, 2, 3, 4, 6, 5)
2122  + bWZ * fGK( 1, 2, 6, 5, 3, 4) );
2123  double xiT = xiGK( tHres, uHres);
2124  double xiU = xiGK( uHres, tHres);
2125  double xjTU = xjGK( tHres, uHres);
2126 
2127  // Weight and maximum weight.
2128  double wt = l4*l4 * fGK135 + r4*r4 * fGK136;
2129  double wtMax = 4. * s3 * s4 * (l4*l4 + r4*r4)
2130  * (aWZ * aWZ * xiT + bWZ * bWZ * xiU + aWZ * bWZ * xjTU);
2131 
2132  // Done.
2133  return wt / wtMax;
2134 
2135 }
2136 
2137 //==========================================================================
2138 
2139 // Sigma2ffbar2WW class.
2140 // Cross section for f fbar -> W- W+ (f is quark or lepton).
2141 
2142 //--------------------------------------------------------------------------
2143 
2144 // Initialize process.
2145 
2146 void Sigma2ffbar2WW::initProc() {
2147 
2148  // Store Z0 mass and width for propagator. Common coupling factor.
2149  mZ = particleDataPtr->m0(23);
2150  widZ = particleDataPtr->mWidth(23);
2151  mZS = mZ*mZ;
2152  mwZS = pow2(mZ * widZ);
2153  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
2154 
2155  // Secondary open width fraction.
2156  openFracPair = particleDataPtr->resOpenFrac(24, -24);
2157 
2158 }
2159 
2160 //--------------------------------------------------------------------------
2161 
2162 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2163 
2164 void Sigma2ffbar2WW::sigmaKin() {
2165 
2166  // Cross section part common for all incoming flavours.
2167  sigma0 = (M_PI / sH2) * pow2(alpEM);
2168 
2169  // Z0 propagator and gamma*/Z0 interference.
2170  double Zprop = sH2 / (pow2(sH - mZS) + mwZS);
2171  double Zint = Zprop * (1. - mZS / sH);
2172 
2173  // Common coupling factors (g = gamma*, Z = Z0, f = t-channel fermion).
2174  cgg = 0.5;
2175  cgZ = thetaWRat * Zint;
2176  cZZ = 0.5 * pow2(thetaWRat) * Zprop;
2177  cfg = thetaWRat;
2178  cfZ = pow2(thetaWRat) * Zint;
2179  cff = pow2(thetaWRat);
2180 
2181  // Kinematical functions.
2182  double rat34 = sH * (2. * (s3 + s4) + pT2) / (s3 * s4);
2183  double lambdaS = pow2(sH - s3 - s4) - 4. * s3 * s4;
2184  double intA = (sH - s3 - s4) * rat34 / sH;
2185  double intB = 4. * (s3 + s4 - pT2);
2186  gSS = (lambdaS * rat34 + 12. * sH * pT2) / sH2;
2187  gTT = rat34 + 4. * sH * pT2 / tH2;
2188  gST = intA + intB / tH;
2189  gUU = rat34 + 4. * sH * pT2 / uH2;
2190  gSU = intA + intB / uH;
2191 
2192 }
2193 
2194 //--------------------------------------------------------------------------
2195 
2196 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2197 
2198 double Sigma2ffbar2WW::sigmaHat() {
2199 
2200  // Flavour-specific couplings.
2201  int idAbs = abs(id1);
2202  double ei = couplingsPtr->ef(idAbs);
2203  double vi = couplingsPtr->vf(idAbs);
2204  double ai = couplingsPtr->af(idAbs);
2205 
2206  // Combine, with different cases for up- and down-type in-flavours.
2207  double sigma = sigma0;
2208  sigma *= (idAbs%2 == 1)
2209  ? (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
2210  + (cfg * ei + cfZ * (vi + ai)) * gST + cff * gTT
2211  : (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
2212  - (cfg * ei + cfZ * (vi + ai)) * gSU + cff * gUU;
2213 
2214  // Initial-state colour factor. Correction for secondary widths. Answer.
2215  if (idAbs < 9) sigma /= 3.;
2216  sigma *= openFracPair;
2217  return sigma;
2218 
2219 }
2220 
2221 //--------------------------------------------------------------------------
2222 
2223 // Select identity, colour and anticolour.
2224 
2225 void Sigma2ffbar2WW::setIdColAcol() {
2226 
2227  // Always order W- W+, i.e. W- first.
2228  setId( id1, id2, -24, 24);
2229 
2230  // tHat is defined between (f, W-) or (fbar, W+),
2231  if (id1 < 0) swapTU = true;
2232 
2233  // Colour flow topologies. Swap when antiquarks.
2234  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2235  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2236  if (id1 < 0) swapColAcol();
2237 
2238 }
2239 
2240 //--------------------------------------------------------------------------
2241 
2242 // Evaluate weight for W+ and W- decay angles.
2243 
2244 double Sigma2ffbar2WW::weightDecay( Event& process, int iResBeg, int iResEnd) {
2245 
2246  // Two resonance decays, but with common weight.
2247  if (iResBeg != 5 || iResEnd != 6) return 1.;
2248 
2249  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
2250  // with f' fbar' from W- and f" fbar" from W+.
2251  int i1 = (process[3].id() < 0) ? 3 : 4;
2252  int i2 = 7 - i1;
2253  int i3 = (process[7].id() > 0) ? 7 : 8;
2254  int i4 = 15 - i3;
2255  int i5 = (process[9].id() > 0) ? 9 : 10;
2256  int i6 = 19 - i5;
2257 
2258  // Set up four-products and internal products.
2259  setupProd( process, i1, i2, i3, i4, i5, i6);
2260 
2261  // tHat and uHat of fbar f -> W- W+ opposite to previous convention.
2262  double tHres = uH;
2263  double uHres = tH;
2264 
2265  // Couplings of incoming (anti)fermion.
2266  int idAbs = process[i1].idAbs();
2267  double ai = couplingsPtr->af(idAbs);
2268  double li = couplingsPtr->lf(idAbs);
2269  double ri = couplingsPtr->rf(idAbs);
2270 
2271  // gamma*/Z0 propagator/interference factor.
2272  double Zint = mZS * (sH - mZS) / (pow2(sH - mZS) + mwZS);
2273 
2274  // Combinations of couplings and kinematics (norm(x) = |x|^2).
2275  double dWW = (li * Zint + ai) / sH;
2276  double aWW = dWW + 0.5 * (ai + 1.) / tHres;
2277  double bWW = dWW + 0.5 * (ai - 1.) / uHres;
2278  double cWW = ri * Zint / sH;
2279  double fGK135 = norm( aWW * fGK( 1, 2, 3, 4, 5, 6)
2280  - bWW * fGK( 1, 2, 5, 6, 3, 4) );
2281  double fGK253 = norm( cWW * ( fGK( 2, 1, 5, 6, 3, 4)
2282  - fGK( 2, 1, 3, 4, 5, 6) ) );
2283  double xiT = xiGK( tHres, uHres);
2284  double xiU = xiGK( uHres, tHres);
2285  double xjTU = xjGK( tHres, uHres);
2286 
2287  // Weight and maximum weight.
2288  double wt = fGK135 + fGK253;
2289  double wtMax = 4. * s3 * s4
2290  * ( aWW * aWW * xiT + bWW * bWW * xiU - aWW * bWW * xjTU
2291  + cWW * cWW * (xiT + xiU - xjTU) );
2292 
2293  // Done.
2294  return wt / wtMax;
2295 }
2296 
2297 //==========================================================================
2298 
2299 // Sigma2ffbargmZggm class.
2300 // Collects common methods for f fbar -> gamma*/Z0 g/gamma and permutations.
2301 
2302 //--------------------------------------------------------------------------
2303 
2304 // Initialize process.
2305 
2306 void Sigma2ffbargmZggm::initProc() {
2307 
2308  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
2309  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
2310 
2311  // Store Z0 mass and width for propagator.
2312  mRes = particleDataPtr->m0(23);
2313  GammaRes = particleDataPtr->mWidth(23);
2314  m2Res = mRes*mRes;
2315  GamMRat = GammaRes / mRes;
2316  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
2317  * couplingsPtr->cos2thetaW());
2318 
2319  // Set pointer to particle properties and decay table.
2320  particlePtr = particleDataPtr->particleDataEntryPtr(23);
2321 
2322 }
2323 
2324 //--------------------------------------------------------------------------
2325 
2326 // Evaluate sum of flavour couplings times phase space.
2327 
2328 void Sigma2ffbargmZggm::flavSum() {
2329 
2330  // Coupling factors for Z0 subsystem.
2331  double alpSZ = couplingsPtr->alphaS(s3);
2332  double colQZ = 3. * (1. + alpSZ / M_PI);
2333 
2334  // Reset quantities to sum. Declare variables in loop.
2335  gamSum = 0.;
2336  intSum = 0.;
2337  resSum = 0.;
2338  int onMode;
2339  double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
2340 
2341  // Loop over all Z0 decay channels.
2342  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
2343  int idAbs = abs( particlePtr->channel(i).product(0) );
2344 
2345  // Only contributions from three fermion generations, except top.
2346  if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
2347  mf = particleDataPtr->m0(idAbs);
2348 
2349  // Check that above threshold. Phase space.
2350  if (m3 > 2. * mf + MASSMARGIN) {
2351  mr = pow2(mf / m3);
2352  betaf = sqrtpos(1. - 4. * mr);
2353  psvec = betaf * (1. + 2. * mr);
2354  psaxi = pow3(betaf);
2355 
2356  // Combine phase space with couplings.
2357  ef2 = couplingsPtr->ef2(idAbs) * psvec;
2358  efvf = couplingsPtr->efvf(idAbs) * psvec;
2359  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
2360  + couplingsPtr->af2(idAbs) * psaxi;
2361  colf = (idAbs < 6) ? colQZ : 1.;
2362 
2363  // Store sum of combinations. For outstate only open channels.
2364  onMode = particlePtr->channel(i).onMode();
2365  if (onMode == 1 || onMode == 2) {
2366  gamSum += colf * ef2;
2367  intSum += colf * efvf;
2368  resSum += colf * vf2af2;
2369  }
2370 
2371  // End loop over fermions.
2372  }
2373  }
2374  }
2375 
2376  // Done. Return values in gamSum, intSum and resSum.
2377 
2378 }
2379 
2380 //--------------------------------------------------------------------------
2381 
2382 // Calculate common parts of gamma/interference/Z0 propagator terms.
2383 
2384 void Sigma2ffbargmZggm::propTerm() {
2385 
2386  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2387  gamProp = 4. * alpEM / (3. * M_PI * s3);
2388  intProp = gamProp * 2. * thetaWRat * s3 * (s3 - m2Res)
2389  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2390  resProp = gamProp * pow2(thetaWRat * s3)
2391  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2392 
2393  // Optionally only keep gamma* or Z0 term.
2394  if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
2395  if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
2396 
2397 }
2398 
2399 //--------------------------------------------------------------------------
2400 
2401 // Evaluate weight for gamma*/Z0 decay angle.
2402 
2403 double Sigma2ffbargmZggm::weightDecay( Event& process, int iResBeg,
2404  int iResEnd) {
2405 
2406  // Z should sit in entry 5 and one more parton in entry 6.
2407  if (iResBeg != 5 || iResEnd != 6) return 1.;
2408 
2409  // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2410  // where f' fbar' come from gamma*/Z0 decay.
2411  int i1, i2;
2412  int i3 = (process[7].id() > 0) ? 7 : 8;
2413  int i4 = 15 - i3;
2414 
2415  // Order so that fbar(1) f(2) -> gamma*/Z0 g/gamma.
2416  if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2417  i1 = (process[3].id() < 0) ? 3 : 4;
2418  i2 = 7 - i1;
2419 
2420  // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) gamma*/Z0.
2421  } else if (process[3].idAbs() < 20) {
2422  i1 = (process[3].id() < 0) ? 3 : 6;
2423  i2 = 9 - i1;
2424  } else {
2425  i1 = (process[4].id() < 0) ? 4 : 6;
2426  i2 = 10 - i1;
2427  }
2428 
2429  // Charge/2, left- and righthanded couplings for in- and out-fermion.
2430  int id1Abs = process[i1].idAbs();
2431  double ei = 0.5 * couplingsPtr->ef(id1Abs);
2432  double li = couplingsPtr->lf(id1Abs);
2433  double ri = couplingsPtr->rf(id1Abs);
2434  int id3Abs = process[i3].idAbs();
2435  double ef = 0.5 * couplingsPtr->ef(id3Abs);
2436  double lf = couplingsPtr->lf(id3Abs);
2437  double rf = couplingsPtr->rf(id3Abs);
2438 
2439  // Combinations of left/right for in/out, gamma*/interference/Z0.
2440  double clilf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*lf
2441  + li*li * resProp * lf*lf;
2442  double clirf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*rf
2443  + li*li * resProp * rf*rf;
2444  double crilf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*lf
2445  + ri*ri * resProp * lf*lf;
2446  double crirf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*rf
2447  + ri*ri * resProp * rf*rf;
2448 
2449  // Evaluate four-vector products.
2450  double p13 = process[i1].p() * process[i3].p();
2451  double p14 = process[i1].p() * process[i4].p();
2452  double p23 = process[i2].p() * process[i3].p();
2453  double p24 = process[i2].p() * process[i4].p();
2454 
2455  // Calculate weight and its maximum.
2456  double wt = (clilf + crirf) * (p13*p13 + p24*p24)
2457  + (clirf + crilf) * (p14*p14 + p23*p23) ;
2458  double wtMax = (clilf + clirf + crilf + crirf)
2459  * (pow2(p13 + p14) + pow2(p23 + p24));
2460 
2461  // Done.
2462  return (wt / wtMax);
2463 
2464 }
2465 
2466 //==========================================================================
2467 
2468 // Sigma2qqbar2gmZg class.
2469 // Cross section for q qbar -> gamma*/Z0 g.
2470 
2471 //--------------------------------------------------------------------------
2472 
2473 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2474 
2475 void Sigma2qqbar2gmZg::sigmaKin() {
2476 
2477  // Cross section part common for all incoming flavours.
2478  sigma0 = (M_PI / sH2) * (alpEM * alpS)
2479  * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2480 
2481  // Calculate flavour sums for final state.
2482  flavSum();
2483 
2484  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2485  propTerm();
2486 
2487 }
2488 
2489 //--------------------------------------------------------------------------
2490 
2491 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2492 
2493 double Sigma2qqbar2gmZg::sigmaHat() {
2494 
2495  // Combine gamma, interference and Z0 parts.
2496  int idAbs = abs(id1);
2497  double sigma = sigma0
2498  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2499  + couplingsPtr->efvf(idAbs) * intProp * intSum
2500  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2501 
2502  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2503  sigma /= runBW3;
2504 
2505  // Answer.
2506  return sigma;
2507 
2508 }
2509 
2510 //--------------------------------------------------------------------------
2511 
2512 // Select identity, colour and anticolour.
2513 
2514 void Sigma2qqbar2gmZg::setIdColAcol() {
2515 
2516  // Flavours trivial.
2517  setId( id1, id2, 23, 21);
2518 
2519  // Colour flow topologies. Swap when antiquarks.
2520  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2521  if (id1 < 0) swapColAcol();
2522 
2523 }
2524 
2525 //==========================================================================
2526 
2527 // Sigma2qg2gmZq class.
2528 // Cross section for q g -> gamma*/Z0 q.
2529 
2530 //--------------------------------------------------------------------------
2531 
2532 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2533 
2534 void Sigma2qg2gmZq::sigmaKin() {
2535 
2536  // Cross section part common for all incoming flavours.
2537  sigma0 = (M_PI / sH2) * (alpEM * alpS)
2538  * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2539 
2540  // Calculate flavour sums for final state.
2541  flavSum();
2542 
2543  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2544  propTerm();
2545 
2546 }
2547 
2548 //--------------------------------------------------------------------------
2549 
2550 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2551 
2552 double Sigma2qg2gmZq::sigmaHat() {
2553 
2554  // Combine gamma, interference and Z0 parts.
2555  int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2556  double sigma = sigma0
2557  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2558  + couplingsPtr->efvf(idAbs) * intProp * intSum
2559  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2560 
2561  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2562  sigma /= runBW3;
2563 
2564  // Answer.
2565  return sigma;
2566 
2567 }
2568 
2569 //--------------------------------------------------------------------------
2570 
2571 // Select identity, colour and anticolour.
2572 
2573 void Sigma2qg2gmZq::setIdColAcol() {
2574 
2575  // Flavour set up for q g -> gamma*/Z0 q.
2576  int idq = (id2 == 21) ? id1 : id2;
2577  setId( id1, id2, 23, idq);
2578 
2579  // tH defined between f and f': must swap tHat <-> uHat if q g in.
2580  swapTU = (id2 == 21);
2581 
2582  // Colour flow topologies. Swap when antiquarks.
2583  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2584  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2585  if (idq < 0) swapColAcol();
2586 
2587 }
2588 
2589 //==========================================================================
2590 
2591 // Sigma2ffbar2gmZgm class.
2592 // Cross section for f fbar -> gamma*/Z0 gamma.
2593 
2594 //--------------------------------------------------------------------------
2595 
2596 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2597 
2598 void Sigma2ffbar2gmZgm::sigmaKin() {
2599 
2600  // Cross section part common for all incoming flavours.
2601  sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2602  * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2603 
2604  // Calculate flavour sums for final state.
2605  flavSum();
2606 
2607  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2608  propTerm();
2609 
2610 
2611 }
2612 
2613 //--------------------------------------------------------------------------
2614 
2615 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2616 
2617 double Sigma2ffbar2gmZgm::sigmaHat() {
2618 
2619  // Combine gamma, interference and Z0 parts.
2620  int idAbs = abs(id1);
2621  double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2622  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2623  + couplingsPtr->efvf(idAbs) * intProp * intSum
2624  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2625 
2626  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2627  sigma /= runBW3;
2628 
2629  // Colour factor. Answer.
2630  if (idAbs < 9) sigma /= 3.;
2631  return sigma;
2632 
2633 }
2634 
2635 //--------------------------------------------------------------------------
2636 
2637 // Select identity, colour and anticolour.
2638 
2639 void Sigma2ffbar2gmZgm::setIdColAcol() {
2640 
2641  // Flavours trivial.
2642  setId( id1, id2, 23, 22);
2643 
2644  // Colour flow topologies. Swap when antiquarks.
2645  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2646  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2647  if (id1 < 0) swapColAcol();
2648 
2649 }
2650 
2651 //==========================================================================
2652 
2653 // Sigma2fgm2gmZf class.
2654 // Cross section for f gamma -> gamma*/Z0 f'.
2655 
2656 //--------------------------------------------------------------------------
2657 
2658 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2659 
2660 void Sigma2fgm2gmZf::sigmaKin() {
2661 
2662  // Cross section part common for all incoming flavours.
2663  sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2664  * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (- sH * uH);
2665 
2666  // Calculate flavour sums for final state.
2667  flavSum();
2668 
2669  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2670  propTerm();
2671 
2672 }
2673 
2674 //--------------------------------------------------------------------------
2675 
2676 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2677 
2678 double Sigma2fgm2gmZf::sigmaHat() {
2679 
2680  // Combine gamma, interference and Z0 parts.
2681  int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
2682  double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2683  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2684  + couplingsPtr->efvf(idAbs) * intProp * intSum
2685  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2686 
2687  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2688  sigma /= runBW3;
2689 
2690  // Answer.
2691  return sigma;
2692 
2693 }
2694 
2695 //--------------------------------------------------------------------------
2696 
2697 // Select identity, colour and anticolour.
2698 
2699 void Sigma2fgm2gmZf::setIdColAcol() {
2700 
2701  // Flavour set up for q gamma -> gamma*/Z0 q.
2702  int idq = (id2 == 22) ? id1 : id2;
2703  setId( id1, id2, 23, idq);
2704 
2705  // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2706  swapTU = (id2 == 22);
2707 
2708  // Colour flow topologies. Swap when antiquarks.
2709  if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2710  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2711  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2712  if (idq < 0) swapColAcol();
2713 
2714 }
2715 
2716 //==========================================================================
2717 
2718 // Sigma2ffbarWggm class.
2719 // Collects common methods for f fbar -> W+- g/gamma and permutations.
2720 
2721 //--------------------------------------------------------------------------
2722 
2723 // Evaluate weight for W+- decay angle.
2724 
2725 double Sigma2ffbarWggm::weightDecay( Event& process, int iResBeg,
2726  int iResEnd) {
2727 
2728  // W should sit in entry 5 and one more parton in entry 6.
2729  if (iResBeg != 5 || iResEnd != 6) return 1.;
2730 
2731  // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2732  // where f' fbar' come from W+- decay.
2733  int i1, i2;
2734  int i3 = (process[7].id() > 0) ? 7 : 8;
2735  int i4 = 15 - i3;
2736 
2737  // Order so that fbar(1) f(2) -> W+- g/gamma.
2738  if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2739  i1 = (process[3].id() < 0) ? 3 : 4;
2740  i2 = 7 - i1;
2741 
2742  // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) W+-.
2743  } else if (process[3].idAbs() < 20) {
2744  i1 = (process[3].id() < 0) ? 3 : 6;
2745  i2 = 9 - i1;
2746  } else {
2747  i1 = (process[4].id() < 0) ? 4 : 6;
2748  i2 = 10 - i1;
2749  }
2750 
2751  // Evaluate four-vector products.
2752  double p13 = process[i1].p() * process[i3].p();
2753  double p14 = process[i1].p() * process[i4].p();
2754  double p23 = process[i2].p() * process[i3].p();
2755  double p24 = process[i2].p() * process[i4].p();
2756 
2757  // Calculate weight and its maximum.
2758  double wt = pow2(p13) + pow2(p24);
2759  double wtMax = pow2(p13 + p14) + pow2(p23 + p24);
2760 
2761  // Done.
2762  return (wt / wtMax);
2763 
2764 }
2765 
2766 //==========================================================================
2767 
2768 // Sigma2qqbar2Wg class.
2769 // Cross section for q qbar' -> W+- g.
2770 
2771 //--------------------------------------------------------------------------
2772 
2773 // Initialize process.
2774 
2775 void Sigma2qqbar2Wg::initProc() {
2776 
2777  // Secondary open width fractions, relevant for top (or heavier).
2778  openFracPos = particleDataPtr->resOpenFrac(24);
2779  openFracNeg = particleDataPtr->resOpenFrac(-24);
2780 
2781 }
2782 
2783 //--------------------------------------------------------------------------
2784 
2785 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2786 
2787 void Sigma2qqbar2Wg::sigmaKin() {
2788 
2789  // Cross section part common for all incoming flavours.
2790  sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2791  * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2792 
2793 }
2794 
2795 //--------------------------------------------------------------------------
2796 
2797 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2798 
2799 double Sigma2qqbar2Wg::sigmaHat() {
2800 
2801  // CKM factor. Secondary width for W+ or W-.
2802  double sigma = sigma0 * couplingsPtr->V2CKMid(abs(id1), abs(id2));
2803  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2804  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2805 
2806  // Answer.
2807  return sigma;
2808 
2809 }
2810 
2811 //--------------------------------------------------------------------------
2812 
2813 // Select identity, colour and anticolour.
2814 
2815 void Sigma2qqbar2Wg::setIdColAcol() {
2816 
2817  // Sign of outgoing W.
2818  int sign = 1 - 2 * (abs(id1)%2);
2819  if (id1 < 0) sign = -sign;
2820  setId( id1, id2, 24 * sign, 21);
2821 
2822  // Colour flow topologies. Swap when antiquarks.
2823  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2824  if (id1 < 0) swapColAcol();
2825 
2826 }
2827 
2828 //==========================================================================
2829 
2830 // Sigma2qg2Wq class.
2831 // Cross section for q g -> W+- q'.
2832 
2833 //--------------------------------------------------------------------------
2834 
2835 // Initialize process.
2836 
2837 void Sigma2qg2Wq::initProc() {
2838 
2839  // Secondary open width fractions, relevant for top (or heavier).
2840  openFracPos = particleDataPtr->resOpenFrac(24);
2841  openFracNeg = particleDataPtr->resOpenFrac(-24);
2842 
2843 }
2844 
2845 //--------------------------------------------------------------------------
2846 
2847 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2848 
2849 void Sigma2qg2Wq::sigmaKin() {
2850 
2851  // Cross section part common for all incoming flavours.
2852  sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2853  * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2854 
2855 }
2856 
2857 //--------------------------------------------------------------------------
2858 
2859 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2860 
2861 double Sigma2qg2Wq::sigmaHat() {
2862 
2863  // CKM factor. Secondary width for W+ or W-.
2864  int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2865  double sigma = sigma0 * couplingsPtr->V2CKMsum(idAbs);
2866  int idUp = (id2 == 21) ? id1 : id2;
2867  if (idAbs%2 == 1) idUp = -idUp;
2868  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2869 
2870  // Answer.
2871  return sigma;
2872 
2873 }
2874 
2875 //--------------------------------------------------------------------------
2876 
2877 // Select identity, colour and anticolour.
2878 
2879 void Sigma2qg2Wq::setIdColAcol() {
2880 
2881  // Sign of outgoing W. Flavour of outgoing quark.
2882  int idq = (id2 == 21) ? id1 : id2;
2883  int sign = 1 - 2 * (abs(idq)%2);
2884  if (idq < 0) sign = -sign;
2885  id4 = couplingsPtr->V2CKMpick(idq);
2886 
2887  // Flavour set up for q g -> W q.
2888  setId( id1, id2, 24 * sign, id4);
2889 
2890  // tH defined between f and f': must swap tHat <-> uHat if q g in.
2891  swapTU = (id2 == 21);
2892 
2893  // Colour flow topologies. Swap when antiquarks.
2894  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2895  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2896  if (idq < 0) swapColAcol();
2897 
2898 }
2899 
2900 //==========================================================================
2901 
2902 // Sigma2ffbar2Wgm class.
2903 // Cross section for f fbar' -> W+- gamma.
2904 
2905 //--------------------------------------------------------------------------
2906 
2907 // Initialize process.
2908 
2909 void Sigma2ffbar2Wgm::initProc() {
2910 
2911  // Secondary open width fractions, relevant for top (or heavier).
2912  openFracPos = particleDataPtr->resOpenFrac(24);
2913  openFracNeg = particleDataPtr->resOpenFrac(-24);
2914 
2915 }
2916 
2917 //--------------------------------------------------------------------------
2918 
2919 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2920 
2921 void Sigma2ffbar2Wgm::sigmaKin() {
2922 
2923  // Cross section part common for all incoming flavours.
2924  sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2925  * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2926 }
2927 
2928 //--------------------------------------------------------------------------
2929 
2930 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2931 
2932 double Sigma2ffbar2Wgm::sigmaHat() {
2933 
2934  // Extrafactor different for e nu and q qbar' instate.
2935  int id1Abs = abs(id1);
2936  int id2Abs = abs(id2);
2937  double chgUp = (id1Abs > 10) ? 0. : 2./3.;
2938  double sigma = sigma0 * pow2( chgUp - tH / (tH + uH) );
2939 
2940  // CKM and colour factors. Secondary width for W+ or W-.
2941  if (id1Abs < 9) sigma *= couplingsPtr->V2CKMid(id1Abs, id2Abs) / 3.;
2942  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2943  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2944 
2945  // Answer.
2946  return sigma;
2947 
2948 }
2949 
2950 //--------------------------------------------------------------------------
2951 
2952 // Select identity, colour and anticolour.
2953 
2954 void Sigma2ffbar2Wgm::setIdColAcol() {
2955 
2956  // Sign of outgoing W.
2957  int sign = 1 - 2 * (abs(id1)%2);
2958  if (id1 < 0) sign = -sign;
2959  setId( id1, id2, 24 * sign, 22);
2960 
2961  // tH defined between (f,W-) or (fbar',W+).
2962  swapTU = (sign * id1 > 0);
2963 
2964  // Colour flow topologies. Swap when antiquarks.
2965  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2966  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2967  if (id1 < 0) swapColAcol();
2968 
2969 }
2970 
2971 //==========================================================================
2972 
2973 // Sigma2fgm2Wf class.
2974 // Cross section for f gamma -> W+- f'.
2975 
2976 //--------------------------------------------------------------------------
2977 
2978 // Initialize process.
2979 
2980 void Sigma2fgm2Wf::initProc() {
2981 
2982  // Secondary open width fractions, relevant for top (or heavier).
2983  openFracPos = particleDataPtr->resOpenFrac(24);
2984  openFracNeg = particleDataPtr->resOpenFrac(-24);
2985 
2986 }
2987 
2988 //--------------------------------------------------------------------------
2989 
2990 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2991 
2992 void Sigma2fgm2Wf::sigmaKin() {
2993 
2994  // Cross section part common for all incoming flavours.
2995  sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2996  * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (pT2 * s3 - sH * uH);
2997 
2998 }
2999 
3000 //--------------------------------------------------------------------------
3001 
3002 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
3003 
3004 double Sigma2fgm2Wf::sigmaHat() {
3005 
3006  // Extrafactor dependent on charge of incoming fermion.
3007  int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
3008  double charge = (idAbs > 10) ? 1. : ( (idAbs%2 == 1) ? 1./3. : 2./3. );
3009  double sigma = sigma0 * pow2( charge - sH / (sH + uH) );
3010 
3011  // CKM factor. Secondary width for W+ or W-.
3012  sigma *= couplingsPtr->V2CKMsum(idAbs);
3013  int idUp = (id2 == 22) ? id1 : id2;
3014  if (idAbs%2 == 1) idUp = -idUp;
3015  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
3016 
3017  // Answer.
3018  return sigma;
3019 
3020 }
3021 
3022 //--------------------------------------------------------------------------
3023 
3024 // Select identity, colour and anticolour.
3025 
3026 void Sigma2fgm2Wf::setIdColAcol() {
3027 
3028  // Sign of outgoing W. Flavour of outgoing fermion.
3029  int idq = (id2 == 22) ? id1 : id2;
3030  int sign = 1 - 2 * (abs(idq)%2);
3031  if (idq < 0) sign = -sign;
3032  id4 = couplingsPtr->V2CKMpick(idq);
3033 
3034  // Flavour set up for q gamma -> W q.
3035  setId( id1, id2, 24 * sign, id4);
3036 
3037  // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
3038  swapTU = (id2 == 22);
3039 
3040  // Colour flow topologies. Swap when antiquarks.
3041  if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
3042  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
3043  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
3044  if (idq < 0) swapColAcol();
3045 
3046 }
3047 
3048 //==========================================================================
3049 
3050 // Sigma2gmgm2ffbar class.
3051 // Cross section for gamma gamma -> l lbar.
3052 
3053 //--------------------------------------------------------------------------
3054 
3055 // Initialize process wrt flavour.
3056 
3057 void Sigma2gmgm2ffbar::initProc() {
3058 
3059  // Process name.
3060  nameSave = "gamma gamma -> f fbar";
3061  if (idNew == 1) nameSave = "gamma gamma -> q qbar (uds)";
3062  if (idNew == 4) nameSave = "gamma gamma -> c cbar";
3063  if (idNew == 5) nameSave = "gamma gamma -> b bbar";
3064  if (idNew == 6) nameSave = "gamma gamma -> t tbar";
3065  if (idNew == 11) nameSave = "gamma gamma -> e+ e-";
3066  if (idNew == 13) nameSave = "gamma gamma -> mu+ mu-";
3067  if (idNew == 15) nameSave = "gamma gamma -> tau+ tau-";
3068 
3069  // Generate massive phase space, except for u+d+s.
3070  idMass = 0;
3071  if (idNew > 3) idMass = idNew;
3072 
3073  // Charge and colour factor.
3074  ef4 = 1.;
3075  if (idNew == 1) ef4 = 3. * (pow4(2./3.) + 2. * pow4(1./3.));
3076  if (idNew == 4 || idNew == 6) ef4 = 3. * pow4(2./3.);
3077  if (idNew == 5) ef4 = 3. * pow4(1./3.);
3078 
3079  // Secondary open width fraction.
3080  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
3081 
3082 }
3083 
3084 //--------------------------------------------------------------------------
3085 
3086 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
3087 
3088 void Sigma2gmgm2ffbar::sigmaKin() {
3089 
3090  // Pick current flavour for u+d+s mix by e_q^4 weights.
3091  if (idNew == 1) {
3092  double rId = 18. * rndmPtr->flat();
3093  idNow = 1;
3094  if (rId > 1.) idNow = 2;
3095  if (rId > 17.) idNow = 3;
3096  s34Avg = pow2(particleDataPtr->m0(idNow));
3097  } else {
3098  idNow = idNew;
3099  s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
3100  }
3101 
3102  // Modified Mandelstam variables for massive kinematics with m3 = m4.
3103  double tHQ = -0.5 * (sH - tH + uH);
3104  double uHQ = -0.5 * (sH + tH - uH);
3105  double tHQ2 = tHQ * tHQ;
3106  double uHQ2 = uHQ * uHQ;
3107 
3108  // Calculate kinematics dependence.
3109  if (sH < 4. * s34Avg) sigTU = 0.;
3110  else sigTU = 2. * (tHQ * uHQ - s34Avg * sH)
3111  * (tHQ2 + uHQ2 + 2. * s34Avg * sH) / (tHQ2 * uHQ2);
3112 
3113  // Answer.
3114  sigma = (M_PI / sH2) * pow2(alpEM) * ef4 * sigTU * openFracPair;
3115 
3116 }
3117 
3118 //--------------------------------------------------------------------------
3119 
3120 // Select identity, colour and anticolour.
3121 
3122 void Sigma2gmgm2ffbar::setIdColAcol() {
3123 
3124  // Flavours are trivial.
3125  setId( id1, id2, idNow, -idNow);
3126 
3127  // Colour flow in singlet state.
3128  if (idNow < 10) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
3129  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
3130 
3131 }
3132 
3133 //==========================================================================
3134 
3135 // Sigma2ggm2qqbar class.
3136 // Cross section for g gamma -> q qbar.
3137 
3138 //--------------------------------------------------------------------------
3139 
3140 // Initialize process wrt flavour.
3141 
3142 void Sigma2ggm2qqbar::initProc() {
3143 
3144  // Initialize the process name according to the incoming partons.
3145  if (inFluxSave == "ggm") {
3146  nameSave = "g gamma -> q qbar";
3147  if (idNew == 1) nameSave = "g gamma -> q qbar (uds)";
3148  if (idNew == 4) nameSave = "g gamma -> c cbar";
3149  if (idNew == 5) nameSave = "g gamma -> b bbar";
3150  if (idNew == 6) nameSave = "g gamma -> t tbar";
3151  } else if (inFluxSave == "gmg") {
3152  nameSave = "gamma g -> q qbar";
3153  if (idNew == 1) nameSave = "gamma g -> q qbar (uds)";
3154  if (idNew == 4) nameSave = "gamma g -> c cbar";
3155  if (idNew == 5) nameSave = "gamma g -> b bbar";
3156  if (idNew == 6) nameSave = "gamma g -> t tbar";
3157  }
3158 
3159  // Generate massive phase space, except for u+d+s.
3160  idMass = 0;
3161  if (idNew > 3) idMass = idNew;
3162 
3163  // Charge and colour factor.
3164  ef2 = 1.;
3165  if (idNew == 1) ef2 = pow2(2./3.) + 2. * pow2(1./3.);
3166  if (idNew == 4 || idNew == 6) ef2 = pow2(2./3.);
3167  if (idNew == 5) ef2 = pow2(1./3.);
3168 
3169  // Secondary open width fraction.
3170  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
3171 
3172 }
3173 
3174 //--------------------------------------------------------------------------
3175 
3176 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
3177 
3178 void Sigma2ggm2qqbar::sigmaKin() {
3179 
3180  // Pick current flavour for u+d+s mix by e_q^2 weights.
3181  if (idNew == 1) {
3182  double rId = 6. * rndmPtr->flat();
3183  idNow = 1;
3184  if (rId > 1.) idNow = 2;
3185  if (rId > 5.) idNow = 3;
3186  s34Avg = pow2(particleDataPtr->m0(idNow));
3187  } else {
3188  idNow = idNew;
3189  s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
3190  }
3191 
3192  // Modified Mandelstam variables for massive kinematics with m3 = m4.
3193  double tHQ = -0.5 * (sH - tH + uH);
3194  double uHQ = -0.5 * (sH + tH - uH);
3195  double tHQ2 = tHQ * tHQ;
3196  double uHQ2 = uHQ * uHQ;
3197 
3198  // Calculate kinematics dependence.
3199  if (sH < 4. * s34Avg) sigTU = 0.;
3200  else sigTU = (tHQ * uHQ - s34Avg * sH)
3201  * (tHQ2 + uHQ2 + 2. * s34Avg * sH) / (tHQ2 * uHQ2);
3202 
3203  // Answer.
3204  sigma = (M_PI/sH2) * alpS * alpEM * ef2 * sigTU * openFracPair;
3205 
3206 }
3207 
3208 //--------------------------------------------------------------------------
3209 
3210 // Select identity, colour and anticolour.
3211 
3212 void Sigma2ggm2qqbar::setIdColAcol() {
3213 
3214  // Construct outgoing flavours.
3215  setId( id1, id2, idNow, -idNow);
3216 
3217  // Colour flow topology. Swap if first is gamma.
3218  setColAcol( 1, 2, 0, 0, 1, 0, 0, 2);
3219  if (id1 == 22) setColAcol( 0, 0, 1, 2, 1, 0, 0, 2);
3220 
3221 }
3222 
3223 //==========================================================================
3224 
3225 // Sigma2qgm2qg class.
3226 // Cross section for q gamma -> q g.
3227 
3228 //--------------------------------------------------------------------------
3229 
3230 // Initialize process wrt incoming particles.
3231 
3232 void Sigma2qgm2qg::initProc() {
3233 
3234  // Initialize the process name according to the incoming partons.
3235  if (inFluxSave == "qgm") nameSave = "q gamma -> q g (udscb)";
3236  if (inFluxSave == "gmq") nameSave = "gamma q -> q g (udscb)";
3237 
3238 }
3239 
3240 //--------------------------------------------------------------------------
3241 
3242 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
3243 
3244 void Sigma2qgm2qg::sigmaKin() {
3245 
3246  // Calculate kinematics dependence.
3247  sigUS = (8./3.) * (sH2 + uH2) / (-sH * uH);
3248 
3249  // Answer.
3250  sigma0 = (M_PI/sH2) * alpS * alpEM * sigUS;
3251 
3252 }
3253 
3254 //--------------------------------------------------------------------------
3255 
3256 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
3257 
3258 double Sigma2qgm2qg::sigmaHat() {
3259 
3260  // Incoming flavour gives charge factor.
3261  int idNow = (id2 == 22) ? id1 : id2;
3262  double eNow = couplingsPtr->ef( abs(idNow) );
3263  return sigma0 * pow2(eNow);
3264 
3265 }
3266 
3267 //--------------------------------------------------------------------------
3268 
3269 // Select identity, colour and anticolour.
3270 
3271 void Sigma2qgm2qg::setIdColAcol() {
3272 
3273  // Construct outgoing flavours.
3274  id3 = (id1 == 22) ? 21 : id1;
3275  id4 = (id2 == 22) ? 21 : id2;
3276  setId( id1, id2, id3, id4);
3277 
3278  // Colour flow topology. Swap if first is gamma, or when antiquark.
3279  setColAcol( 1, 0, 0, 0, 2, 0, 1, 2);
3280  if (id1 == 22) swapCol1234();
3281  if (id1 < 0 || id2 < 0) swapColAcol();
3282 
3283 }
3284 
3285 //==========================================================================
3286 
3287 // Sigma2qgm2qgm class.
3288 // Cross section for q gamma -> q gamma.
3289 
3290 //--------------------------------------------------------------------------
3291 
3292 // Initialize process wrt incoming particles.
3293 
3294 void Sigma2qgm2qgm::initProc() {
3295 
3296  // Initialize the process name according to the incoming partons.
3297  if (inFluxSave == "qgm") nameSave = "q gamma -> q gamma (udscb)";
3298  if (inFluxSave == "gmq") nameSave = "gamma q -> q gamma (udscb)";
3299 
3300 }
3301 
3302 //--------------------------------------------------------------------------
3303 
3304 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
3305 
3306 void Sigma2qgm2qgm::sigmaKin() {
3307 
3308  // Calculate kinematics dependence.
3309  sigUS = 2. * (sH2 + uH2) / (-sH * uH);
3310 
3311  // Answer.
3312  sigma0 = (M_PI/sH2) * pow2(alpEM) * sigUS;
3313 
3314 }
3315 
3316 //--------------------------------------------------------------------------
3317 
3318 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
3319 
3320 double Sigma2qgm2qgm::sigmaHat() {
3321 
3322  // Incoming flavour gives charge factor.
3323  int idNow = (id2 == 22) ? id1 : id2;
3324  double eNow = couplingsPtr->ef( abs(idNow) );
3325  return sigma0 * pow4(eNow);
3326 
3327 }
3328 
3329 //--------------------------------------------------------------------------
3330 
3331 // Select identity, colour and anticolour.
3332 
3333 void Sigma2qgm2qgm::setIdColAcol() {
3334 
3335  // Construct outgoing flavours.
3336  id3 = (id1 == 22) ? 22 : id1;
3337  id4 = (id2 == 22) ? 22 : id2;
3338  setId( id1, id2, id3, id4);
3339 
3340  // Colour flow topology. Swap if first is gamma, or when antiquark.
3341  if (id2 == 22) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
3342  if (id1 == 22) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
3343  if (id1 < 0 || id2 < 0) swapColAcol();
3344 
3345 }
3346 
3347 //==========================================================================
3348 
3349 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5