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) 2012 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) for the
7 // electroweak simulation classes (including photon processes).
8 
9 #include "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 
908 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
909 
910 void Sigma2ffbar2ffbarsgm::sigmaKin() {
911 
912  // Pick new flavour. Allow three leptons and five quarks.
913  double colQ = 1. + (alpS / M_PI);
914  double flavWt = 3. + colQ * 11. / 3.;
915  double flavRndm = rndmPtr->flat() * flavWt;
916  if (flavRndm < 3.) {
917  if (flavRndm < 1.) idNew = 11;
918  else if (flavRndm < 2.) idNew = 13;
919  else idNew = 15;
920  } else {
921  flavRndm = 3. * (flavRndm - 3.) / colQ;
922  if (flavRndm < 4.) idNew = 2;
923  else if (flavRndm < 8.) idNew = 4;
924  else if (flavRndm < 9.) idNew = 1;
925  else if (flavRndm < 10.) idNew = 3;
926  else idNew = 5;
927  }
928  double mNew = particleDataPtr->m0(idNew);
929  double m2New = mNew*mNew;
930 
931  // Calculate kinematics dependence. Give correct mass factors for
932  // tHat, uHat defined as if massless kinematics, d(sigma)/d(Omega)
933  // = beta (1 + cos^2(theta) + (1 - beta^2) sin^2(theta)).
934  // Special case related to phase space form in multiparton interactions.
935  double sigS = 0.;
936  if (sH > 4. * m2New) {
937  double beta = sqrt(1. - 4. * m2New / sH);
938  sigS = beta * (2.* (tH2 + uH2) + 4. * (1. - beta * beta) * tH * uH)
939  / sH2;
940  }
941 
942  // Answer is proportional to number of outgoing flavours.
943  sigma0 = (M_PI/sH2) * pow2(alpEM) * sigS * flavWt;
944 
945 }
946 
947 //--------------------------------------------------------------------------
948 
949 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
950 
951 double Sigma2ffbar2ffbarsgm::sigmaHat() {
952 
953  // Charge and colour factors.
954  double eNow = couplingsPtr->ef( abs(id1) );
955  double sigma = sigma0 * pow2(eNow);
956  if (abs(id1) < 9) sigma /= 3.;
957 
958  // Answer.
959  return sigma;
960 
961 }
962 
963 //--------------------------------------------------------------------------
964 
965 // Select identity, colour and anticolour.
966 
967 void Sigma2ffbar2ffbarsgm::setIdColAcol() {
968 
969  // Set outgoing flavours.
970  id3 = (id1 > 0) ? idNew : -idNew;
971  setId( id1, id2, id3, -id3);
972 
973  // Colour flow topologies. Swap when antiquarks.
974  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
975  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
976  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
977  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
978  if (id1 < 0) swapColAcol();
979 
980 }
981 
982 //==========================================================================
983 
984 // Sigma2ffbar2FFbarsgmZ class.
985 // Cross section f fbar -> gamma*/Z0 -> F Fbar.
986 
987 //--------------------------------------------------------------------------
988 
989 // Initialize process.
990 
991 void Sigma2ffbar2FFbarsgmZ::initProc() {
992 
993  // Process name.
994  nameSave = "f fbar -> F Fbar (s-channel gamma*/Z0)";
995  if (idNew == 4) nameSave = "f fbar -> c cbar (s-channel gamma*/Z0)";
996  if (idNew == 5) nameSave = "f fbar -> b bbar (s-channel gamma*/Z0)";
997  if (idNew == 6) nameSave = "f fbar -> t tbar (s-channel gamma*/Z0)";
998  if (idNew == 7) nameSave = "f fbar -> b' b'bar (s-channel gamma*/Z0)";
999  if (idNew == 8) nameSave = "f fbar -> t' t'bar (s-channel gamma*/Z0)";
1000  if (idNew == 15) nameSave = "f fbar -> tau+ tau- (s-channel gamma*/Z0)";
1001  if (idNew == 17) nameSave = "f fbar -> tau'+ tau'- (s-channel gamma*/Z0)";
1002  if (idNew == 18)
1003  nameSave = "f fbar -> nu'_tau nu'bar_tau (s-channel gamma*/Z0)";
1004 
1005  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1006  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
1007 
1008  // Store Z0 mass and width for propagator.
1009  mRes = particleDataPtr->m0(23);
1010  GammaRes = particleDataPtr->mWidth(23);
1011  m2Res = mRes*mRes;
1012  GamMRat = GammaRes / mRes;
1013  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1014  * couplingsPtr->cos2thetaW());
1015 
1016  // Store couplings of F.
1017  ef = couplingsPtr->ef(idNew);
1018  vf = couplingsPtr->vf(idNew);
1019  af = couplingsPtr->af(idNew);
1020 
1021  // Secondary open width fraction, relevant for top (or heavier).
1022  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
1023 
1024 }
1025 
1026 //--------------------------------------------------------------------------
1027 
1028 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1029 
1030 void Sigma2ffbar2FFbarsgmZ::sigmaKin() {
1031 
1032  // Check that above threshold.
1033  isPhysical = true;
1034  if (mH < m3 + m4 + MASSMARGIN) {
1035  isPhysical = false;
1036  return;
1037  }
1038 
1039  // Define average F, Fbar mass so same beta. Phase space.
1040  double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1041  mr = s34Avg / sH;
1042  betaf = sqrtpos(1. - 4. * mr);
1043 
1044  // Final-state colour factor.
1045  double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) : 1.;
1046 
1047  // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1048  cosThe = (tH - uH) / (betaf * sH);
1049 
1050  // Calculate prefactors for gamma/interference/Z0 cross section terms.
1051  gamProp = colF * M_PI * pow2(alpEM) / sH2;
1052  intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1053  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1054  resProp = gamProp * pow2(thetaWRat * sH)
1055  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1056 
1057  // Optionally only keep gamma* or Z0 term.
1058  if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1059  if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1060 
1061 }
1062 
1063 //--------------------------------------------------------------------------
1064 
1065 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1066 
1067 double Sigma2ffbar2FFbarsgmZ::sigmaHat() {
1068 
1069  // Fail if below threshold.
1070  if (!isPhysical) return 0.;
1071 
1072  // Couplings for in-flavours.
1073  int idAbs = abs(id1);
1074  double ei = couplingsPtr->ef(idAbs);
1075  double vi = couplingsPtr->vf(idAbs);
1076  double ai = couplingsPtr->af(idAbs);
1077 
1078  // Coefficients of angular expression.
1079  double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
1080  + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
1081  double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
1082  + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
1083  double coefAsym = betaf * ( ei * ai * intProp * ef * af
1084  + 4. * vi * ai * resProp * vf * af );
1085 
1086  // Combine gamma, interference and Z0 parts.
1087  double sigma = coefTran * (1. + pow2(cosThe))
1088  + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
1089 
1090  // Top: corrections for closed decay channels.
1091  sigma *= openFracPair;
1092 
1093  // Initial-state colour factor. Answer.
1094  if (idAbs < 9) sigma /= 3.;
1095  return sigma;
1096 
1097 }
1098 
1099 //--------------------------------------------------------------------------
1100 
1101 // Select identity, colour and anticolour.
1102 
1103 void Sigma2ffbar2FFbarsgmZ::setIdColAcol() {
1104 
1105  // Set outgoing flavours.
1106  id3 = (id1 > 0) ? idNew : -idNew;
1107  setId( id1, id2, id3, -id3);
1108 
1109  // Colour flow topologies. Swap when antiquarks.
1110  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1111  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1112  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1113  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1114  if (id1 < 0) swapColAcol();
1115 
1116 }
1117 
1118 //--------------------------------------------------------------------------
1119 
1120 // Evaluate weight for decay angles of W in top decay.
1121 
1122 double Sigma2ffbar2FFbarsgmZ::weightDecay( Event& process, int iResBeg,
1123  int iResEnd) {
1124 
1125  // For top decay hand over to standard routine, else done.
1126  if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1127  return weightTopDecay( process, iResBeg, iResEnd);
1128  else return 1.;
1129 
1130 }
1131 
1132 //==========================================================================
1133 
1134 // Sigma2ffbar2FfbarsW class.
1135 // Cross section f fbar' -> W+- -> F fbar".
1136 
1137 //--------------------------------------------------------------------------
1138 
1139 // Initialize process.
1140 
1141 void Sigma2ffbar2FfbarsW::initProc() {
1142 
1143  // Process name.
1144  nameSave = "f fbar -> F fbar (s-channel W+-)";
1145  if (idNew == 4) nameSave = "f fbar -> c qbar (s-channel W+-)";
1146  if (idNew == 5) nameSave = "f fbar -> b qbar (s-channel W+-)";
1147  if (idNew == 6) nameSave = "f fbar -> t qbar (s-channel W+-)";
1148  if (idNew == 7) nameSave = "f fbar -> b' qbar (s-channel W+-)";
1149  if (idNew == 8) nameSave = "f fbar -> t' qbar (s-channel W+-)";
1150  if (idNew == 7 && idNew2 == 6)
1151  nameSave = "f fbar -> b' tbar (s-channel W+-)";
1152  if (idNew == 8 && idNew2 == 7)
1153  nameSave = "f fbar -> t' b'bar (s-channel W+-)";
1154  if (idNew == 15 || idNew == 16)
1155  nameSave = "f fbar -> tau nu_taubar (s-channel W+-)";
1156  if (idNew == 17 || idNew == 18)
1157  nameSave = "f fbar -> tau' nu'_taubar (s-channel W+-)";
1158 
1159  // Store W+- mass and width for propagator.
1160  mRes = particleDataPtr->m0(24);
1161  GammaRes = particleDataPtr->mWidth(24);
1162  m2Res = mRes*mRes;
1163  GamMRat = GammaRes / mRes;
1164  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1165 
1166  // For t/t' want to use at least b mass.
1167  idPartner = idNew2;
1168  if ( (idNew == 6 || idNew == 8) && idNew2 == 0 ) idPartner = 5;
1169 
1170  // Sum of CKM weights for quarks.
1171  V2New = (idNew < 9) ? couplingsPtr->V2CKMsum(idNew) : 1.;
1172  if (idNew2 != 0) V2New = couplingsPtr->V2CKMid(idNew, idNew2);
1173 
1174  // Secondary open width fractions, relevant for top or heavier.
1175  openFracPos = particleDataPtr->resOpenFrac( idNew, -idNew2);
1176  openFracNeg = particleDataPtr->resOpenFrac(-idNew, idNew2);
1177 
1178 }
1179 
1180 //--------------------------------------------------------------------------
1181 
1182 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1183 
1184 void Sigma2ffbar2FfbarsW::sigmaKin() {
1185 
1186  // Check that above threshold.
1187  isPhysical = true;
1188  if (mH < m3 + m4 + MASSMARGIN) {
1189  isPhysical = false;
1190  return;
1191  }
1192 
1193  // Phase space factors.
1194  double mr1 = s3 / sH;
1195  double mr2 = s4 / sH;
1196  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
1197 
1198  // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1199  double cosThe = (tH - uH) / (betaf * sH);
1200 
1201  // Set up Breit-Wigner and in- and out-widths.
1202  double sigBW = 9. * M_PI * pow2(alpEM * thetaWRat)
1203  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1204 
1205  // Initial-state colour factor.
1206  double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) * V2New : 1.;
1207 
1208  // Angular dependence.
1209  double wt = pow2(1. + betaf * cosThe) - pow2(mr1 - mr2);
1210 
1211  // Temporary answer.
1212  sigma0 = sigBW * colF * wt;
1213 
1214 }
1215 
1216 //--------------------------------------------------------------------------
1217 
1218 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1219 
1220 double Sigma2ffbar2FfbarsW::sigmaHat() {
1221 
1222  // Fail if below threshold.
1223  if (!isPhysical) return 0.;
1224 
1225  // CKM and colour factors.
1226  double sigma = sigma0;
1227  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1228 
1229  // Correction for secondary width in top (or heavier) decay.
1230  int idSame = ((abs(id1) + idNew)%2 == 0) ? id1 : id2;
1231  sigma *= (idSame > 0) ? openFracPos : openFracNeg;
1232 
1233  // Answer.
1234  return sigma;
1235 
1236 }
1237 
1238 //--------------------------------------------------------------------------
1239 
1240 // Select identity, colour and anticolour.
1241 
1242 void Sigma2ffbar2FfbarsW::setIdColAcol() {
1243 
1244  // Set outgoing flavours.
1245  id3 = idNew;
1246  id4 = (idNew2 != 0) ? idNew2 : couplingsPtr->V2CKMpick(idNew);
1247  if (idNew%2 == 0) {
1248  int idInUp = (abs(id1)%2 == 0) ? id1 : id2;
1249  if (idInUp > 0) id4 = -id4;
1250  else id3 = -id3;
1251  } else {
1252  int idInDn = (abs(id1)%2 == 1) ? id1 : id2;
1253  if (idInDn > 0) id4 = -id4;
1254  else id3 = -id3;
1255  }
1256  setId( id1, id2, id3, id4);
1257 
1258  // Swap tHat and uHat for fbar' f -> F f".
1259  if (id1 * id3 < 0) swapTU = true;
1260 
1261  // Colour flow topologies. Swap when antiquarks.
1262  if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1263  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1264  else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1265  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1266  if (id1 < 0) swapCol12();
1267  if (id3 < 0) swapCol34();
1268 
1269 }
1270 
1271 //--------------------------------------------------------------------------
1272 
1273 // Evaluate weight for decay angles of W in top decay.
1274 
1275 double Sigma2ffbar2FfbarsW::weightDecay( Event& process, int iResBeg,
1276  int iResEnd) {
1277 
1278  // For top decay hand over to standard routine, else done.
1279  if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1280  return weightTopDecay( process, iResBeg, iResEnd);
1281  else return 1.;
1282 
1283 }
1284 
1285 //==========================================================================
1286 
1287 // Sigma2ffbargmZWgmZW class.
1288 // Collects common methods for f fbar -> gamma*/Z0/W+- gamma*/Z0/W-+.
1289 
1290 //--------------------------------------------------------------------------
1291 
1292 // Calculate and store internal products.
1293 
1294 void Sigma2ffbargmZWgmZW::setupProd( Event& process, int i1, int i2,
1295  int i3, int i4, int i5, int i6) {
1296 
1297  // Store incoming and outgoing momenta,
1298  pRot[1] = process[i1].p();
1299  pRot[2] = process[i2].p();
1300  pRot[3] = process[i3].p();
1301  pRot[4] = process[i4].p();
1302  pRot[5] = process[i5].p();
1303  pRot[6] = process[i6].p();
1304 
1305  // Do random rotation to avoid accidental zeroes in HA expressions.
1306  bool smallPT = false;
1307  do {
1308  smallPT = false;
1309  double thetaNow = acos(2. * rndmPtr->flat() - 1.);
1310  double phiNow = 2. * M_PI * rndmPtr->flat();
1311  for (int i = 1; i <= 6; ++i) {
1312  pRot[i].rot( thetaNow, phiNow);
1313  if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
1314  }
1315  } while (smallPT);
1316 
1317  // Calculate internal products.
1318  for (int i = 1; i < 6; ++i) {
1319  for (int j = i + 1; j <= 6; ++j) {
1320  hA[i][j] =
1321  sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
1322  / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
1323  - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
1324  / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
1325  hC[i][j] = conj( hA[i][j] );
1326  if (i <= 2) {
1327  hA[i][j] *= complex( 0., 1.);
1328  hC[i][j] *= complex( 0., 1.);
1329  }
1330  hA[j][i] = - hA[i][j];
1331  hC[j][i] = - hC[i][j];
1332  }
1333  }
1334 
1335 }
1336 
1337 //--------------------------------------------------------------------------
1338 
1339 // Evaluate the F function of Gunion and Kunszt.
1340 
1341 complex Sigma2ffbargmZWgmZW::fGK(int j1, int j2, int j3, int j4, int j5,
1342  int j6) {
1343 
1344  return 4. * hA[j1][j3] * hC[j2][j6]
1345  * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
1346 
1347 }
1348 
1349 //--------------------------------------------------------------------------
1350 
1351 // Evaluate the Xi function of Gunion and Kunszt.
1352 
1353 double Sigma2ffbargmZWgmZW::xiGK( double tHnow, double uHnow) {
1354 
1355  return - 4. * s3 * s4 + tHnow * (3. * tHnow + 4. * uHnow)
1356  + tHnow * tHnow * ( tHnow * uHnow / (s3 * s4)
1357  - 2. * (1. / s3 + 1./s4) * (tHnow + uHnow)
1358  + 2. * (s3 / s4 + s4 / s3) );
1359 
1360 }
1361 
1362 //--------------------------------------------------------------------------
1363 
1364 // Evaluate the Xj function of Gunion and Kunszt.
1365 
1366 double Sigma2ffbargmZWgmZW::xjGK( double tHnow, double uHnow) {
1367 
1368  return 8. * pow2(s3 + s4) - 8. * (s3 + s4) * (tHnow + uHnow)
1369  - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
1370  / (s3 * s4) - 2. * (1. / s3 + 1. / s4) * (tHnow + uHnow)
1371  + 2. * (s3 / s4 + s4 / s3) );
1372 
1373 }
1374 
1375 //==========================================================================
1376 
1377 // Sigma2ffbar2gmZgmZ class.
1378 // Cross section for f fbar -> gamma*/Z0 gamma*/Z0 (f is quark or lepton).
1379 
1380 //--------------------------------------------------------------------------
1381 
1382 // Initialize process.
1383 
1384 void Sigma2ffbar2gmZgmZ::initProc() {
1385 
1386  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1387  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
1388 
1389  // Store Z0 mass and width for propagator.
1390  mRes = particleDataPtr->m0(23);
1391  GammaRes = particleDataPtr->mWidth(23);
1392  m2Res = mRes*mRes;
1393  GamMRat = GammaRes / mRes;
1394  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1395  * couplingsPtr->cos2thetaW());
1396 
1397  // Set pointer to particle properties and decay table.
1398  particlePtr = particleDataPtr->particleDataEntryPtr(23);
1399 
1400 }
1401 
1402 //--------------------------------------------------------------------------
1403 
1404 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1405 
1406 void Sigma2ffbar2gmZgmZ::sigmaKin() {
1407 
1408  // Cross section part common for all incoming flavours.
1409  sigma0 = (M_PI / sH2) * pow2(alpEM) * 0.5
1410  * ( (tH2 + uH2 + 2. * (s3 + s4) * sH) / (tH * uH)
1411  - s3 * s4 * (1./tH2 + 1./uH2) );
1412 
1413  // Common coupling factors at the resonance masses
1414  double alpEM3 = couplingsPtr->alphaEM(s3);
1415  double alpS3 = couplingsPtr->alphaS(s3);
1416  double colQ3 = 3. * (1. + alpS3 / M_PI);
1417  double alpEM4 = couplingsPtr->alphaEM(s4);
1418  double alpS4 = couplingsPtr->alphaS(s4);
1419  double colQ4 = 3. * (1. + alpS4 / M_PI);
1420 
1421  // Reset quantities to sum. Declare variables in loop.
1422  gamSum3 = 0.;
1423  intSum3 = 0.;
1424  resSum3 = 0.;
1425  gamSum4 = 0.;
1426  intSum4 = 0.;
1427  resSum4 = 0.;
1428  int onMode;
1429  double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
1430 
1431  // Loop over all Z0 decay channels.
1432  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1433  int idAbs = abs( particlePtr->channel(i).product(0) );
1434 
1435  // Only contributions from three fermion generations, except top.
1436  if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
1437  mf = particleDataPtr->m0(idAbs);
1438  onMode = particlePtr->channel(i).onMode();
1439 
1440  // First Z0: check that above threshold. Phase space.
1441  if (m3 > 2. * mf + MASSMARGIN) {
1442  mr = pow2(mf / m3);
1443  betaf = sqrtpos(1. - 4. * mr);
1444  psvec = betaf * (1. + 2. * mr);
1445  psaxi = pow3(betaf);
1446 
1447  // First Z0: combine phase space with couplings.
1448  ef2 = couplingsPtr->ef2(idAbs) * psvec;
1449  efvf = couplingsPtr->efvf(idAbs) * psvec;
1450  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1451  + couplingsPtr->af2(idAbs) * psaxi;
1452  colf = (idAbs < 6) ? colQ3 : 1.;
1453 
1454  // First Z0: store sum of combinations for open outstate channels.
1455  if (onMode == 1 || onMode == 2) {
1456  gamSum3 += colf * ef2;
1457  intSum3 += colf * efvf;
1458  resSum3 += colf * vf2af2;
1459  }
1460  }
1461 
1462  // Second Z0: check that above threshold. Phase space.
1463  if (m4 > 2. * mf + MASSMARGIN) {
1464  mr = pow2(mf / m4);
1465  betaf = sqrtpos(1. - 4. * mr);
1466  psvec = betaf * (1. + 2. * mr);
1467  psaxi = pow3(betaf);
1468 
1469  // Second Z0: combine phase space with couplings.
1470  ef2 = couplingsPtr->ef2(idAbs) * psvec;
1471  efvf = couplingsPtr->efvf(idAbs) * psvec;
1472  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1473  + couplingsPtr->af2(idAbs) * psaxi;
1474  colf = (idAbs < 6) ? colQ4 : 1.;
1475 
1476  // Second Z0: store sum of combinations for open outstate channels.
1477  if (onMode == 1 || onMode == 2) {
1478  gamSum4 += colf * ef2;
1479  intSum4 += colf * efvf;
1480  resSum4 += colf * vf2af2;
1481  }
1482  }
1483 
1484  // End loop over fermions.
1485  }
1486  }
1487 
1488  // First Z0: calculate prefactors for gamma/interference/Z0 terms.
1489  gamProp3 = 4. * alpEM3 / (3. * M_PI * s3);
1490  intProp3 = gamProp3 * 2. * thetaWRat * s3 * (s3 - m2Res)
1491  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1492  resProp3 = gamProp3 * pow2(thetaWRat * s3)
1493  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1494 
1495  // First Z0: optionally only keep gamma* or Z0 term.
1496  if (gmZmode == 1) {intProp3 = 0.; resProp3 = 0.;}
1497  if (gmZmode == 2) {gamProp3 = 0.; intProp3 = 0.;}
1498 
1499  // Second Z0: calculate prefactors for gamma/interference/Z0 terms.
1500  gamProp4 = 4. * alpEM4 / (3. * M_PI * s4);
1501  intProp4 = gamProp4 * 2. * thetaWRat * s4 * (s4 - m2Res)
1502  / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1503  resProp4 = gamProp4 * pow2(thetaWRat * s4)
1504  / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1505 
1506  // Second Z0: optionally only keep gamma* or Z0 term.
1507  if (gmZmode == 1) {intProp4 = 0.; resProp4 = 0.;}
1508  if (gmZmode == 2) {gamProp4 = 0.; intProp4 = 0.;}
1509 
1510 }
1511 
1512 //--------------------------------------------------------------------------
1513 
1514 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1515 
1516 double Sigma2ffbar2gmZgmZ::sigmaHat() {
1517 
1518  // Charge/2, left- and righthanded couplings for in-fermion.
1519  int idAbs = abs(id1);
1520  double ei = 0.5 * couplingsPtr->ef(idAbs);
1521  double li = couplingsPtr->lf(idAbs);
1522  double ri = couplingsPtr->rf(idAbs);
1523 
1524  // Combine left/right gamma, interference and Z0 parts for each Z0.
1525  double left3 = ei * ei * gamProp3 * gamSum3
1526  + ei * li * intProp3 * intSum3
1527  + li * li * resProp3 * resSum3;
1528  double right3 = ei * ei * gamProp3 * gamSum3
1529  + ei * ri * intProp3 * intSum3
1530  + ri * ri * resProp3 * resSum3;
1531  double left4 = ei * ei * gamProp4 * gamSum4
1532  + ei * li * intProp4 * intSum4
1533  + li * li * resProp4 * resSum4;
1534  double right4 = ei * ei * gamProp4 * gamSum4
1535  + ei * ri * intProp4 * intSum4
1536  + ri * ri * resProp4 * resSum4;
1537 
1538  // Combine left- and right-handed couplings for the two Z0's.
1539  double sigma = sigma0 * (left3 * left4 + right3 * right4);
1540 
1541  // Correct for the running-width Z0 propagators weight in PhaseSpace.
1542  sigma /= (runBW3 * runBW4);
1543 
1544  // Initial-state colour factor. Answer.
1545  if (idAbs < 9) sigma /= 3.;
1546  return sigma;
1547 
1548 }
1549 
1550 //--------------------------------------------------------------------------
1551 
1552 // Select identity, colour and anticolour.
1553 
1554 void Sigma2ffbar2gmZgmZ::setIdColAcol() {
1555 
1556  // Flavours trivial.
1557  setId( id1, id2, 23, 23);
1558 
1559  // Colour flow topologies. Swap when antiquarks.
1560  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1561  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1562  if (id1 < 0) swapColAcol();
1563 
1564 }
1565 
1566 //--------------------------------------------------------------------------
1567 
1568 // Evaluate correlated decay flavours of the two gamma*/Z0.
1569 // Unique complication, caused by gamma*/Z0 mix different left/right.
1570 
1571 double Sigma2ffbar2gmZgmZ::weightDecayFlav( Event& process) {
1572 
1573  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1574  i1 = (process[3].id() < 0) ? 3 : 4;
1575  i2 = 7 - i1;
1576  i3 = (process[7].id() > 0) ? 7 : 8;
1577  i4 = 15 - i3;
1578  i5 = (process[9].id() > 0) ? 9 : 10;
1579  i6 = 19 - i5;
1580 
1581  // Charge/2, left- and righthanded couplings for in- and out-fermions.
1582  int idAbs = process[i1].idAbs();
1583  double ei = 0.5 * couplingsPtr->ef(idAbs);
1584  double li = couplingsPtr->lf(idAbs);
1585  double ri = couplingsPtr->rf(idAbs);
1586  idAbs = process[i3].idAbs();
1587  double e3 = 0.5 * couplingsPtr->ef(idAbs);
1588  double l3 = couplingsPtr->lf(idAbs);
1589  double r3 = couplingsPtr->rf(idAbs);
1590  idAbs = process[i5].idAbs();
1591  double e4 = 0.5 * couplingsPtr->ef(idAbs);
1592  double l4 = couplingsPtr->lf(idAbs);
1593  double r4 = couplingsPtr->rf(idAbs);
1594 
1595  // Left- and righthanded couplings combined with propagators.
1596  c3LL = ei * ei * gamProp3 * e3 * e3
1597  + ei * li * intProp3 * e3 * l3
1598  + li * li * resProp3 * l3 * l3;
1599  c3LR = ei * ei * gamProp3 * e3 * e3
1600  + ei * li * intProp3 * e3 * r3
1601  + li * li * resProp3 * r3 * r3;
1602  c3RL = ei * ei * gamProp3 * e3 * e3
1603  + ei * ri * intProp3 * e3 * l3
1604  + ri * ri * resProp3 * l3 * l3;
1605  c3RR = ei * ei * gamProp3 * e3 * e3
1606  + ei * ri * intProp3 * e3 * r3
1607  + ri * ri * resProp3 * r3 * r3;
1608  c4LL = ei * ei * gamProp4 * e4 * e4
1609  + ei * li * intProp4 * e4 * l4
1610  + li * li * resProp4 * l4 * l4;
1611  c4LR = ei * ei * gamProp4 * e4 * e4
1612  + ei * li * intProp4 * e4 * r4
1613  + li * li * resProp4 * r4 * r4;
1614  c4RL = ei * ei * gamProp4 * e4 * e4
1615  + ei * ri * intProp4 * e4 * l4
1616  + ri * ri * resProp4 * l4 * l4;
1617  c4RR = ei * ei * gamProp4 * e4 * e4
1618  + ei * ri * intProp4 * e4 * r4
1619  + ri * ri * resProp4 * r4 * r4;
1620 
1621  // Flavour weight and maximum.
1622  flavWt = (c3LL + c3LR) * (c4LL + c4LR) + (c3RL + c3RR) * (c4RL + c4RR);
1623  double flavWtMax = (c3LL + c3LR + c3RL + c3RR) * (c4LL + c4LR + c4RL + c4RR);
1624 
1625  // Done.
1626  return flavWt / flavWtMax;
1627 
1628 }
1629 
1630 //--------------------------------------------------------------------------
1631 
1632 // Evaluate weight for decay angles of the two gamma*/Z0.
1633 
1634 double Sigma2ffbar2gmZgmZ::weightDecay( Event& process, int iResBeg,
1635  int iResEnd) {
1636 
1637  // Two resonance decays, but with common weight.
1638  if (iResBeg != 5 || iResEnd != 6) return 1.;
1639 
1640  // Set up four-products and internal products.
1641  setupProd( process, i1, i2, i3, i4, i5, i6);
1642 
1643  // Flip tHat and uHat if first incoming is fermion.
1644  double tHres = tH;
1645  double uHres = uH;
1646  if (process[3].id() > 0) swap( tHres, uHres);
1647 
1648  // Kinematics factors (norm(x) = |x|^2).
1649  double fGK135 = norm( fGK( 1, 2, 3, 4, 5, 6) / tHres
1650  + fGK( 1, 2, 5, 6, 3, 4) / uHres );
1651  double fGK145 = norm( fGK( 1, 2, 4, 3, 5, 6) / tHres
1652  + fGK( 1, 2, 5, 6, 4, 3) / uHres );
1653  double fGK136 = norm( fGK( 1, 2, 3, 4, 6, 5) / tHres
1654  + fGK( 1, 2, 6, 5, 3, 4) / uHres );
1655  double fGK146 = norm( fGK( 1, 2, 4, 3, 6, 5) / tHres
1656  + fGK( 1, 2, 6, 5, 4, 3) / uHres );
1657  double fGK253 = norm( fGK( 2, 1, 5, 6, 3, 4) / tHres
1658  + fGK( 2, 1, 3, 4, 5, 6) / uHres );
1659  double fGK263 = norm( fGK( 2, 1, 6, 5, 3, 4) / tHres
1660  + fGK( 2, 1, 3, 4, 6, 5) / uHres );
1661  double fGK254 = norm( fGK( 2, 1, 5, 6, 4, 3) / tHres
1662  + fGK( 2, 1, 4, 3, 5, 6) / uHres );
1663  double fGK264 = norm( fGK( 2, 1, 6, 5, 4, 3) / tHres
1664  + fGK( 2, 1, 4, 3, 6, 5) / uHres );
1665 
1666  // Weight and maximum.
1667  double wt = c3LL * c4LL * fGK135 + c3LR * c4LL * fGK145
1668  + c3LL * c4LR * fGK136 + c3LR * c4LR * fGK146
1669  + c3RL * c4RL * fGK253 + c3RR * c4RL * fGK263
1670  + c3RL * c4RR * fGK254 + c3RR * c4RR * fGK264;
1671  double wtMax = 16. * s3 * s4 * flavWt
1672  * ( (tHres*tHres + uHres*uHres + 2. * sH * (s3 + s4)) / (tHres * uHres)
1673  - s3 * s4 * (1. / (tHres*tHres) + 1. / (uHres*uHres)) );
1674 
1675  // Done.
1676  return wt / wtMax;
1677 
1678 }
1679 
1680 //==========================================================================
1681 
1682 // Sigma2ffbar2ZW class.
1683 // Cross section for f fbar' -> Z0 W+- (f is quark or lepton).
1684 
1685 //--------------------------------------------------------------------------
1686 
1687 // Initialize process.
1688 
1689 void Sigma2ffbar2ZW::initProc() {
1690 
1691  // Store W+- mass and width for propagator.
1692  mW = particleDataPtr->m0(24);
1693  widW = particleDataPtr->mWidth(24);
1694  mWS = mW*mW;
1695  mwWS = pow2(mW * widW);
1696 
1697  // Left-handed couplings for up/nu- and down/e-type quarks.
1698  lun = (hasLeptonBeams) ? couplingsPtr->lf(12) : couplingsPtr->lf(2);
1699  lde = (hasLeptonBeams) ? couplingsPtr->lf(11) : couplingsPtr->lf(1);
1700 
1701  // Common weak coupling factor.
1702  sin2thetaW = couplingsPtr->sin2thetaW();
1703  cos2thetaW = couplingsPtr->cos2thetaW();
1704  thetaWRat = 1. / (4. * cos2thetaW);
1705  cotT = sqrt(cos2thetaW / sin2thetaW);
1706  thetaWpt = (9. - 8. * sin2thetaW) / 4.;
1707  thetaWmm = (8. * sin2thetaW - 6.) / 4.;
1708 
1709  // Secondary open width fractions.
1710  openFracPos = particleDataPtr->resOpenFrac(23, 24);
1711  openFracNeg = particleDataPtr->resOpenFrac(23, -24);
1712 
1713 }
1714 
1715 //--------------------------------------------------------------------------
1716 
1717 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1718 
1719 void Sigma2ffbar2ZW::sigmaKin() {
1720 
1721  // Evaluate cross section, as programmed by Merlin Kole (after tidying),
1722  // based on Brown, Sahdev, Mikaelian, Phys Rev. D20 (1979) 1069.
1723  /*
1724  double resBW = 1. / (pow2(sH - mWS) + mwWS);
1725  double prefac = 12.0 * M_PI * pow2(alpEM) / (sH2 * 8. * sin2thetaW);
1726  double temp1 = tH * uH - s3 * s4;
1727  double temp2 = temp1 / (s3 * s4);
1728  double temp3 = (s3 + s4) / sH;
1729  double temp4 = s3 * s4 / sH;
1730  double partA = temp2 * (0.25 - 0.5 * temp3
1731  + (pow2(s3 + s4) + 8. * s3 * s4)/(4. * sH2) )
1732  + (s3 + s4)/(s3 * s4) * (sH/2. - s3 - s4 + pow2(s3 - s4)/(2. * sH));
1733  double partB1 = lun * (0.25 * temp2 * (1. - temp3 - 4. * temp4 / tH)
1734  + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / tH) );
1735  double partB2 = lde * ( 0.25 * temp2 * (1.- temp3 - 4. * temp4 / uH)
1736  + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / uH) );
1737  double partE = 0.25 * temp2 + 0.5 *(s3 + s4) / temp4;
1738  sigma0 = prefac * (pow2(cotT) * sH2 * resBW * partA
1739  + 2.* sH * cotT * resBW * (sH - mWS) * (partB2 - partB1)
1740  + pow2(lun - lde) * partE + pow2(lde) * temp1/uH2
1741  + pow2(lun) * temp1/tH2 + 2. * lun * lde * sH * (s3 + s4) / (uH * tH));
1742  */
1743 
1744  // Evaluate cross section. Expression from EHLQ, with bug fix,
1745  // but can still give negative cross section so suspect.
1746  double resBW = 1. / (pow2(sH - mWS) + mwWS);
1747  sigma0 = (M_PI / sH2) * 0.5 * pow2(alpEM / sin2thetaW);
1748  sigma0 *= sH * resBW * (thetaWpt * pT2 + thetaWmm * (s3 + s4))
1749  + (sH - mWS) * resBW * sH * (pT2 - s3 - s4) * (lun / tH - lde / uH)
1750  + thetaWRat * sH * pT2 * ( lun*lun / tH2 + lde*lde / uH2 )
1751  + 2. * thetaWRat * sH * (s3 + s4) * lun * lde / (tH * uH);
1752  // Need to protect against negative cross sections at times.
1753  sigma0 = max(0., sigma0);
1754 
1755 }
1756 
1757 //--------------------------------------------------------------------------
1758 
1759 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1760 
1761 double Sigma2ffbar2ZW::sigmaHat() {
1762 
1763  // CKM and colour factors.
1764  double sigma = sigma0;
1765  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1766 
1767  // Corrections for secondary widths in Z0 and W+- decays.
1768  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
1769  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
1770 
1771  // Answer.
1772  return sigma;
1773 
1774 }
1775 
1776 //--------------------------------------------------------------------------
1777 
1778 // Select identity, colour and anticolour.
1779 
1780 void Sigma2ffbar2ZW::setIdColAcol() {
1781 
1782  // Sign of outgoing W.
1783  int sign = 1 - 2 * (abs(id1)%2);
1784  if (id1 < 0) sign = -sign;
1785  setId( id1, id2, 23, 24 * sign);
1786 
1787  // tHat is defined between (f, W-) or (fbar, W+),
1788  // so OK for u/ubar on side 1, but must swap tHat <-> uHat if d/dbar.
1789  if (abs(id1)%2 == 1) swapTU = true;
1790 
1791  // Colour flow topologies. Swap when antiquarks.
1792  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1793  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1794  if (id1 < 0) swapColAcol();
1795 
1796 }
1797 
1798 //--------------------------------------------------------------------------
1799 
1800 // Evaluate weight for Z0 and W+- decay angles.
1801 
1802 double Sigma2ffbar2ZW::weightDecay( Event& process, int iResBeg, int iResEnd) {
1803 
1804  // Two resonance decays, but with common weight.
1805  if (iResBeg != 5 || iResEnd != 6) return 1.;
1806 
1807  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6)
1808  // with f' fbar' from W+- and f" fbar" from Z0 (note flip Z0 <-> W+-).
1809  int i1 = (process[3].id() < 0) ? 3 : 4;
1810  int i2 = 7 - i1;
1811  int i3 = (process[9].id() > 0) ? 9 : 10;
1812  int i4 = 19 - i3;
1813  int i5 = (process[7].id() > 0) ? 7 : 8;
1814  int i6 = 15 - i5;
1815 
1816  // Set up four-products and internal products.
1817  setupProd( process, i1, i2, i3, i4, i5, i6);
1818 
1819  // Swap tHat and uHat if incoming fermion is downtype.
1820  double tHres = tH;
1821  double uHres = uH;
1822  if (process[i2].id()%2 == 1) swap( tHres, uHres);
1823 
1824  // Couplings of incoming (anti)fermions and outgoing from Z0.
1825  int idAbs = process[i1].idAbs();
1826  double ai = couplingsPtr->af(idAbs);
1827  double li1 = couplingsPtr->lf(idAbs);
1828  idAbs = process[i2].idAbs();
1829  double li2 = couplingsPtr->lf(idAbs);
1830  idAbs = process[i5].idAbs();
1831  double l4 = couplingsPtr->lf(idAbs);
1832  double r4 = couplingsPtr->rf(idAbs);
1833 
1834  // W propagator/interference factor.
1835  double Wint = cos2thetaW * (sH - mWS) / (pow2(sH - mWS) + mwWS);
1836 
1837  // Combinations of couplings and kinematics (norm(x) = |x|^2).
1838  double aWZ = li2 / tHres - 2. * Wint * ai;
1839  double bWZ = li1 / uHres + 2. * Wint * ai;
1840  double fGK135 = norm( aWZ * fGK( 1, 2, 3, 4, 5, 6)
1841  + bWZ * fGK( 1, 2, 5, 6, 3, 4) );
1842  double fGK136 = norm( aWZ * fGK( 1, 2, 3, 4, 6, 5)
1843  + bWZ * fGK( 1, 2, 6, 5, 3, 4) );
1844  double xiT = xiGK( tHres, uHres);
1845  double xiU = xiGK( uHres, tHres);
1846  double xjTU = xjGK( tHres, uHres);
1847 
1848  // Weight and maximum weight.
1849  double wt = l4*l4 * fGK135 + r4*r4 * fGK136;
1850  double wtMax = 4. * s3 * s4 * (l4*l4 + r4*r4)
1851  * (aWZ * aWZ * xiT + bWZ * bWZ * xiU + aWZ * bWZ * xjTU);
1852 
1853  // Done.
1854  return wt / wtMax;
1855 
1856 }
1857 
1858 //==========================================================================
1859 
1860 // Sigma2ffbar2WW class.
1861 // Cross section for f fbar -> W- W+ (f is quark or lepton).
1862 
1863 //--------------------------------------------------------------------------
1864 
1865 // Initialize process.
1866 
1867 void Sigma2ffbar2WW::initProc() {
1868 
1869  // Store Z0 mass and width for propagator. Common coupling factor.
1870  mZ = particleDataPtr->m0(23);
1871  widZ = particleDataPtr->mWidth(23);
1872  mZS = mZ*mZ;
1873  mwZS = pow2(mZ * widZ);
1874  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
1875 
1876  // Secondary open width fraction.
1877  openFracPair = particleDataPtr->resOpenFrac(24, -24);
1878 
1879 }
1880 
1881 //--------------------------------------------------------------------------
1882 
1883 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1884 
1885 void Sigma2ffbar2WW::sigmaKin() {
1886 
1887  // Cross section part common for all incoming flavours.
1888  sigma0 = (M_PI / sH2) * pow2(alpEM);
1889 
1890  // Z0 propagator and gamma*/Z0 interference.
1891  double Zprop = sH2 / (pow2(sH - mZS) + mwZS);
1892  double Zint = Zprop * (1. - mZS / sH);
1893 
1894  // Common coupling factors (g = gamma*, Z = Z0, f = t-channel fermion).
1895  cgg = 0.5;
1896  cgZ = thetaWRat * Zint;
1897  cZZ = 0.5 * pow2(thetaWRat) * Zprop;
1898  cfg = thetaWRat;
1899  cfZ = pow2(thetaWRat) * Zint;
1900  cff = pow2(thetaWRat);
1901 
1902  // Kinematical functions.
1903  double rat34 = sH * (2. * (s3 + s4) + pT2) / (s3 * s4);
1904  double lambdaS = pow2(sH - s3 - s4) - 4. * s3 * s4;
1905  double intA = (sH - s3 - s4) * rat34 / sH;
1906  double intB = 4. * (s3 + s4 - pT2);
1907  gSS = (lambdaS * rat34 + 12. * sH * pT2) / sH2;
1908  gTT = rat34 + 4. * sH * pT2 / tH2;
1909  gST = intA + intB / tH;
1910  gUU = rat34 + 4. * sH * pT2 / uH2;
1911  gSU = intA + intB / uH;
1912 
1913 }
1914 
1915 //--------------------------------------------------------------------------
1916 
1917 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1918 
1919 double Sigma2ffbar2WW::sigmaHat() {
1920 
1921  // Flavour-specific couplings.
1922  int idAbs = abs(id1);
1923  double ei = couplingsPtr->ef(idAbs);
1924  double vi = couplingsPtr->vf(idAbs);
1925  double ai = couplingsPtr->af(idAbs);
1926 
1927  // Combine, with different cases for up- and down-type in-flavours.
1928  double sigma = sigma0;
1929  sigma *= (idAbs%2 == 1)
1930  ? (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
1931  + (cfg * ei + cfZ * (vi + ai)) * gST + cff * gTT
1932  : (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
1933  - (cfg * ei + cfZ * (vi + ai)) * gSU + cff * gUU;
1934 
1935  // Initial-state colour factor. Correction for secondary widths. Answer.
1936  if (idAbs < 9) sigma /= 3.;
1937  sigma *= openFracPair;
1938  return sigma;
1939 
1940 }
1941 
1942 //--------------------------------------------------------------------------
1943 
1944 // Select identity, colour and anticolour.
1945 
1946 void Sigma2ffbar2WW::setIdColAcol() {
1947 
1948  // Always order W- W+, i.e. W- first.
1949  setId( id1, id2, -24, 24);
1950 
1951  // tHat is defined between (f, W-) or (fbar, W+),
1952  if (id1 < 0) swapTU = true;
1953 
1954  // Colour flow topologies. Swap when antiquarks.
1955  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1956  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1957  if (id1 < 0) swapColAcol();
1958 
1959 }
1960 
1961 //--------------------------------------------------------------------------
1962 
1963 // Evaluate weight for W+ and W- decay angles.
1964 
1965 double Sigma2ffbar2WW::weightDecay( Event& process, int iResBeg, int iResEnd) {
1966 
1967  // Two resonance decays, but with common weight.
1968  if (iResBeg != 5 || iResEnd != 6) return 1.;
1969 
1970  // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1971  // with f' fbar' from W- and f" fbar" from W+.
1972  int i1 = (process[3].id() < 0) ? 3 : 4;
1973  int i2 = 7 - i1;
1974  int i3 = (process[7].id() > 0) ? 7 : 8;
1975  int i4 = 15 - i3;
1976  int i5 = (process[9].id() > 0) ? 9 : 10;
1977  int i6 = 19 - i5;
1978 
1979  // Set up four-products and internal products.
1980  setupProd( process, i1, i2, i3, i4, i5, i6);
1981 
1982  // tHat and uHat of fbar f -> W- W+ opposite to previous convention.
1983  double tHres = uH;
1984  double uHres = tH;
1985 
1986  // Couplings of incoming (anti)fermion.
1987  int idAbs = process[i1].idAbs();
1988  double ai = couplingsPtr->af(idAbs);
1989  double li = couplingsPtr->lf(idAbs);
1990  double ri = couplingsPtr->rf(idAbs);
1991 
1992  // gamma*/Z0 propagator/interference factor.
1993  double Zint = mZS * (sH - mZS) / (pow2(sH - mZS) + mwZS);
1994 
1995  // Combinations of couplings and kinematics (norm(x) = |x|^2).
1996  double dWW = (li * Zint + ai) / sH;
1997  double aWW = dWW + 0.5 * (ai + 1.) / tHres;
1998  double bWW = dWW + 0.5 * (ai - 1.) / uHres;
1999  double cWW = ri * Zint / sH;
2000  double fGK135 = norm( aWW * fGK( 1, 2, 3, 4, 5, 6)
2001  - bWW * fGK( 1, 2, 5, 6, 3, 4) );
2002  double fGK253 = norm( cWW * ( fGK( 2, 1, 5, 6, 3, 4)
2003  - fGK( 2, 1, 3, 4, 5, 6) ) );
2004  double xiT = xiGK( tHres, uHres);
2005  double xiU = xiGK( uHres, tHres);
2006  double xjTU = xjGK( tHres, uHres);
2007 
2008  // Weight and maximum weight.
2009  double wt = fGK135 + fGK253;
2010  double wtMax = 4. * s3 * s4
2011  * ( aWW * aWW * xiT + bWW * bWW * xiU - aWW * bWW * xjTU
2012  + cWW * cWW * (xiT + xiU - xjTU) );
2013 
2014  // Done.
2015  return wt / wtMax;
2016 }
2017 
2018 //==========================================================================
2019 
2020 // Sigma2ffbargmZggm class.
2021 // Collects common methods for f fbar -> gamma*/Z0 g/gamma and permutations.
2022 
2023 //--------------------------------------------------------------------------
2024 
2025 // Initialize process.
2026 
2027 void Sigma2ffbargmZggm::initProc() {
2028 
2029  // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
2030  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
2031 
2032  // Store Z0 mass and width for propagator.
2033  mRes = particleDataPtr->m0(23);
2034  GammaRes = particleDataPtr->mWidth(23);
2035  m2Res = mRes*mRes;
2036  GamMRat = GammaRes / mRes;
2037  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
2038  * couplingsPtr->cos2thetaW());
2039 
2040  // Set pointer to particle properties and decay table.
2041  particlePtr = particleDataPtr->particleDataEntryPtr(23);
2042 
2043 }
2044 
2045 //--------------------------------------------------------------------------
2046 
2047 // Evaluate sum of flavour couplings times phase space.
2048 
2049 void Sigma2ffbargmZggm::flavSum() {
2050 
2051  // Coupling factors for Z0 subsystem.
2052  double alpSZ = couplingsPtr->alphaS(s3);
2053  double colQZ = 3. * (1. + alpSZ / M_PI);
2054 
2055  // Reset quantities to sum. Declare variables in loop.
2056  gamSum = 0.;
2057  intSum = 0.;
2058  resSum = 0.;
2059  int onMode;
2060  double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
2061 
2062  // Loop over all Z0 decay channels.
2063  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
2064  int idAbs = abs( particlePtr->channel(i).product(0) );
2065 
2066  // Only contributions from three fermion generations, except top.
2067  if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
2068  mf = particleDataPtr->m0(idAbs);
2069 
2070  // Check that above threshold. Phase space.
2071  if (m3 > 2. * mf + MASSMARGIN) {
2072  mr = pow2(mf / m3);
2073  betaf = sqrtpos(1. - 4. * mr);
2074  psvec = betaf * (1. + 2. * mr);
2075  psaxi = pow3(betaf);
2076 
2077  // Combine phase space with couplings.
2078  ef2 = couplingsPtr->ef2(idAbs) * psvec;
2079  efvf = couplingsPtr->efvf(idAbs) * psvec;
2080  vf2af2 = couplingsPtr->vf2(idAbs) * psvec
2081  + couplingsPtr->af2(idAbs) * psaxi;
2082  colf = (idAbs < 6) ? colQZ : 1.;
2083 
2084  // Store sum of combinations. For outstate only open channels.
2085  onMode = particlePtr->channel(i).onMode();
2086  if (onMode == 1 || onMode == 2) {
2087  gamSum += colf * ef2;
2088  intSum += colf * efvf;
2089  resSum += colf * vf2af2;
2090  }
2091 
2092  // End loop over fermions.
2093  }
2094  }
2095  }
2096 
2097  // Done. Return values in gamSum, intSum and resSum.
2098 
2099 }
2100 
2101 //--------------------------------------------------------------------------
2102 
2103 // Calculate common parts of gamma/interference/Z0 propagator terms.
2104 
2105 void Sigma2ffbargmZggm::propTerm() {
2106 
2107  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2108  gamProp = 4. * alpEM / (3. * M_PI * s3);
2109  intProp = gamProp * 2. * thetaWRat * s3 * (s3 - m2Res)
2110  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2111  resProp = gamProp * pow2(thetaWRat * s3)
2112  / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2113 
2114  // Optionally only keep gamma* or Z0 term.
2115  if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
2116  if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
2117 
2118 }
2119 
2120 //--------------------------------------------------------------------------
2121 
2122 // Evaluate weight for gamma*/Z0 decay angle.
2123 
2124 double Sigma2ffbargmZggm::weightDecay( Event& process, int iResBeg,
2125  int iResEnd) {
2126 
2127  // Z should sit in entry 5 and one more parton in entry 6.
2128  if (iResBeg != 5 || iResEnd != 6) return 1.;
2129 
2130  // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2131  // where f' fbar' come from gamma*/Z0 decay.
2132  int i1, i2;
2133  int i3 = (process[7].id() > 0) ? 7 : 8;
2134  int i4 = 15 - i3;
2135 
2136  // Order so that fbar(1) f(2) -> gamma*/Z0 g/gamma.
2137  if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2138  i1 = (process[3].id() < 0) ? 3 : 4;
2139  i2 = 7 - i1;
2140 
2141  // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) gamma*/Z0.
2142  } else if (process[3].idAbs() < 20) {
2143  i1 = (process[3].id() < 0) ? 3 : 6;
2144  i2 = 9 - i1;
2145  } else {
2146  i1 = (process[4].id() < 0) ? 4 : 6;
2147  i2 = 10 - i1;
2148  }
2149 
2150  // Charge/2, left- and righthanded couplings for in- and out-fermion.
2151  int id1Abs = process[i1].idAbs();
2152  double ei = 0.5 * couplingsPtr->ef(id1Abs);
2153  double li = couplingsPtr->lf(id1Abs);
2154  double ri = couplingsPtr->rf(id1Abs);
2155  int id3Abs = process[i3].idAbs();
2156  double ef = 0.5 * couplingsPtr->ef(id3Abs);
2157  double lf = couplingsPtr->lf(id3Abs);
2158  double rf = couplingsPtr->rf(id3Abs);
2159 
2160  // Combinations of left/right for in/out, gamma*/interference/Z0.
2161  double clilf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*lf
2162  + li*li * resProp * lf*lf;
2163  double clirf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*rf
2164  + li*li * resProp * rf*rf;
2165  double crilf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*lf
2166  + ri*ri * resProp * lf*lf;
2167  double crirf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*rf
2168  + ri*ri * resProp * rf*rf;
2169 
2170  // Evaluate four-vector products.
2171  double p13 = process[i1].p() * process[i3].p();
2172  double p14 = process[i1].p() * process[i4].p();
2173  double p23 = process[i2].p() * process[i3].p();
2174  double p24 = process[i2].p() * process[i4].p();
2175 
2176  // Calculate weight and its maximum.
2177  double wt = (clilf + crirf) * (p13*p13 + p24*p24)
2178  + (clirf + crilf) * (p14*p14 + p23*p23) ;
2179  double wtMax = (clilf + clirf + crilf + crirf)
2180  * (pow2(p13 + p14) + pow2(p23 + p24));
2181 
2182  // Done.
2183  return (wt / wtMax);
2184 
2185 }
2186 
2187 //==========================================================================
2188 
2189 // Sigma2qqbar2gmZg class.
2190 // Cross section for q qbar -> gamma*/Z0 g.
2191 
2192 //--------------------------------------------------------------------------
2193 
2194 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2195 
2196 void Sigma2qqbar2gmZg::sigmaKin() {
2197 
2198  // Cross section part common for all incoming flavours.
2199  sigma0 = (M_PI / sH2) * (alpEM * alpS)
2200  * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2201 
2202  // Calculate flavour sums for final state.
2203  flavSum();
2204 
2205  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2206  propTerm();
2207 
2208 }
2209 
2210 //--------------------------------------------------------------------------
2211 
2212 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2213 
2214 double Sigma2qqbar2gmZg::sigmaHat() {
2215 
2216  // Combine gamma, interference and Z0 parts.
2217  int idAbs = abs(id1);
2218  double sigma = sigma0
2219  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2220  + couplingsPtr->efvf(idAbs) * intProp * intSum
2221  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2222 
2223  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2224  sigma /= runBW3;
2225 
2226  // Answer.
2227  return sigma;
2228 
2229 }
2230 
2231 //--------------------------------------------------------------------------
2232 
2233 // Select identity, colour and anticolour.
2234 
2235 void Sigma2qqbar2gmZg::setIdColAcol() {
2236 
2237  // Flavours trivial.
2238  setId( id1, id2, 23, 21);
2239 
2240  // Colour flow topologies. Swap when antiquarks.
2241  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2242  if (id1 < 0) swapColAcol();
2243 
2244 }
2245 
2246 //==========================================================================
2247 
2248 // Sigma2qg2gmZq class.
2249 // Cross section for q g -> gamma*/Z0 q.
2250 
2251 //--------------------------------------------------------------------------
2252 
2253 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2254 
2255 void Sigma2qg2gmZq::sigmaKin() {
2256 
2257  // Cross section part common for all incoming flavours.
2258  sigma0 = (M_PI / sH2) * (alpEM * alpS)
2259  * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2260 
2261  // Calculate flavour sums for final state.
2262  flavSum();
2263 
2264  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2265  propTerm();
2266 
2267 }
2268 
2269 //--------------------------------------------------------------------------
2270 
2271 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2272 
2273 double Sigma2qg2gmZq::sigmaHat() {
2274 
2275  // Combine gamma, interference and Z0 parts.
2276  int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2277  double sigma = sigma0
2278  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2279  + couplingsPtr->efvf(idAbs) * intProp * intSum
2280  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2281 
2282  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2283  sigma /= runBW3;
2284 
2285  // Answer.
2286  return sigma;
2287 
2288 }
2289 
2290 //--------------------------------------------------------------------------
2291 
2292 // Select identity, colour and anticolour.
2293 
2294 void Sigma2qg2gmZq::setIdColAcol() {
2295 
2296  // Flavour set up for q g -> gamma*/Z0 q.
2297  int idq = (id2 == 21) ? id1 : id2;
2298  setId( id1, id2, 23, idq);
2299 
2300  // tH defined between f and f': must swap tHat <-> uHat if q g in.
2301  swapTU = (id2 == 21);
2302 
2303  // Colour flow topologies. Swap when antiquarks.
2304  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2305  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2306  if (idq < 0) swapColAcol();
2307 
2308 }
2309 
2310 //==========================================================================
2311 
2312 // Sigma2ffbar2gmZgm class.
2313 // Cross section for f fbar -> gamma*/Z0 gamma.
2314 
2315 //--------------------------------------------------------------------------
2316 
2317 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2318 
2319 void Sigma2ffbar2gmZgm::sigmaKin() {
2320 
2321  // Cross section part common for all incoming flavours.
2322  sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2323  * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2324 
2325  // Calculate flavour sums for final state.
2326  flavSum();
2327 
2328  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2329  propTerm();
2330 
2331 
2332 }
2333 
2334 //--------------------------------------------------------------------------
2335 
2336 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2337 
2338 double Sigma2ffbar2gmZgm::sigmaHat() {
2339 
2340  // Combine gamma, interference and Z0 parts.
2341  int idAbs = abs(id1);
2342  double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2343  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2344  + couplingsPtr->efvf(idAbs) * intProp * intSum
2345  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2346 
2347  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2348  sigma /= runBW3;
2349 
2350  // Colour factor. Answer.
2351  if (idAbs < 9) sigma /= 3.;
2352  return sigma;
2353 
2354 }
2355 
2356 //--------------------------------------------------------------------------
2357 
2358 // Select identity, colour and anticolour.
2359 
2360 void Sigma2ffbar2gmZgm::setIdColAcol() {
2361 
2362  // Flavours trivial.
2363  setId( id1, id2, 23, 22);
2364 
2365  // Colour flow topologies. Swap when antiquarks.
2366  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2367  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2368  if (id1 < 0) swapColAcol();
2369 
2370 }
2371 
2372 //==========================================================================
2373 
2374 // Sigma2fgm2gmZf class.
2375 // Cross section for f gamma -> gamma*/Z0 f'.
2376 
2377 //--------------------------------------------------------------------------
2378 
2379 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2380 
2381 void Sigma2fgm2gmZf::sigmaKin() {
2382 
2383  // Cross section part common for all incoming flavours.
2384  sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2385  * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (- sH * uH);
2386 
2387  // Calculate flavour sums for final state.
2388  flavSum();
2389 
2390  // Calculate prefactors for gamma/interference/Z0 cross section terms.
2391  propTerm();
2392 
2393 }
2394 
2395 //--------------------------------------------------------------------------
2396 
2397 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2398 
2399 double Sigma2fgm2gmZf::sigmaHat() {
2400 
2401  // Combine gamma, interference and Z0 parts.
2402  int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
2403  double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2404  * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2405  + couplingsPtr->efvf(idAbs) * intProp * intSum
2406  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2407 
2408  // Correct for the running-width Z0 propagater weight in PhaseSpace.
2409  sigma /= runBW3;
2410 
2411  // Answer.
2412  return sigma;
2413 
2414 }
2415 
2416 //--------------------------------------------------------------------------
2417 
2418 // Select identity, colour and anticolour.
2419 
2420 void Sigma2fgm2gmZf::setIdColAcol() {
2421 
2422  // Flavour set up for q gamma -> gamma*/Z0 q.
2423  int idq = (id2 == 22) ? id1 : id2;
2424  setId( id1, id2, 23, idq);
2425 
2426  // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2427  swapTU = (id2 == 22);
2428 
2429  // Colour flow topologies. Swap when antiquarks.
2430  if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2431  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2432  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2433  if (idq < 0) swapColAcol();
2434 
2435 }
2436 
2437 //==========================================================================
2438 
2439 // Sigma2ffbarWggm class.
2440 // Collects common methods for f fbar -> W+- g/gamma and permutations.
2441 
2442 //--------------------------------------------------------------------------
2443 
2444 // Evaluate weight for W+- decay angle.
2445 
2446 double Sigma2ffbarWggm::weightDecay( Event& process, int iResBeg,
2447  int iResEnd) {
2448 
2449  // W should sit in entry 5 and one more parton in entry 6.
2450  if (iResBeg != 5 || iResEnd != 6) return 1.;
2451 
2452  // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2453  // where f' fbar' come from W+- decay.
2454  int i1, i2;
2455  int i3 = (process[7].id() > 0) ? 7 : 8;
2456  int i4 = 15 - i3;
2457 
2458  // Order so that fbar(1) f(2) -> W+- g/gamma.
2459  if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2460  i1 = (process[3].id() < 0) ? 3 : 4;
2461  i2 = 7 - i1;
2462 
2463  // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) W+-.
2464  } else if (process[3].idAbs() < 20) {
2465  i1 = (process[3].id() < 0) ? 3 : 6;
2466  i2 = 9 - i1;
2467  } else {
2468  i1 = (process[4].id() < 0) ? 4 : 6;
2469  i2 = 10 - i1;
2470  }
2471 
2472  // Evaluate four-vector products.
2473  double p13 = process[i1].p() * process[i3].p();
2474  double p14 = process[i1].p() * process[i4].p();
2475  double p23 = process[i2].p() * process[i3].p();
2476  double p24 = process[i2].p() * process[i4].p();
2477 
2478  // Calculate weight and its maximum.
2479  double wt = pow2(p13) + pow2(p24);
2480  double wtMax = pow2(p13 + p14) + pow2(p23 + p24);
2481 
2482  // Done.
2483  return (wt / wtMax);
2484 
2485 }
2486 
2487 //==========================================================================
2488 
2489 // Sigma2qqbar2Wg class.
2490 // Cross section for q qbar' -> W+- g.
2491 
2492 //--------------------------------------------------------------------------
2493 
2494 // Initialize process.
2495 
2496 void Sigma2qqbar2Wg::initProc() {
2497 
2498  // Secondary open width fractions, relevant for top (or heavier).
2499  openFracPos = particleDataPtr->resOpenFrac(24);
2500  openFracNeg = particleDataPtr->resOpenFrac(-24);
2501 
2502 }
2503 
2504 //--------------------------------------------------------------------------
2505 
2506 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2507 
2508 void Sigma2qqbar2Wg::sigmaKin() {
2509 
2510  // Cross section part common for all incoming flavours.
2511  sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2512  * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2513 
2514 }
2515 
2516 //--------------------------------------------------------------------------
2517 
2518 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2519 
2520 double Sigma2qqbar2Wg::sigmaHat() {
2521 
2522  // CKM factor. Secondary width for W+ or W-.
2523  double sigma = sigma0 * couplingsPtr->V2CKMid(abs(id1), abs(id2));
2524  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2525  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2526 
2527  // Answer.
2528  return sigma;
2529 
2530 }
2531 
2532 //--------------------------------------------------------------------------
2533 
2534 // Select identity, colour and anticolour.
2535 
2536 void Sigma2qqbar2Wg::setIdColAcol() {
2537 
2538  // Sign of outgoing W.
2539  int sign = 1 - 2 * (abs(id1)%2);
2540  if (id1 < 0) sign = -sign;
2541  setId( id1, id2, 24 * sign, 21);
2542 
2543  // Colour flow topologies. Swap when antiquarks.
2544  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2545  if (id1 < 0) swapColAcol();
2546 
2547 }
2548 
2549 //==========================================================================
2550 
2551 // Sigma2qg2Wq class.
2552 // Cross section for q g -> W+- q'.
2553 
2554 //--------------------------------------------------------------------------
2555 
2556 // Initialize process.
2557 
2558 void Sigma2qg2Wq::initProc() {
2559 
2560  // Secondary open width fractions, relevant for top (or heavier).
2561  openFracPos = particleDataPtr->resOpenFrac(24);
2562  openFracNeg = particleDataPtr->resOpenFrac(-24);
2563 
2564 }
2565 
2566 //--------------------------------------------------------------------------
2567 
2568 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2569 
2570 void Sigma2qg2Wq::sigmaKin() {
2571 
2572  // Cross section part common for all incoming flavours.
2573  sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2574  * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2575 
2576 }
2577 
2578 //--------------------------------------------------------------------------
2579 
2580 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2581 
2582 double Sigma2qg2Wq::sigmaHat() {
2583 
2584  // CKM factor. Secondary width for W+ or W-.
2585  int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2586  double sigma = sigma0 * couplingsPtr->V2CKMsum(idAbs);
2587  int idUp = (id2 == 21) ? id1 : id2;
2588  if (idAbs%2 == 1) idUp = -idUp;
2589  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2590 
2591  // Answer.
2592  return sigma;
2593 
2594 }
2595 
2596 //--------------------------------------------------------------------------
2597 
2598 // Select identity, colour and anticolour.
2599 
2600 void Sigma2qg2Wq::setIdColAcol() {
2601 
2602  // Sign of outgoing W. Flavour of outgoing quark.
2603  int idq = (id2 == 21) ? id1 : id2;
2604  int sign = 1 - 2 * (abs(idq)%2);
2605  if (idq < 0) sign = -sign;
2606  id4 = couplingsPtr->V2CKMpick(idq);
2607 
2608  // Flavour set up for q g -> W q.
2609  setId( id1, id2, 24 * sign, id4);
2610 
2611  // tH defined between f and f': must swap tHat <-> uHat if q g in.
2612  swapTU = (id2 == 21);
2613 
2614  // Colour flow topologies. Swap when antiquarks.
2615  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2616  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2617  if (idq < 0) swapColAcol();
2618 
2619 }
2620 
2621 //==========================================================================
2622 
2623 // Sigma2ffbar2Wgm class.
2624 // Cross section for f fbar' -> W+- gamma.
2625 
2626 //--------------------------------------------------------------------------
2627 
2628 // Initialize process.
2629 
2630 void Sigma2ffbar2Wgm::initProc() {
2631 
2632  // Secondary open width fractions, relevant for top (or heavier).
2633  openFracPos = particleDataPtr->resOpenFrac(24);
2634  openFracNeg = particleDataPtr->resOpenFrac(-24);
2635 
2636 }
2637 
2638 //--------------------------------------------------------------------------
2639 
2640 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2641 
2642 void Sigma2ffbar2Wgm::sigmaKin() {
2643 
2644  // Cross section part common for all incoming flavours.
2645  sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2646  * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2647 }
2648 
2649 //--------------------------------------------------------------------------
2650 
2651 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2652 
2653 double Sigma2ffbar2Wgm::sigmaHat() {
2654 
2655  // Extrafactor different for e nu and q qbar' instate.
2656  int id1Abs = abs(id1);
2657  int id2Abs = abs(id2);
2658  double chgUp = (id1Abs > 10) ? 0. : 2./3.;
2659  double sigma = sigma0 * pow2( chgUp - tH / (tH + uH) );
2660 
2661  // CKM and colour factors. Secondary width for W+ or W-.
2662  if (id1Abs < 9) sigma *= couplingsPtr->V2CKMid(id1Abs, id2Abs) / 3.;
2663  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2664  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2665 
2666  // Answer.
2667  return sigma;
2668 
2669 }
2670 
2671 //--------------------------------------------------------------------------
2672 
2673 // Select identity, colour and anticolour.
2674 
2675 void Sigma2ffbar2Wgm::setIdColAcol() {
2676 
2677  // Sign of outgoing W.
2678  int sign = 1 - 2 * (abs(id1)%2);
2679  if (id1 < 0) sign = -sign;
2680  setId( id1, id2, 24 * sign, 22);
2681 
2682  // tH defined between (f,W-) or (fbar',W+).
2683  swapTU = (sign * id1 > 0);
2684 
2685  // Colour flow topologies. Swap when antiquarks.
2686  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2687  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2688  if (id1 < 0) swapColAcol();
2689 
2690 }
2691 
2692 //==========================================================================
2693 
2694 // Sigma2fgm2Wf class.
2695 // Cross section for f gamma -> W+- f'.
2696 
2697 //--------------------------------------------------------------------------
2698 
2699 // Initialize process.
2700 
2701 void Sigma2fgm2Wf::initProc() {
2702 
2703  // Secondary open width fractions, relevant for top (or heavier).
2704  openFracPos = particleDataPtr->resOpenFrac(24);
2705  openFracNeg = particleDataPtr->resOpenFrac(-24);
2706 
2707 }
2708 
2709 //--------------------------------------------------------------------------
2710 
2711 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2712 
2713 void Sigma2fgm2Wf::sigmaKin() {
2714 
2715  // Cross section part common for all incoming flavours.
2716  sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2717  * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (pT2 * s3 - sH * uH);
2718 
2719 }
2720 
2721 //--------------------------------------------------------------------------
2722 
2723 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2724 
2725 double Sigma2fgm2Wf::sigmaHat() {
2726 
2727  // Extrafactor dependent on charge of incoming fermion.
2728  int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
2729  double charge = (idAbs > 10) ? 1. : ( (idAbs%2 == 1) ? 1./3. : 2./3. );
2730  double sigma = sigma0 * pow2( charge - sH / (sH + uH) );
2731 
2732  // CKM factor. Secondary width for W+ or W-.
2733  sigma *= couplingsPtr->V2CKMsum(idAbs);
2734  int idUp = (id2 == 22) ? id1 : id2;
2735  if (idAbs%2 == 1) idUp = -idUp;
2736  sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2737 
2738  // Answer.
2739  return sigma;
2740 
2741 }
2742 
2743 //--------------------------------------------------------------------------
2744 
2745 // Select identity, colour and anticolour.
2746 
2747 void Sigma2fgm2Wf::setIdColAcol() {
2748 
2749  // Sign of outgoing W. Flavour of outgoing fermion.
2750  int idq = (id2 == 22) ? id1 : id2;
2751  int sign = 1 - 2 * (abs(idq)%2);
2752  if (idq < 0) sign = -sign;
2753  id4 = couplingsPtr->V2CKMpick(idq);
2754 
2755  // Flavour set up for q gamma -> W q.
2756  setId( id1, id2, 24 * sign, id4);
2757 
2758  // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2759  swapTU = (id2 == 22);
2760 
2761  // Colour flow topologies. Swap when antiquarks.
2762  if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2763  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2764  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2765  if (idq < 0) swapColAcol();
2766 
2767 }
2768 
2769 //==========================================================================
2770 
2771 // Sigma2gmgm2ffbar class.
2772 // Cross section for gamma gamma -> l lbar.
2773 
2774 //--------------------------------------------------------------------------
2775 
2776 // Initialize process wrt flavour.
2777 
2778 void Sigma2gmgm2ffbar::initProc() {
2779 
2780  // Process name.
2781  nameSave = "gamma gamma -> f fbar";
2782  if (idNew == 1) nameSave = "gamma gamma -> q qbar (uds)";
2783  if (idNew == 4) nameSave = "gamma gamma -> c cbar";
2784  if (idNew == 5) nameSave = "gamma gamma -> b bbar";
2785  if (idNew == 6) nameSave = "gamma gamma -> t tbar";
2786  if (idNew == 11) nameSave = "gamma gamma -> e+ e-";
2787  if (idNew == 13) nameSave = "gamma gamma -> mu+ mu-";
2788  if (idNew == 15) nameSave = "gamma gamma -> tau+ tau-";
2789 
2790  // Generate massive phase space, except for u+d+s.
2791  idMass = 0;
2792  if (idNew > 3) idMass = idNew;
2793 
2794  // Charge and colour factor.
2795  ef4 = 1.;
2796  if (idNew == 1) ef4 = 3. * (pow4(2./3.) + 2. * pow4(1./3.));
2797  if (idNew == 4 || idNew == 6) ef4 = 3. * pow4(2./3.);
2798  if (idNew == 5) ef4 = 3. * pow4(1./3.);
2799 
2800  // Secondary open width fraction.
2801  openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
2802 
2803 }
2804 
2805 //--------------------------------------------------------------------------
2806 
2807 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2808 
2809 void Sigma2gmgm2ffbar::sigmaKin() {
2810 
2811  // Pick current flavour for u+d+s mix by e_q^4 weights.
2812  if (idNew == 1) {
2813  double rId = 18. * rndmPtr->flat();
2814  idNow = 1;
2815  if (rId > 1.) idNow = 2;
2816  if (rId > 17.) idNow = 3;
2817  s34Avg = pow2(particleDataPtr->m0(idNow));
2818  } else {
2819  idNow = idNew;
2820  s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
2821  }
2822 
2823  // Modified Mandelstam variables for massive kinematics with m3 = m4.
2824  double tHQ = -0.5 * (sH - tH + uH);
2825  double uHQ = -0.5 * (sH + tH - uH);
2826  double tHQ2 = tHQ * tHQ;
2827  double uHQ2 = uHQ * uHQ;
2828 
2829  // Calculate kinematics dependence.
2830  if (sH < 4. * s34Avg) sigTU = 0.;
2831  else sigTU = 2. * (tHQ * uHQ - s34Avg * sH)
2832  * (tHQ2 + uHQ2 + 2. * s34Avg * sH) / (tHQ2 * uHQ2);
2833 
2834  // Answer.
2835  sigma = (M_PI / sH2) * pow2(alpEM) * ef4 * sigTU * openFracPair;
2836 
2837 }
2838 
2839 //--------------------------------------------------------------------------
2840 
2841 // Select identity, colour and anticolour.
2842 
2843 void Sigma2gmgm2ffbar::setIdColAcol() {
2844 
2845  // Flavours are trivial.
2846  setId( id1, id2, idNow, -idNow);
2847 
2848  // Colour flow in singlet state.
2849  if (idNow < 10) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
2850  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2851 
2852 }
2853 
2854 //==========================================================================
2855 
2856 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5