PYTHIA8/pythia8145/src/SigmaEW.cxx

   1 // SigmaEW.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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() * couplingsPtr->cos2thetaW());
 319
 320 }
 321
 322 //--------------------------------------------------------------------------
 323
 324 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 325
 326 void Sigma2ff2fftgmZ::sigmaKin() {
 327
 328   // Cross section part common for all incoming flavours.
 329   double sigma0 = (M_PI / sH2) * pow2(alpEM);
 330
 331   // Kinematical functions for gamma-gamma, gamma-Z and Z-Z parts.
 332   sigmagmgm = sigma0 * 2. * (sH2 + uH2) / tH2;
 333   sigmagmZ  = sigma0 * 4. * thetaWRat * sH2 / (tH * (tH - mZS));
 334   sigmaZZ   = sigma0 * 2. * pow2(thetaWRat) * sH2 / pow2(tH - mZS);
 335   if (gmZmode == 1) {sigmagmZ = 0.; sigmaZZ = 0.;}
 336   if (gmZmode == 2) {sigmagmgm = 0.; sigmagmZ = 0.;}
 337
 338 }
 339
 340 //--------------------------------------------------------------------------
 341
 342 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 343
 344 double Sigma2ff2fftgmZ::sigmaHat() {
 345
 346   // Couplings for current flavour combination.
 347   int id1Abs = abs(id1);
 348   double  e1 = couplingsPtr->ef(id1Abs);
 349   double  v1 = couplingsPtr->vf(id1Abs);
 350   double  a1 = couplingsPtr->af(id1Abs);
 351   int id2Abs = abs(id2);
 352   double  e2 = couplingsPtr->ef(id2Abs);
 353   double  v2 = couplingsPtr->vf(id2Abs);
 354   double  a2 = couplingsPtr->af(id2Abs);
 355
 356   // Distinguish same-sign and opposite-sign fermions.
 357   double epsi = (id1 * id2 > 0) ? 1. : -1.;
 358
 359   // Flavour-dependent cross section.
 360   double sigma = sigmagmgm * pow2(e1 * e2)
 361     + sigmagmZ * e1 * e2 * (v1 * v2 * (1. + uH2 / sH2)
 362       + a1 * a2 * epsi * (1. - uH2 / sH2))
 363     + sigmaZZ * ((v1*v1 + a1*a1) * (v2*v2 + a2*a2) * (1. + uH2 / sH2)
 364       + 4. * v1 * a1 * v2 * a2 * epsi * (1. - uH2 / sH2));
 365
 366   // Spin-state extra factor 2 per incoming neutrino.
 367   if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
 368   if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
 369
 370   // Answer.
 371   return sigma;
 372
 373 }
 374
 375 //--------------------------------------------------------------------------
 376
 377 // Select identity, colour and anticolour.
 378
 379 void Sigma2ff2fftgmZ::setIdColAcol() {
 380
 381   // Trivial flavours: out = in.
 382   setId( id1, id2, id1, id2);
 383
 384   // Colour flow topologies. Swap when antiquarks.
 385   if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
 386                          setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
 387   else if (abs(id1) < 9 && abs(id2) < 9)
 388                          setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 389   else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
 390   else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
 391   else                   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 392   if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
 393     swapColAcol();
 394
 395 }
 396
 397 //==========================================================================
 398
 399 // Sigma2ff2fftW class.
 400 // Cross section for f_1 f_2 -> f_3 f_4 via t-channel W+- exchange
 401 // (f is quark or lepton).
 402
 403 //--------------------------------------------------------------------------
 404
 405 // Initialize process.
 406
 407 void Sigma2ff2fftW::initProc() {
 408
 409   // Store W+- mass for propagator. Common coupling factor.
 410   mW        = particleDataPtr->m0(24);
 411   mWS       = mW*mW;
 412   thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
 413
 414 }
 415
 416 //--------------------------------------------------------------------------
 417
 418 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 419
 420 void Sigma2ff2fftW::sigmaKin() {
 421
 422   // Cross section part common for all incoming flavours.
 423   sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat)
 424     * 4. * sH2 / pow2(tH - mWS);
 425
 426 }
 427
 428 //--------------------------------------------------------------------------
 429
 430 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 431
 432 double Sigma2ff2fftW::sigmaHat() {
 433
 434   // Some flavour combinations not possible.
 435   int id1Abs = abs(id1);
 436   int id2Abs = abs(id2);
 437   if ( (id1Abs%2 == id2Abs%2 && id1 * id2 > 0)
 438     || (id1Abs%2 != id2Abs%2 && id1 * id2 < 0) ) return 0.;
 439
 440   // Basic cross section.
 441   double sigma = sigma0;
 442   if (id1 * id2 < 0) sigma *= uH2 / sH2;
 443
 444   // CKM factors for final states.
 445   sigma *= couplingsPtr->V2CKMsum(id1Abs) *  couplingsPtr->V2CKMsum(id2Abs);
 446
 447   // Spin-state extra factor 2 per incoming neutrino.
 448   if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
 449   if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
 450
 451   // Answer.
 452   return sigma;
 453
 454 }
 455
 456 //--------------------------------------------------------------------------
 457
 458 // Select identity, colour and anticolour.
 459
 460 void Sigma2ff2fftW::setIdColAcol() {
 461
 462   // Pick out-flavours by relative CKM weights.
 463   id3 = couplingsPtr->V2CKMpick(id1);
 464   id4 = couplingsPtr->V2CKMpick(id2);
 465   setId( id1, id2, id3, id4);
 466
 467   // Colour flow topologies. Swap when antiquarks.
 468   if      (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
 469                          setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
 470   else if (abs(id1) < 9 && abs(id2) < 9)
 471                          setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 472   else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
 473   else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
 474   else                   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 475   if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
 476     swapColAcol();
 477
 478 }
 479
 480
 481 //==========================================================================
 482
 483 // Sigma2qq2QqtW class.
 484 // Cross section for q q' -> Q q" via t-channel W+- exchange.
 485 // Related to Sigma2ff2ffViaW class, but with massive matrix elements.
 486
 487 //--------------------------------------------------------------------------
 488
 489 // Initialize process.
 490
 491 void Sigma2qq2QqtW::initProc() {
 492
 493   // Process name.
 494   nameSave                 = "q q -> Q q (t-channel W+-)";
 495   if (idNew == 4) nameSave = "q q -> c q (t-channel W+-)";
 496   if (idNew == 5) nameSave = "q q -> b q (t-channel W+-)";
 497   if (idNew == 6) nameSave = "q q -> t q (t-channel W+-)";
 498   if (idNew == 7) nameSave = "q q -> b' q (t-channel W+-)";
 499   if (idNew == 8) nameSave = "q q -> t' q (t-channel W+-)";
 500
 501   // Store W+- mass for propagator. Common coupling factor.
 502   mW        = particleDataPtr->m0(24);
 503   mWS       = mW*mW;
 504   thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
 505
 506   // Secondary open width fractions, relevant for top (or heavier).
 507   openFracPos = particleDataPtr->resOpenFrac(idNew);
 508   openFracNeg = particleDataPtr->resOpenFrac(-idNew);
 509
 510 }
 511
 512 //--------------------------------------------------------------------------
 513
 514 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 515
 516 void Sigma2qq2QqtW::sigmaKin() {
 517
 518   // Cross section part common for all incoming flavours.
 519   sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat) * 4. / pow2(tH - mWS);
 520
 521 }
 522
 523 //--------------------------------------------------------------------------
 524
 525 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 526
 527 double Sigma2qq2QqtW::sigmaHat() {
 528
 529   // Some flavour combinations not possible.
 530   int id1Abs  = abs(id1);
 531   int id2Abs  = abs(id2);
 532   bool diff12 = (id1Abs%2 != id2Abs%2);
 533   if ( (!diff12 && id1 * id2 > 0)
 534     || ( diff12 && id1 * id2 < 0) ) return 0.;
 535
 536   // Basic cross section.
 537   double sigma = sigma0;
 538   sigma *= (id1 * id2 > 0) ? sH * (sH - s3) : uH * (uH - s3);
 539
 540   // Secondary width if t or tbar produced on either side.
 541   double openFrac1 = (id1 > 0) ? openFracPos : openFracNeg;
 542   double openFrac2 = (id2 > 0) ? openFracPos : openFracNeg;
 543
 544   // CKM factors for final states; further impossible case.
 545   bool diff1N = (id1Abs%2 != idNew%2);
 546   bool diff2N = (id2Abs%2 != idNew%2);
 547   if (diff1N && diff2N)
 548     sigma *= ( couplingsPtr->V2CKMid(id1Abs, idNew) * openFrac1
 549              * couplingsPtr->V2CKMsum(id2Abs) + couplingsPtr->V2CKMsum(id1Abs)
 550              * couplingsPtr->V2CKMid(id2Abs, idNew) * openFrac2 );
 551   else if (diff1N)
 552     sigma *= couplingsPtr->V2CKMid(id1Abs, idNew) * openFrac1
 553              * couplingsPtr->V2CKMsum(id2Abs);
 554   else if (diff2N)
 555     sigma *= couplingsPtr->V2CKMsum(id1Abs)
 556              * couplingsPtr->V2CKMid(id2Abs, idNew) * openFrac2;
 557   else sigma = 0.;
 558
 559   // Spin-state extra factor 2 per incoming neutrino.
 560   if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
 561   if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
 562
 563   // Answer.
 564   return  sigma;
 565
 566 }
 567
 568 //--------------------------------------------------------------------------
 569
 570 // Select identity, colour and anticolour.
 571
 572 void Sigma2qq2QqtW::setIdColAcol() {
 573
 574   // For topologies like d dbar -> (t/c/u) (t/c/u)bar pick side.
 575   int id1Abs = abs(id1);
 576   int id2Abs = abs(id2);
 577   int side   = 1;
 578   if ( (id1Abs + idNew)%2 == 1 && (id2Abs + idNew)%2 == 1 ) {
 579     double prob1 = couplingsPtr->V2CKMid(id1Abs, idNew)
 580                    * couplingsPtr->V2CKMsum(id2Abs);
 581     prob1       *= (id1 > 0) ? openFracPos : openFracNeg;
 582     double prob2 = couplingsPtr->V2CKMid(id2Abs, idNew)
 583                    * couplingsPtr->V2CKMsum(id1Abs);
 584     prob2       *= (id2 > 0) ? openFracPos : openFracNeg;
 585     if (prob2 > rndmPtr->flat() * (prob1 + prob2)) side = 2;
 586   }
 587   else if ((id2Abs + idNew)%2 == 1) side = 2;
 588
 589   // Pick out-flavours by relative CKM weights.
 590   if (side == 1) {
 591     // q q' -> t q" : correct order from start.
 592     id3 = (id1 > 0) ? idNew : -idNew;
 593     id4 = couplingsPtr->V2CKMpick(id2);
 594     setId( id1, id2, id3, id4);
 595   } else {
 596     // q q' -> q" t : stored as t q" so swap tHat <-> uHat.
 597     swapTU = true;
 598     id3 = couplingsPtr->V2CKMpick(id1);
 599     id4 = (id2 > 0) ? idNew : -idNew;
 600     setId( id1, id2, id4, id3);
 601   }
 602
 603   // Colour flow topologies. Swap when antiquarks on side 1.
 604   if      (side == 1 && id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
 605   else if              (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
 606   else if (side == 1)                  setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 607   else                                 setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
 608   if (id1 < 0) swapColAcol();
 609
 610 }
 611
 612 //--------------------------------------------------------------------------
 613
 614 // Evaluate weight for decay angles of W in top decay.
 615
 616 double Sigma2qq2QqtW::weightDecay( Event& process, int iResBeg,
 617   int iResEnd) {
 618
 619   // For top decay hand over to standard routine, else done.
 620   if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
 621        return weightTopDecay( process, iResBeg, iResEnd);
 622   else return 1.;
 623
 624 }
 625
 626 //==========================================================================
 627
 628 // Sigma1ffbar2gmZ class.
 629 // Cross section for f fbar -> gamma*/Z0 (f is quark or lepton).
 630
 631 //--------------------------------------------------------------------------
 632
 633 // Initialize process.
 634
 635 void Sigma1ffbar2gmZ::initProc() {
 636
 637   // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
 638   gmZmode     = settingsPtr->mode("WeakZ0:gmZmode");
 639
 640   // Store Z0 mass and width for propagator.
 641   mRes        = particleDataPtr->m0(23);
 642   GammaRes    = particleDataPtr->mWidth(23);
 643   m2Res       = mRes*mRes;
 644   GamMRat     = GammaRes / mRes;
 645   thetaWRat   = 1. / (16. * couplingsPtr->sin2thetaW() * couplingsPtr->cos2thetaW());
 646
 647   // Set pointer to particle properties and decay table.
 648   particlePtr = particleDataPtr->particleDataEntryPtr(23);
 649
 650 }
 651
 652 //--------------------------------------------------------------------------
 653
 654 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 655
 656 void Sigma1ffbar2gmZ::sigmaKin() {
 657
 658   // Common coupling factors.
 659   double colQ = 3. * (1. + alpS / M_PI);
 660
 661   // Reset quantities to sum. Declare variables in loop.
 662   gamSum = 0.;
 663   intSum = 0.;
 664   resSum = 0.;
 665   int    idAbs, onMode;
 666   double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
 667
 668   // Loop over all Z0 decay channels.
 669   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
 670     idAbs = abs( particlePtr->channel(i).product(0) );
 671
 672     // Only contributions from three fermion generations, except top.
 673     if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
 674       mf = particleDataPtr->m0(idAbs);
 675
 676       // Check that above threshold. Phase space.
 677       if (mH > 2. * mf + MASSMARGIN) {
 678         mr    = pow2(mf / mH);
 679         betaf = sqrtpos(1. - 4. * mr);
 680         psvec = betaf * (1. + 2. * mr);
 681         psaxi = pow3(betaf);
 682
 683         // Combine phase space with couplings.
 684         ef2    = couplingsPtr->ef2(idAbs) * psvec;
 685         efvf   = couplingsPtr->efvf(idAbs) * psvec;
 686         vf2af2 = couplingsPtr->vf2(idAbs) * psvec
 687                + couplingsPtr->af2(idAbs) * psaxi;
 688         colf   = (idAbs < 6) ? colQ : 1.;
 689
 690         // Store sum of combinations. For outstate only open channels.
 691         onMode = particlePtr->channel(i).onMode();
 692         if (onMode == 1 || onMode == 2) {
 693           gamSum += colf * ef2;
 694           intSum += colf * efvf;
 695           resSum += colf * vf2af2;
 696         }
 697
 698       // End loop over fermions.
 699       }
 700     }
 701   }
 702
 703   // Calculate prefactors for gamma/interference/Z0 cross section terms.
 704   gamProp = 4. * M_PI * pow2(alpEM) / (3. * sH);
 705   intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
 706           / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 707   resProp = gamProp * pow2(thetaWRat * sH)
 708           / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 709
 710   // Optionally only keep gamma* or Z0 term.
 711   if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
 712   if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
 713
 714 }
 715
 716 //--------------------------------------------------------------------------
 717
 718 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 719
 720 double Sigma1ffbar2gmZ::sigmaHat() {
 721
 722   // Combine gamma, interference and Z0 parts.
 723   int idAbs = abs(id1);
 724   double sigma = couplingsPtr->ef2(idAbs)    * gamProp * gamSum
 725                + couplingsPtr->efvf(idAbs)   * intProp * intSum
 726                + couplingsPtr->vf2af2(idAbs) * resProp * resSum;
 727
 728   // Colour factor. Answer.
 729   if (idAbs < 9) sigma /= 3.;
 730   return sigma;
 731
 732 }
 733
 734 //--------------------------------------------------------------------------
 735
 736 // Select identity, colour and anticolour.
 737
 738 void Sigma1ffbar2gmZ::setIdColAcol() {
 739
 740   // Flavours trivial.
 741   setId( id1, id2, 23);
 742
 743   // Colour flow topologies. Swap when antiquarks.
 744   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
 745   else              setColAcol( 0, 0, 0, 0, 0, 0);
 746   if (id1 < 0) swapColAcol();
 747
 748 }
 749
 750 //--------------------------------------------------------------------------
 751
 752 // Evaluate weight for gamma*/Z0 decay angle.
 753
 754 double Sigma1ffbar2gmZ::weightDecay( Event& process, int iResBeg,
 755   int iResEnd) {
 756
 757   // Z should sit in entry 5.
 758   if (iResBeg != 5 || iResEnd != 5) return 1.;
 759
 760   // Couplings for in- and out-flavours.
 761   int idInAbs  = process[3].idAbs();
 762   double ei    = couplingsPtr->ef(idInAbs);
 763   double vi    = couplingsPtr->vf(idInAbs);
 764   double ai    = couplingsPtr->af(idInAbs);
 765   int idOutAbs = process[6].idAbs();
 766   double ef    = couplingsPtr->ef(idOutAbs);
 767   double vf    = couplingsPtr->vf(idOutAbs);
 768   double af    = couplingsPtr->af(idOutAbs);
 769
 770   // Phase space factors. (One power of beta left out in formulae.)
 771   double mf    = process[6].m();
 772   double mr    = mf*mf / sH;
 773   double betaf = sqrtpos(1. - 4. * mr);
 774
 775   // Coefficients of angular expression.
 776   double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
 777     + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
 778   double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
 779     + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
 780   double coefAsym = betaf * ( ei * ai * intProp * ef * af
 781     + 4. * vi * ai * resProp * vf * af );
 782
 783   // Flip asymmetry for in-fermion + out-antifermion.
 784   if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
 785
 786   // Reconstruct decay angle and weight for it.
 787   double cosThe = (process[3].p() - process[4].p())
 788     * (process[7].p() - process[6].p()) / (sH * betaf);
 789   double wtMax = 2. * (coefTran + abs(coefAsym));
 790   double wt    = coefTran * (1. + pow2(cosThe))
 791      + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
 792
 793   // Done.
 794   return (wt / wtMax);
 795
 796 }
 797
 798 //==========================================================================
 799
 800 // Sigma1ffbar2W class.
 801 // Cross section for f fbar' -> W+- (f is quark or lepton).
 802
 803 //--------------------------------------------------------------------------
 804
 805 // Initialize process.
 806
 807 void Sigma1ffbar2W::initProc() {
 808
 809   // Store W+- mass and width for propagator.
 810   mRes     = particleDataPtr->m0(24);
 811   GammaRes = particleDataPtr->mWidth(24);
 812   m2Res    = mRes*mRes;
 813   GamMRat  = GammaRes / mRes;
 814   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
 815
 816   // Set pointer to particle properties and decay table.
 817   particlePtr = particleDataPtr->particleDataEntryPtr(24);
 818
 819 }
 820
 821 //--------------------------------------------------------------------------
 822
 823 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 824
 825 void Sigma1ffbar2W::sigmaKin() {
 826
 827   // Set up Breit-Wigner. Cross section for W+ and W- separately.
 828   double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 829   double preFac = alpEM * thetaWRat * mH;
 830   sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(24, mH);
 831   sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-24, mH);
 832
 833 }
 834
 835 //--------------------------------------------------------------------------
 836
 837 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 838
 839 double Sigma1ffbar2W::sigmaHat() {
 840
 841   // Secondary width for W+ or W-. CKM and colour factors.
 842   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
 843   double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
 844   if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
 845
 846   // Answer.
 847   return sigma;
 848
 849 }
 850
 851 //--------------------------------------------------------------------------
 852
 853 // Select identity, colour and anticolour.
 854
 855 void Sigma1ffbar2W::setIdColAcol() {
 856
 857   // Sign of outgoing W.
 858   int sign          = 1 - 2 * (abs(id1)%2);
 859   if (id1 < 0) sign = -sign;
 860   setId( id1, id2, 24 * sign);
 861
 862   // Colour flow topologies. Swap when antiquarks.
 863   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
 864   else              setColAcol( 0, 0, 0, 0, 0, 0);
 865   if (id1 < 0) swapColAcol();
 866
 867 }
 868
 869 //--------------------------------------------------------------------------
 870
 871 // Evaluate weight for W decay angle.
 872
 873 double Sigma1ffbar2W::weightDecay( Event& process, int iResBeg,
 874   int iResEnd) {
 875
 876   // W should sit in entry 5.
 877   if (iResBeg != 5 || iResEnd != 5) return 1.;
 878
 879   // Phase space factors.
 880   double mr1    = pow2(process[6].m()) / sH;
 881   double mr2    = pow2(process[7].m()) / sH;
 882   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 883
 884   // Sign of asymmetry.
 885   double eps    = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
 886
 887   // Reconstruct decay angle and weight for it.
 888   double cosThe = (process[3].p() - process[4].p())
 889     * (process[7].p() - process[6].p()) / (sH * betaf);
 890   double wtMax  = 4.;
 891   double wt     = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
 892
 893   // Done.
 894   return (wt / wtMax);
 895
 896 }
 897
 898 //==========================================================================
 899
 900 // Sigma2ffbar2ffbarsgm class.
 901 // Cross section f fbar -> gamma* -> f' fbar', for multiple interactions.
 902
 903
 904 //--------------------------------------------------------------------------
 905
 906 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 907
 908 void Sigma2ffbar2ffbarsgm::sigmaKin() {
 909
 910   // Pick new flavour. Allow three leptons and five quarks.
 911   double colQ     = 1. + (alpS / M_PI);
 912   double flavWt   = 3. + colQ * 11. / 3.;
 913   double flavRndm = rndmPtr->flat() * flavWt;
 914   if (flavRndm < 3.) {
 915     if      (flavRndm < 1.) idNew = 11;
 916     else if (flavRndm < 2.) idNew = 13;
 917     else                    idNew = 15;
 918   } else {
 919     flavRndm = 3. * (flavWt - 3.) / colQ;
 920     if      (flavRndm <  4.) idNew = 2;
 921     else if (flavRndm <  8.) idNew = 4;
 922     else if (flavRndm <  9.) idNew = 1;
 923     else if (flavRndm < 10.) idNew = 3;
 924     else                     idNew = 5;
 925   }
 926   double mNew  = particleDataPtr->m0(idNew);
 927   double m2New = mNew*mNew;
 928
 929   // Calculate kinematics dependence. Give correct mass factors for
 930   // tHat, uHat defined as if massless kinematics, d(sigma)/d(Omega)
 931   // = beta (1 + cos^2(theta) + (1 - beta^2) sin^2(theta)).
 932   // Special case related to phase space form in multiple interactions.
 933   double sigS = 0.;
 934   if (sH > 4. * m2New) {
 935     double beta = sqrt(1. - 4. * m2New / sH);
 936     sigS = beta * (2.* (tH2 + uH2) + 4. * (1. - beta * beta) * tH * uH)
 937       / sH2;
 938   }
 939
 940   // Answer is proportional to number of outgoing flavours.
 941   sigma0 = (M_PI/sH2) * pow2(alpEM) * sigS * flavWt;
 942
 943 }
 944
 945 //--------------------------------------------------------------------------
 946
 947 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 948
 949 double Sigma2ffbar2ffbarsgm::sigmaHat() {
 950
 951   // Charge and colour factors.
 952   double eNow  = couplingsPtr->ef( abs(id1) );
 953   double sigma = sigma0 * pow2(eNow);
 954   if (abs(id1) < 9) sigma /= 3.;
 955
 956   // Answer.
 957   return sigma;
 958
 959 }
 960
 961 //--------------------------------------------------------------------------
 962
 963 // Select identity, colour and anticolour.
 964
 965 void Sigma2ffbar2ffbarsgm::setIdColAcol() {
 966
 967   // Set outgoing flavours.
 968   id3 = (id1 > 0) ? idNew : -idNew;
 969   setId( id1, id2, id3, -id3);
 970
 971   // Colour flow topologies. Swap when antiquarks.
 972   if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
 973   else if (abs(id1) < 9)         setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
 974   else if (idNew < 9)            setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
 975   else                           setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 976   if (id1 < 0) swapColAcol();
 977
 978 }
 979
 980 //==========================================================================
 981
 982 // Sigma2ffbar2FFbarsgmZ class.
 983 // Cross section f fbar -> gamma*/Z0 -> F Fbar.
 984
 985 //--------------------------------------------------------------------------
 986
 987 // Initialize process.
 988
 989 void Sigma2ffbar2FFbarsgmZ::initProc() {
 990
 991   // Process name.
 992   nameSave                 = "f fbar -> F Fbar (s-channel gamma*/Z0)";
 993   if (idNew == 4) nameSave = "f fbar -> c cbar (s-channel gamma*/Z0)";
 994   if (idNew == 5) nameSave = "f fbar -> b bbar (s-channel gamma*/Z0)";
 995   if (idNew == 6) nameSave = "f fbar -> t tbar (s-channel gamma*/Z0)";
 996   if (idNew == 7) nameSave = "f fbar -> b' b'bar (s-channel gamma*/Z0)";
 997   if (idNew == 8) nameSave = "f fbar -> t' t'bar (s-channel gamma*/Z0)";
 998   if (idNew == 15) nameSave = "f fbar -> tau+ tau- (s-channel gamma*/Z0)";
 999   if (idNew == 17) nameSave = "f fbar -> tau'+ tau'- (s-channel gamma*/Z0)";
1000   if (idNew == 18)
1001     nameSave   = "f fbar -> nu'_tau nu'bar_tau (s-channel gamma*/Z0)";
1002
1003   // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1004   gmZmode      = settingsPtr->mode("WeakZ0:gmZmode");
1005
1006   // Store Z0 mass and width for propagator.
1007   mRes         = particleDataPtr->m0(23);
1008   GammaRes     = particleDataPtr->mWidth(23);
1009   m2Res        = mRes*mRes;
1010   GamMRat      = GammaRes / mRes;
1011   thetaWRat    = 1. / (16. * couplingsPtr->sin2thetaW()
1012                  * couplingsPtr->cos2thetaW());
1013
1014   // Store couplings of F.
1015   ef           = couplingsPtr->ef(idNew);
1016   vf           = couplingsPtr->vf(idNew);
1017   af           = couplingsPtr->af(idNew);
1018
1019   // Secondary open width fraction, relevant for top (or heavier).
1020   openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
1021
1022 }
1023
1024 //--------------------------------------------------------------------------
1025
1026 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1027
1028 void Sigma2ffbar2FFbarsgmZ::sigmaKin() {
1029
1030   // Check that above threshold.
1031   isPhysical     = true;
1032   if (mH < m3 + m4 + MASSMARGIN) {
1033     isPhysical   = false;
1034     return;
1035   }
1036
1037   // Define average F, Fbar mass so same beta. Phase space.
1038   double s34Avg  = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1039   mr             = s34Avg / sH;
1040   betaf          = sqrtpos(1. - 4. * mr);
1041
1042   // Final-state colour factor.
1043   double colF    = (idNew < 9) ? 3. * (1. + alpS / M_PI) : 1.;
1044
1045   // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1046   cosThe         = (tH - uH) / (betaf * sH);
1047
1048   // Calculate prefactors for gamma/interference/Z0 cross section terms.
1049   gamProp = colF * M_PI * pow2(alpEM) / sH2;
1050   intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1051           / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1052   resProp = gamProp * pow2(thetaWRat * sH)
1053           / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1054
1055   // Optionally only keep gamma* or Z0 term.
1056   if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1057   if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1058
1059 }
1060
1061 //--------------------------------------------------------------------------
1062
1063 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1064
1065 double Sigma2ffbar2FFbarsgmZ::sigmaHat() {
1066
1067   // Fail if below threshold.
1068   if (!isPhysical) return 0.;
1069
1070   // Couplings for in-flavours.
1071   int idAbs       = abs(id1);
1072   double ei       = couplingsPtr->ef(idAbs);
1073   double vi       = couplingsPtr->vf(idAbs);
1074   double ai       = couplingsPtr->af(idAbs);
1075
1076   // Coefficients of angular expression.
1077   double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
1078     + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
1079   double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
1080     + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
1081   double coefAsym = betaf * ( ei * ai * intProp * ef * af
1082     + 4. * vi * ai * resProp * vf * af );
1083
1084   // Combine gamma, interference and Z0 parts.
1085   double sigma    = coefTran * (1. + pow2(cosThe))
1086    + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
1087
1088   // Top: corrections for closed decay channels.
1089   sigma *= openFracPair;
1090
1091   // Initial-state colour factor. Answer.
1092   if (idAbs < 9) sigma /= 3.;
1093   return sigma;
1094
1095 }
1096
1097 //--------------------------------------------------------------------------
1098
1099 // Select identity, colour and anticolour.
1100
1101 void Sigma2ffbar2FFbarsgmZ::setIdColAcol() {
1102
1103   // Set outgoing flavours.
1104   id3 = (id1 > 0) ? idNew : -idNew;
1105   setId( id1, id2, id3, -id3);
1106
1107   // Colour flow topologies. Swap when antiquarks.
1108   if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1109   else if (abs(id1) < 9)         setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1110   else if (idNew < 9)            setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1111   else                           setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1112   if (id1 < 0) swapColAcol();
1113
1114 }
1115
1116 //--------------------------------------------------------------------------
1117
1118 // Evaluate weight for decay angles of W in top decay.
1119
1120 double Sigma2ffbar2FFbarsgmZ::weightDecay( Event& process, int iResBeg,
1121   int iResEnd) {
1122
1123   // For top decay hand over to standard routine, else done.
1124   if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1125        return weightTopDecay( process, iResBeg, iResEnd);
1126   else return 1.;
1127
1128 }
1129
1130 //==========================================================================
1131
1132 // Sigma2ffbar2FfbarsW class.
1133 // Cross section f fbar' -> W+- -> F fbar".
1134
1135 //--------------------------------------------------------------------------
1136
1137 // Initialize process.
1138
1139 void Sigma2ffbar2FfbarsW::initProc() {
1140
1141   // Process name.
1142   nameSave                 = "f fbar -> F fbar (s-channel W+-)";
1143   if (idNew == 4) nameSave = "f fbar -> c qbar (s-channel W+-)";
1144   if (idNew == 5) nameSave = "f fbar -> b qbar (s-channel W+-)";
1145   if (idNew == 6) nameSave = "f fbar -> t qbar (s-channel W+-)";
1146   if (idNew == 7) nameSave = "f fbar -> b' qbar (s-channel W+-)";
1147   if (idNew == 8) nameSave = "f fbar -> t' qbar (s-channel W+-)";
1148   if (idNew == 7 && idNew2 == 6)
1149     nameSave = "f fbar -> b' tbar (s-channel W+-)";
1150   if (idNew == 8 && idNew2 == 7)
1151     nameSave = "f fbar -> t' b'bar (s-channel W+-)";
1152   if (idNew == 15 || idNew == 16)
1153     nameSave = "f fbar -> tau nu_taubar (s-channel W+-)";
1154   if (idNew == 17 || idNew == 18)
1155     nameSave = "f fbar -> tau'  nu'_taubar (s-channel W+-)";
1156
1157   // Store W+- mass and width for propagator.
1158   mRes      = particleDataPtr->m0(24);
1159   GammaRes  = particleDataPtr->mWidth(24);
1160   m2Res     = mRes*mRes;
1161   GamMRat   = GammaRes / mRes;
1162   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1163
1164   // For t/t' want to use at least b mass.
1165   idPartner = idNew2;
1166   if ( (idNew == 6 || idNew == 8) && idNew2 == 0 ) idPartner = 5;
1167
1168   // Sum of CKM weights for quarks.
1169   V2New     = (idNew < 9) ? couplingsPtr->V2CKMsum(idNew) : 1.;
1170   if (idNew2 != 0) V2New = couplingsPtr->V2CKMid(idNew, idNew2);
1171
1172   // Secondary open width fractions, relevant for top or heavier.
1173   openFracPos = particleDataPtr->resOpenFrac( idNew, -idNew2);
1174   openFracNeg = particleDataPtr->resOpenFrac(-idNew,  idNew2);
1175
1176 }
1177
1178 //--------------------------------------------------------------------------
1179
1180 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1181
1182 void Sigma2ffbar2FfbarsW::sigmaKin() {
1183
1184   // Check that above threshold.
1185   isPhysical    = true;
1186   if (mH < m3 + m4 + MASSMARGIN) {
1187     isPhysical  = false;
1188     return;
1189   }
1190
1191   // Phase space factors.
1192   double mr1    = s3 / sH;
1193   double mr2    = s4 / sH;
1194   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
1195
1196   // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1197   double cosThe = (tH - uH) / (betaf * sH);
1198
1199   // Set up Breit-Wigner and in- and out-widths.
1200   double sigBW  = 9. * M_PI * pow2(alpEM * thetaWRat)
1201                 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1202
1203   // Initial-state colour factor.
1204   double colF   = (idNew < 9) ? 3. * (1. + alpS / M_PI) * V2New : 1.;
1205
1206   // Angular dependence.
1207   double wt     = pow2(1. + betaf * cosThe) - pow2(mr1 - mr2);
1208
1209   // Temporary answer.
1210   sigma0        = sigBW * colF * wt;
1211
1212 }
1213
1214 //--------------------------------------------------------------------------
1215
1216 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1217
1218 double Sigma2ffbar2FfbarsW::sigmaHat() {
1219
1220   // Fail if below threshold.
1221   if (!isPhysical) return 0.;
1222
1223   // CKM and colour factors.
1224   double sigma = sigma0;
1225   if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1226
1227   // Correction for secondary width in top (or heavier) decay.
1228   int idSame = ((abs(id1) + idNew)%2 == 0) ? id1 : id2;
1229   sigma *= (idSame > 0) ? openFracPos : openFracNeg;
1230
1231   // Answer.
1232   return sigma;
1233
1234 }
1235
1236 //--------------------------------------------------------------------------
1237
1238 // Select identity, colour and anticolour.
1239
1240 void Sigma2ffbar2FfbarsW::setIdColAcol() {
1241
1242   // Set outgoing flavours.
1243   id3 = idNew;
1244   id4 = (idNew2 != 0) ? idNew2 : couplingsPtr->V2CKMpick(idNew);
1245   if (idNew%2 == 0) {
1246     int idInUp = (abs(id1)%2 == 0) ? id1 : id2;
1247     if (idInUp > 0) id4 = -id4;
1248     else            id3 = -id3;
1249   } else {
1250     int idInDn = (abs(id1)%2 == 1) ? id1 : id2;
1251     if (idInDn > 0) id4 = -id4;
1252     else            id3 = -id3;
1253   }
1254   setId( id1, id2, id3, id4);
1255
1256   // Swap tHat and uHat for fbar' f -> F f".
1257   if (id1 * id3 < 0) swapTU = true;
1258
1259   // Colour flow topologies. Swap when antiquarks.
1260   if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1261   else if (abs(id1) < 9)         setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1262   else if (idNew < 9)            setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1263   else                           setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1264   if (id1 < 0) swapCol12();
1265   if (id3 < 0) swapCol34();
1266
1267 }
1268
1269 //--------------------------------------------------------------------------
1270
1271 // Evaluate weight for decay angles of W in top decay.
1272
1273 double Sigma2ffbar2FfbarsW::weightDecay( Event& process, int iResBeg,
1274   int iResEnd) {
1275
1276   // For top decay hand over to standard routine, else done.
1277   if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1278        return weightTopDecay( process, iResBeg, iResEnd);
1279   else return 1.;
1280
1281 }
1282
1283 //==========================================================================
1284
1285 // Sigma2ffbargmZWgmZW class.
1286 // Collects common methods for f fbar ->  gamma*/Z0/W+- gamma*/Z0/W-+.
1287
1288 //--------------------------------------------------------------------------
1289
1290 // Calculate and store internal products.
1291
1292 void Sigma2ffbargmZWgmZW::setupProd( Event& process, int i1, int i2,
1293   int i3, int i4, int i5, int i6) {
1294
1295   // Store incoming and outgoing momenta,
1296   pRot[1] = process[i1].p();
1297   pRot[2] = process[i2].p();
1298   pRot[3] = process[i3].p();
1299   pRot[4] = process[i4].p();
1300   pRot[5] = process[i5].p();
1301   pRot[6] = process[i6].p();
1302
1303   // Do random rotation to avoid accidental zeroes in HA expressions.
1304   bool smallPT = false;
1305   do {
1306     smallPT = false;
1307     double thetaNow = acos(2. * rndmPtr->flat() - 1.);
1308     double phiNow   = 2. * M_PI * rndmPtr->flat();
1309     for (int i = 1; i <= 6; ++i) {
1310       pRot[i].rot( thetaNow, phiNow);
1311       if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
1312     }
1313   } while (smallPT);
1314
1315   // Calculate internal products.
1316   for (int i = 1; i < 6; ++i) {
1317     for (int j = i + 1; j <= 6; ++j) {
1318       hA[i][j] =
1319           sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
1320         / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
1321         - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
1322         / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
1323       hC[i][j] = conj( hA[i][j] );
1324       if (i <= 2) {
1325         hA[i][j] *= complex( 0., 1.);
1326         hC[i][j] *= complex( 0., 1.);
1327       }
1328       hA[j][i] = - hA[i][j];
1329       hC[j][i] = - hC[i][j];
1330     }
1331   }
1332
1333 }
1334
1335 //--------------------------------------------------------------------------
1336
1337 // Evaluate the F function of Gunion and Kunszt.
1338
1339 complex Sigma2ffbargmZWgmZW::fGK(int j1, int j2, int j3, int j4, int j5,
1340   int j6) {
1341
1342   return 4. * hA[j1][j3] * hC[j2][j6]
1343          * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
1344
1345 }
1346
1347 //--------------------------------------------------------------------------
1348
1349 // Evaluate the Xi function of Gunion and Kunszt.
1350
1351 double Sigma2ffbargmZWgmZW::xiGK( double tHnow, double uHnow) {
1352
1353   return - 4. * s3 * s4 + tHnow * (3. * tHnow + 4. * uHnow)
1354          + tHnow * tHnow * ( tHnow * uHnow / (s3 * s4)
1355            - 2. * (1. / s3 + 1./s4) * (tHnow + uHnow)
1356            + 2. * (s3 / s4 + s4 / s3) );
1357
1358 }
1359
1360 //--------------------------------------------------------------------------
1361
1362 // Evaluate the Xj function of Gunion and Kunszt.
1363
1364 double Sigma2ffbargmZWgmZW::xjGK( double tHnow, double uHnow) {
1365
1366   return 8. * pow2(s3 + s4) - 8. * (s3 + s4) * (tHnow + uHnow)
1367          - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
1368            / (s3 * s4) - 2. * (1. / s3 + 1. / s4) * (tHnow + uHnow)
1369            + 2. * (s3 / s4 + s4 / s3) );
1370
1371 }
1372
1373 //==========================================================================
1374
1375 // Sigma2ffbar2gmZgmZ class.
1376 // Cross section for f fbar -> gamma*/Z0 gamma*/Z0 (f is quark or lepton).
1377
1378 //--------------------------------------------------------------------------
1379
1380 // Initialize process.
1381
1382 void Sigma2ffbar2gmZgmZ::initProc() {
1383
1384   // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1385   gmZmode     = settingsPtr->mode("WeakZ0:gmZmode");
1386
1387   // Store Z0 mass and width for propagator.
1388   mRes        = particleDataPtr->m0(23);
1389   GammaRes    = particleDataPtr->mWidth(23);
1390   m2Res       = mRes*mRes;
1391   GamMRat     = GammaRes / mRes;
1392   thetaWRat   = 1. / (16. * couplingsPtr->sin2thetaW()
1393                 * couplingsPtr->cos2thetaW());
1394
1395   // Set pointer to particle properties and decay table.
1396   particlePtr = particleDataPtr->particleDataEntryPtr(23);
1397
1398 }
1399
1400 //--------------------------------------------------------------------------
1401
1402 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1403
1404 void Sigma2ffbar2gmZgmZ::sigmaKin() {
1405
1406   // Cross section part common for all incoming flavours.
1407   sigma0 = (M_PI / sH2) * pow2(alpEM) * 0.5
1408     * ( (tH2 + uH2 + 2. * (s3 + s4) * sH) / (tH * uH)
1409     - s3 * s4 * (1./tH2 + 1./uH2) );
1410
1411   // Common coupling factors at the resonance masses
1412   double alpEM3 = couplingsPtr->alphaEM(s3);
1413   double alpS3  = couplingsPtr->alphaS(s3);
1414   double colQ3  = 3. * (1. + alpS3 / M_PI);
1415   double alpEM4 = couplingsPtr->alphaEM(s4);
1416   double alpS4  = couplingsPtr->alphaS(s4);
1417   double colQ4  = 3. * (1. + alpS4 / M_PI);
1418
1419   // Reset quantities to sum. Declare variables in loop.
1420   gamSum3 = 0.;
1421   intSum3 = 0.;
1422   resSum3 = 0.;
1423   gamSum4 = 0.;
1424   intSum4 = 0.;
1425   resSum4 = 0.;
1426   int    onMode;
1427   double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
1428
1429   // Loop over all Z0 decay channels.
1430   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1431     int idAbs = abs( particlePtr->channel(i).product(0) );
1432
1433     // Only contributions from three fermion generations, except top.
1434     if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
1435       mf     = particleDataPtr->m0(idAbs);
1436       onMode = particlePtr->channel(i).onMode();
1437
1438       // First Z0: check that above threshold. Phase space.
1439       if (m3 > 2. * mf + MASSMARGIN) {
1440         mr    = pow2(mf / m3);
1441         betaf = sqrtpos(1. - 4. * mr);
1442         psvec = betaf * (1. + 2. * mr);
1443         psaxi  = pow3(betaf);
1444
1445         // First Z0: combine phase space with couplings.
1446         ef2    = couplingsPtr->ef2(idAbs) * psvec;
1447         efvf   = couplingsPtr->efvf(idAbs) * psvec;
1448         vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1449                + couplingsPtr->af2(idAbs) * psaxi;
1450         colf   = (idAbs < 6) ? colQ3 : 1.;
1451
1452         // First Z0: store sum of combinations for open outstate channels.
1453         if (onMode == 1 || onMode == 2) {
1454           gamSum3 += colf * ef2;
1455           intSum3 += colf * efvf;
1456           resSum3 += colf * vf2af2;
1457         }
1458       }
1459
1460       // Second Z0: check that above threshold. Phase space.
1461       if (m4 > 2. * mf + MASSMARGIN) {
1462         mr    = pow2(mf / m4);
1463         betaf = sqrtpos(1. - 4. * mr);
1464         psvec = betaf * (1. + 2. * mr);
1465         psaxi = pow3(betaf);
1466
1467         // Second Z0: combine phase space with couplings.
1468         ef2    = couplingsPtr->ef2(idAbs) * psvec;
1469         efvf   = couplingsPtr->efvf(idAbs) * psvec;
1470         vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1471                + couplingsPtr->af2(idAbs) * psaxi;
1472         colf   = (idAbs < 6) ? colQ4 : 1.;
1473
1474         // Second Z0: store sum of combinations for open outstate channels.
1475         if (onMode == 1 || onMode == 2) {
1476           gamSum4 += colf * ef2;
1477           intSum4 += colf * efvf;
1478           resSum4 += colf * vf2af2;
1479         }
1480       }
1481
1482     // End loop over fermions.
1483     }
1484   }
1485
1486   // First Z0: calculate prefactors for gamma/interference/Z0 terms.
1487   gamProp3 = 4. * alpEM3 / (3. * M_PI * s3);
1488   intProp3 = gamProp3 * 2. * thetaWRat * s3 * (s3 - m2Res)
1489            / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1490   resProp3 = gamProp3 * pow2(thetaWRat * s3)
1491            / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1492
1493   // First Z0: optionally only keep gamma* or Z0 term.
1494   if (gmZmode == 1) {intProp3 = 0.; resProp3 = 0.;}
1495   if (gmZmode == 2) {gamProp3 = 0.; intProp3 = 0.;}
1496
1497   // Second Z0: calculate prefactors for gamma/interference/Z0 terms.
1498   gamProp4 = 4. * alpEM4 / (3. * M_PI * s4);
1499   intProp4 = gamProp4 * 2. * thetaWRat * s4 * (s4 - m2Res)
1500            / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1501   resProp4 = gamProp4 * pow2(thetaWRat * s4)
1502            / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1503
1504   // Second Z0: optionally only keep gamma* or Z0 term.
1505   if (gmZmode == 1) {intProp4 = 0.; resProp4 = 0.;}
1506   if (gmZmode == 2) {gamProp4 = 0.; intProp4 = 0.;}
1507
1508 }
1509
1510 //--------------------------------------------------------------------------
1511
1512 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1513
1514 double Sigma2ffbar2gmZgmZ::sigmaHat() {
1515
1516   // Charge/2, left- and righthanded couplings for in-fermion.
1517   int idAbs = abs(id1);
1518   double ei = 0.5 * couplingsPtr->ef(idAbs);
1519   double li =       couplingsPtr->lf(idAbs);
1520   double ri =       couplingsPtr->rf(idAbs);
1521
1522   // Combine left/right gamma, interference and Z0 parts for each Z0.
1523   double left3  = ei * ei * gamProp3 * gamSum3
1524                 + ei * li * intProp3 * intSum3
1525                 + li * li * resProp3 * resSum3;
1526   double right3 = ei * ei * gamProp3 * gamSum3
1527                 + ei * ri * intProp3 * intSum3
1528                 + ri * ri * resProp3 * resSum3;
1529   double left4  = ei * ei * gamProp4 * gamSum4
1530                 + ei * li * intProp4 * intSum4
1531                 + li * li * resProp4 * resSum4;
1532   double right4 = ei * ei * gamProp4 * gamSum4
1533                 + ei * ri * intProp4 * intSum4
1534                 + ri * ri * resProp4 * resSum4;
1535
1536   // Combine left- and right-handed couplings for the two Z0's.
1537   double sigma = sigma0 * (left3 * left4 + right3 * right4);
1538
1539   // Correct for the running-width Z0 propagators weight in PhaseSpace.
1540   sigma /= (runBW3 * runBW4);
1541
1542   // Initial-state colour factor. Answer.
1543   if (idAbs < 9) sigma /= 3.;
1544   return  sigma;
1545
1546 }
1547
1548 //--------------------------------------------------------------------------
1549
1550 // Select identity, colour and anticolour.
1551
1552 void Sigma2ffbar2gmZgmZ::setIdColAcol() {
1553
1554   // Flavours trivial.
1555   setId( id1, id2, 23, 23);
1556
1557   // Colour flow topologies. Swap when antiquarks.
1558   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1559   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1560   if (id1 < 0) swapColAcol();
1561
1562 }
1563
1564 //--------------------------------------------------------------------------
1565
1566 // Evaluate correlated decay flavours of the two gamma*/Z0.
1567 // Unique complication, caused by gamma*/Z0 mix different left/right.
1568
1569 double Sigma2ffbar2gmZgmZ::weightDecayFlav( Event& process) {
1570
1571   // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1572   i1 = (process[3].id() < 0) ? 3 : 4;
1573   i2 = 7 - i1;
1574   i3 = (process[7].id() > 0) ? 7 : 8;
1575   i4 = 15 - i3;
1576   i5 = (process[9].id() > 0) ? 9 : 10;
1577   i6 = 19 - i5;
1578
1579   // Charge/2, left- and righthanded couplings for in- and out-fermions.
1580   int idAbs = process[i1].idAbs();
1581   double ei = 0.5 * couplingsPtr->ef(idAbs);
1582   double li =       couplingsPtr->lf(idAbs);
1583   double ri =       couplingsPtr->rf(idAbs);
1584   idAbs     = process[i3].idAbs();
1585   double e3  = 0.5 * couplingsPtr->ef(idAbs);
1586   double l3  =       couplingsPtr->lf(idAbs);
1587   double r3  =       couplingsPtr->rf(idAbs);
1588   idAbs      = process[i5].idAbs();
1589   double e4  = 0.5 * couplingsPtr->ef(idAbs);
1590   double l4  =       couplingsPtr->lf(idAbs);
1591   double r4  =       couplingsPtr->rf(idAbs);
1592
1593   // Left- and righthanded couplings combined with propagators.
1594   c3LL = ei * ei * gamProp3 * e3 * e3
1595        + ei * li * intProp3 * e3 * l3
1596        + li * li * resProp3 * l3 * l3;
1597   c3LR = ei * ei * gamProp3 * e3 * e3
1598        + ei * li * intProp3 * e3 * r3
1599        + li * li * resProp3 * r3 * r3;
1600   c3RL = ei * ei * gamProp3 * e3 * e3
1601        + ei * ri * intProp3 * e3 * l3
1602        + ri * ri * resProp3 * l3 * l3;
1603   c3RR = ei * ei * gamProp3 * e3 * e3
1604        + ei * ri * intProp3 * e3 * r3
1605        + ri * ri * resProp3 * r3 * r3;
1606   c4LL = ei * ei * gamProp4 * e4 * e4
1607        + ei * li * intProp4 * e4 * l4
1608        + li * li * resProp4 * l4 * l4;
1609   c4LR = ei * ei * gamProp4 * e4 * e4
1610        + ei * li * intProp4 * e4 * r4
1611        + li * li * resProp4 * r4 * r4;
1612   c4RL = ei * ei * gamProp4 * e4 * e4
1613        + ei * ri * intProp4 * e4 * l4
1614        + ri * ri * resProp4 * l4 * l4;
1615   c4RR = ei * ei * gamProp4 * e4 * e4
1616        + ei * ri * intProp4 * e4 * r4
1617        + ri * ri * resProp4 * r4 * r4;
1618
1619   // Flavour weight and maximum.
1620   flavWt = (c3LL + c3LR) * (c4LL + c4LR) + (c3RL + c3RR) * (c4RL + c4RR);
1621   double flavWtMax = (c3LL + c3LR + c3RL + c3RR) * (c4LL + c4LR + c4RL + c4RR);
1622
1623   // Done.
1624   return flavWt / flavWtMax;
1625
1626 }
1627
1628 //--------------------------------------------------------------------------
1629
1630 // Evaluate weight for decay angles of the two gamma*/Z0.
1631
1632 double Sigma2ffbar2gmZgmZ::weightDecay( Event& process, int iResBeg,
1633   int iResEnd) {
1634
1635   // Two resonance decays, but with common weight.
1636   if (iResBeg != 5 || iResEnd != 6) return 1.;
1637
1638   // Set up four-products and internal products.
1639   setupProd( process, i1, i2, i3, i4, i5, i6);
1640
1641   // Flip tHat and uHat if first incoming is fermion.
1642   double tHres   = tH;
1643   double uHres   = uH;
1644   if (process[3].id() > 0) swap( tHres, uHres);
1645
1646   // Kinematics factors (norm(x) = |x|^2).
1647   double fGK135 = norm( fGK( 1, 2, 3, 4, 5, 6) / tHres
1648                       + fGK( 1, 2, 5, 6, 3, 4) / uHres );
1649   double fGK145 = norm( fGK( 1, 2, 4, 3, 5, 6) / tHres
1650                       + fGK( 1, 2, 5, 6, 4, 3) / uHres );
1651   double fGK136 = norm( fGK( 1, 2, 3, 4, 6, 5) / tHres
1652                       + fGK( 1, 2, 6, 5, 3, 4) / uHres );
1653   double fGK146 = norm( fGK( 1, 2, 4, 3, 6, 5) / tHres
1654                       + fGK( 1, 2, 6, 5, 4, 3) / uHres );
1655   double fGK253 = norm( fGK( 2, 1, 5, 6, 3, 4) / tHres
1656                       + fGK( 2, 1, 3, 4, 5, 6) / uHres );
1657   double fGK263 = norm( fGK( 2, 1, 6, 5, 3, 4) / tHres
1658                       + fGK( 2, 1, 3, 4, 6, 5) / uHres );
1659   double fGK254 = norm( fGK( 2, 1, 5, 6, 4, 3) / tHres
1660                       + fGK( 2, 1, 4, 3, 5, 6) / uHres );
1661   double fGK264 = norm( fGK( 2, 1, 6, 5, 4, 3) / tHres
1662                       + fGK( 2, 1, 4, 3, 6, 5) / uHres );
1663
1664   // Weight and maximum.
1665   double wt     = c3LL * c4LL * fGK135 + c3LR * c4LL * fGK145
1666                 + c3LL * c4LR * fGK136 + c3LR * c4LR * fGK146
1667                 + c3RL * c4RL * fGK253 + c3RR * c4RL * fGK263
1668                 + c3RL * c4RR * fGK254 + c3RR * c4RR * fGK264;
1669   double wtMax  = 16. * s3 * s4 * flavWt
1670     * ( (tHres*tHres + uHres*uHres + 2. * sH * (s3 + s4)) / (tHres * uHres)
1671       - s3 * s4 * (1. / (tHres*tHres) + 1. / (uHres*uHres)) );
1672
1673   // Done.
1674   return wt / wtMax;
1675
1676 }
1677
1678 //==========================================================================
1679
1680 // Sigma2ffbar2ZW class.
1681 // Cross section for f fbar' -> Z0 W+- (f is quark or lepton).
1682
1683 //--------------------------------------------------------------------------
1684
1685 // Initialize process.
1686
1687 void Sigma2ffbar2ZW::initProc() {
1688
1689   // Store W+- mass and width for propagator.
1690   mW   = particleDataPtr->m0(24);
1691   widW = particleDataPtr->mWidth(24);
1692   mWS  = mW*mW;
1693   mwWS = pow2(mW * widW);
1694
1695   // Left-handed couplings for up/nu- and down/e-type quarks.
1696   lun   = (hasLeptonBeams) ? couplingsPtr->lf(12) : couplingsPtr->lf(2);
1697   lde   = (hasLeptonBeams) ? couplingsPtr->lf(11) : couplingsPtr->lf(1);
1698
1699   // Common weak coupling factor.
1700   sin2thetaW = couplingsPtr->sin2thetaW();
1701   cos2thetaW = couplingsPtr->cos2thetaW();
1702   thetaWRat  = 1. / (4. * cos2thetaW);
1703   cotT       = sqrt(cos2thetaW / sin2thetaW);
1704   thetaWpt   = (9. - 8. * sin2thetaW) / 4.;
1705   thetaWmm   = (8. * sin2thetaW - 6.) / 4.;
1706
1707   // Secondary open width fractions.
1708   openFracPos = particleDataPtr->resOpenFrac(23,  24);
1709   openFracNeg = particleDataPtr->resOpenFrac(23, -24);
1710
1711 }
1712
1713 //--------------------------------------------------------------------------
1714
1715 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1716
1717 void Sigma2ffbar2ZW::sigmaKin() {
1718
1719   // Evaluate cross section, as programmed by Merlin Kole (after tidying),
1720   // based on Brown, Sahdev, Mikaelian, Phys Rev. D20 (1979) 1069.
1721   /*
1722   double resBW  = 1. / (pow2(sH - mWS) + mwWS);
1723   double prefac = 12.0 * M_PI * pow2(alpEM) / (sH2 * 8. * sin2thetaW);
1724   double temp1  = tH * uH - s3 * s4;
1725   double temp2  = temp1 / (s3 * s4);
1726   double temp3  = (s3 + s4) / sH;
1727   double temp4  = s3 * s4 / sH;
1728   double partA  = temp2 * (0.25 - 0.5 * temp3
1729     + (pow2(s3 + s4) + 8. * s3 * s4)/(4. * sH2) )
1730     + (s3 + s4)/(s3 * s4) * (sH/2. - s3 - s4 + pow2(s3 - s4)/(2. * sH));
1731   double partB1 = lun * (0.25 * temp2 * (1. - temp3 - 4. * temp4 / tH)
1732     + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / tH) );
1733   double partB2 = lde * ( 0.25 * temp2 * (1.- temp3 - 4. * temp4 / uH)
1734     + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / uH) );
1735   double partE = 0.25 * temp2 + 0.5 *(s3 + s4) / temp4;
1736   sigma0 = prefac * (pow2(cotT) * sH2 * resBW * partA
1737     + 2.* sH * cotT * resBW * (sH - mWS) * (partB2 - partB1)
1738     + pow2(lun - lde) * partE + pow2(lde) * temp1/uH2
1739     + pow2(lun) * temp1/tH2 + 2. * lun * lde * sH * (s3 + s4) / (uH * tH));
1740   */
1741
1742   // Evaluate cross section. Expression from EHLQ, with bug fix,
1743   // but can still give negative cross section so suspect.
1744   double resBW  = 1. / (pow2(sH - mWS) + mwWS);
1745   sigma0  = (M_PI / sH2) * 0.5 * pow2(alpEM / sin2thetaW);
1746   sigma0 *= sH * resBW * (thetaWpt * pT2 + thetaWmm * (s3 + s4))
1747     + (sH - mWS) * resBW * sH * (pT2 - s3 - s4) * (lun / tH - lde / uH)
1748     + thetaWRat * sH * pT2 * ( lun*lun / tH2 + lde*lde / uH2 )
1749     + 2. * thetaWRat * sH * (s3 + s4) * lun * lde / (tH * uH);
1750   // Need to protect against negative cross sections at times.
1751   sigma0 = max(0., sigma0);
1752
1753 }
1754
1755 //--------------------------------------------------------------------------
1756
1757 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1758
1759 double Sigma2ffbar2ZW::sigmaHat() {
1760
1761   // CKM and colour factors.
1762   double sigma = sigma0;
1763   if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1764
1765   // Corrections for secondary widths in Z0 and W+- decays.
1766   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
1767   sigma *= (idUp > 0) ? openFracPos : openFracNeg;
1768
1769   // Answer.
1770   return sigma;
1771
1772 }
1773
1774 //--------------------------------------------------------------------------
1775
1776 // Select identity, colour and anticolour.
1777
1778 void Sigma2ffbar2ZW::setIdColAcol() {
1779
1780   // Sign of outgoing W.
1781   int sign = 1 - 2 * (abs(id1)%2);
1782   if (id1 < 0) sign = -sign;
1783   setId( id1, id2, 23, 24 * sign);
1784
1785   // tHat is defined between (f, W-) or (fbar, W+),
1786   // so OK for u/ubar on side 1, but must swap tHat <-> uHat if d/dbar.
1787   if (abs(id1)%2 == 1) swapTU = true;
1788
1789   // Colour flow topologies. Swap when antiquarks.
1790   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1791   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1792   if (id1 < 0) swapColAcol();
1793
1794 }
1795
1796 //--------------------------------------------------------------------------
1797
1798 // Evaluate weight for Z0 and W+- decay angles.
1799
1800 double Sigma2ffbar2ZW::weightDecay( Event& process, int iResBeg, int iResEnd) {
1801
1802   // Two resonance decays, but with common weight.
1803   if (iResBeg != 5 || iResEnd != 6) return 1.;
1804
1805   // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6)
1806   // with f' fbar' from W+- and f" fbar" from Z0 (note flip Z0 <-> W+-).
1807   int i1 = (process[3].id() < 0) ? 3 : 4;
1808   int i2 = 7 - i1;
1809   int i3 = (process[9].id() > 0) ? 9 : 10;
1810   int i4 = 19 - i3;
1811   int i5 = (process[7].id() > 0) ? 7 : 8;
1812   int i6 = 15 - i5;
1813
1814   // Set up four-products and internal products.
1815   setupProd( process, i1, i2, i3, i4, i5, i6);
1816
1817   // Swap tHat and uHat if incoming fermion is downtype.
1818   double tHres   = tH;
1819   double uHres   = uH;
1820   if (process[i2].id()%2 == 1) swap( tHres, uHres);
1821
1822   //  Couplings of incoming (anti)fermions and outgoing from Z0.
1823   int idAbs     = process[i1].idAbs();
1824   double ai     = couplingsPtr->af(idAbs);
1825   double li1    = couplingsPtr->lf(idAbs);
1826   idAbs         = process[i2].idAbs();
1827   double li2    = couplingsPtr->lf(idAbs);
1828   idAbs         = process[i5].idAbs();
1829   double l4     = couplingsPtr->lf(idAbs);
1830   double r4     = couplingsPtr->rf(idAbs);
1831
1832   // W propagator/interference factor.
1833   double Wint   = cos2thetaW * (sH - mWS) / (pow2(sH - mWS) + mwWS);
1834
1835   // Combinations of couplings and kinematics (norm(x) = |x|^2).
1836   double aWZ    = li2 / tHres - 2. * Wint * ai;
1837   double bWZ    = li1 / uHres + 2. * Wint * ai;
1838   double fGK135 = norm( aWZ * fGK( 1, 2, 3, 4, 5, 6)
1839                       + bWZ * fGK( 1, 2, 5, 6, 3, 4) );
1840   double fGK136 = norm( aWZ * fGK( 1, 2, 3, 4, 6, 5)
1841                       + bWZ * fGK( 1, 2, 6, 5, 3, 4) );
1842   double xiT    = xiGK( tHres, uHres);
1843   double xiU    = xiGK( uHres, tHres);
1844   double xjTU   = xjGK( tHres, uHres);
1845
1846   // Weight and maximum weight.
1847   double wt     = l4*l4 * fGK135 + r4*r4 * fGK136;
1848   double wtMax  = 4. * s3 * s4 * (l4*l4 + r4*r4)
1849                 * (aWZ * aWZ * xiT + bWZ * bWZ * xiU + aWZ * bWZ * xjTU);
1850
1851   // Done.
1852   return wt / wtMax;
1853
1854 }
1855
1856 //==========================================================================
1857
1858 // Sigma2ffbar2WW class.
1859 // Cross section for f fbar -> W- W+ (f is quark or lepton).
1860
1861 //--------------------------------------------------------------------------
1862
1863 // Initialize process.
1864
1865 void Sigma2ffbar2WW::initProc() {
1866
1867   // Store Z0 mass and width for propagator. Common coupling factor.
1868   mZ           = particleDataPtr->m0(23);
1869   widZ         = particleDataPtr->mWidth(23);
1870   mZS          = mZ*mZ;
1871   mwZS         = pow2(mZ * widZ);
1872   thetaWRat    = 1. / (4. * couplingsPtr->sin2thetaW());
1873
1874   // Secondary open width fraction.
1875   openFracPair = particleDataPtr->resOpenFrac(24, -24);
1876
1877 }
1878
1879 //--------------------------------------------------------------------------
1880
1881 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1882
1883 void Sigma2ffbar2WW::sigmaKin() {
1884
1885   // Cross section part common for all incoming flavours.
1886   sigma0 =  (M_PI / sH2) * pow2(alpEM);
1887
1888   // Z0 propagator and gamma*/Z0 interference.
1889   double Zprop   = sH2 / (pow2(sH - mZS) + mwZS);
1890   double Zint    = Zprop * (1. - mZS / sH);
1891
1892   // Common coupling factors (g = gamma*, Z = Z0, f = t-channel fermion).
1893   cgg            = 0.5;
1894   cgZ            = thetaWRat * Zint;
1895   cZZ            = 0.5 * pow2(thetaWRat) * Zprop;
1896   cfg            = thetaWRat;
1897   cfZ            = pow2(thetaWRat) * Zint;
1898   cff            = pow2(thetaWRat);
1899
1900   // Kinematical functions.
1901   double rat34   = sH * (2. * (s3 + s4) + pT2) / (s3 * s4);
1902   double lambdaS = pow2(sH - s3 - s4) - 4. * s3 * s4;
1903   double intA    = (sH - s3 - s4) * rat34 / sH;
1904   double intB    = 4. * (s3 + s4 - pT2);
1905   gSS            = (lambdaS * rat34 + 12. * sH * pT2) / sH2;
1906   gTT            = rat34 + 4. * sH * pT2 / tH2;
1907   gST            = intA + intB / tH;
1908   gUU            = rat34 + 4. * sH * pT2 / uH2;
1909   gSU            = intA + intB / uH;
1910
1911 }
1912
1913 //--------------------------------------------------------------------------
1914
1915 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1916
1917 double Sigma2ffbar2WW::sigmaHat() {
1918
1919   // Flavour-specific couplings.
1920   int idAbs = abs(id1);
1921   double ei = couplingsPtr->ef(idAbs);
1922   double vi = couplingsPtr->vf(idAbs);
1923   double ai = couplingsPtr->af(idAbs);
1924
1925   // Combine, with different cases for up- and down-type in-flavours.
1926   double sigma = sigma0;
1927   sigma *= (idAbs%2 == 1)
1928     ? (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
1929       + (cfg * ei + cfZ * (vi + ai)) * gST + cff * gTT
1930     : (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
1931       - (cfg * ei + cfZ * (vi + ai)) * gSU + cff * gUU;
1932
1933   // Initial-state colour factor. Correction for secondary widths. Answer.
1934   if (idAbs < 9) sigma /= 3.;
1935   sigma *= openFracPair;
1936   return sigma;
1937
1938 }
1939
1940 //--------------------------------------------------------------------------
1941
1942 // Select identity, colour and anticolour.
1943
1944 void Sigma2ffbar2WW::setIdColAcol() {
1945
1946   // Always order W- W+, i.e. W- first.
1947   setId( id1, id2, -24, 24);
1948
1949   // tHat is defined between (f, W-) or (fbar, W+),
1950   if (id1 < 0) swapTU = true;
1951
1952   // Colour flow topologies. Swap when antiquarks.
1953   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1954   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1955   if (id1 < 0) swapColAcol();
1956
1957 }
1958
1959 //--------------------------------------------------------------------------
1960
1961 // Evaluate weight for W+ and W- decay angles.
1962
1963 double Sigma2ffbar2WW::weightDecay( Event& process, int iResBeg, int iResEnd) {
1964
1965   // Two resonance decays, but with common weight.
1966   if (iResBeg != 5 || iResEnd != 6) return 1.;
1967
1968   // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1969   // with f' fbar' from W- and f" fbar" from W+.
1970   int i1 = (process[3].id() < 0) ? 3 : 4;
1971   int i2 = 7 - i1;
1972   int i3 = (process[7].id() > 0) ? 7 : 8;
1973   int i4 = 15 - i3;
1974   int i5 = (process[9].id() > 0) ? 9 : 10;
1975   int i6 = 19 - i5;
1976
1977   // Set up four-products and internal products.
1978   setupProd( process, i1, i2, i3, i4, i5, i6);
1979
1980   // tHat and uHat of fbar f -> W- W+ opposite to previous convention.
1981   double tHres   = uH;
1982   double uHres   = tH;
1983
1984   //  Couplings of incoming (anti)fermion.
1985   int idAbs     = process[i1].idAbs();
1986   double ai     = couplingsPtr->af(idAbs);
1987   double li     = couplingsPtr->lf(idAbs);
1988   double ri     = couplingsPtr->rf(idAbs);
1989
1990   // gamma*/Z0 propagator/interference factor.
1991   double Zint   = mZS * (sH - mZS) / (pow2(sH - mZS) + mwZS);
1992
1993   // Combinations of couplings and kinematics (norm(x) = |x|^2).
1994   double dWW    = (li * Zint + ai) / sH;
1995   double aWW    = dWW + 0.5 * (ai + 1.) / tHres;
1996   double bWW    = dWW + 0.5 * (ai - 1.) / uHres;
1997   double cWW    = ri * Zint / sH;
1998   double fGK135 = norm( aWW * fGK( 1, 2, 3, 4, 5, 6)
1999                       - bWW * fGK( 1, 2, 5, 6, 3, 4) );
2000   double fGK253 = norm( cWW * ( fGK( 2, 1, 5, 6, 3, 4)
2001                               - fGK( 2, 1, 3, 4, 5, 6) ) );
2002   double xiT    = xiGK( tHres, uHres);
2003   double xiU    = xiGK( uHres, tHres);
2004   double xjTU   = xjGK( tHres, uHres);
2005
2006   // Weight and maximum weight.
2007   double wt     = fGK135 + fGK253;
2008   double wtMax  = 4. * s3 * s4
2009                 * ( aWW * aWW * xiT + bWW * bWW * xiU - aWW * bWW * xjTU
2010                   + cWW * cWW * (xiT + xiU - xjTU) );
2011
2012   // Done.
2013   return wt / wtMax;
2014 }
2015
2016 //==========================================================================
2017
2018 // Sigma2ffbargmZggm class.
2019 // Collects common methods for f fbar -> gamma*/Z0 g/gamma and permutations.
2020
2021 //--------------------------------------------------------------------------
2022
2023 // Initialize process.
2024
2025 void Sigma2ffbargmZggm::initProc() {
2026
2027   // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
2028   gmZmode     = settingsPtr->mode("WeakZ0:gmZmode");
2029
2030   // Store Z0 mass and width for propagator.
2031   mRes        = particleDataPtr->m0(23);
2032   GammaRes    = particleDataPtr->mWidth(23);
2033   m2Res       = mRes*mRes;
2034   GamMRat     = GammaRes / mRes;
2035   thetaWRat   = 1. / (16. * couplingsPtr->sin2thetaW()
2036                 * couplingsPtr->cos2thetaW());
2037
2038   // Set pointer to particle properties and decay table.
2039   particlePtr = particleDataPtr->particleDataEntryPtr(23);
2040
2041 }
2042
2043 //--------------------------------------------------------------------------
2044
2045 // Evaluate sum of flavour couplings times phase space.
2046
2047 void Sigma2ffbargmZggm::flavSum() {
2048
2049   // Coupling factors for Z0 subsystem.
2050   double alpSZ = couplingsPtr->alphaS(s3);
2051   double colQZ = 3. * (1. + alpSZ / M_PI);
2052
2053   // Reset quantities to sum. Declare variables in loop.
2054   gamSum = 0.;
2055   intSum = 0.;
2056   resSum = 0.;
2057   int    onMode;
2058   double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
2059
2060   // Loop over all Z0 decay channels.
2061   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
2062     int idAbs = abs( particlePtr->channel(i).product(0) );
2063
2064     // Only contributions from three fermion generations, except top.
2065     if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
2066       mf = particleDataPtr->m0(idAbs);
2067
2068       // Check that above threshold. Phase space.
2069       if (m3 > 2. * mf + MASSMARGIN) {
2070         mr    = pow2(mf / m3);
2071         betaf = sqrtpos(1. - 4. * mr);
2072         psvec = betaf * (1. + 2. * mr);
2073         psaxi = pow3(betaf);
2074
2075         // Combine phase space with couplings.
2076         ef2    = couplingsPtr->ef2(idAbs) * psvec;
2077         efvf   = couplingsPtr->efvf(idAbs) * psvec;
2078         vf2af2 = couplingsPtr->vf2(idAbs) * psvec
2079                + couplingsPtr->af2(idAbs) * psaxi;
2080         colf   = (idAbs < 6) ? colQZ : 1.;
2081
2082         // Store sum of combinations. For outstate only open channels.
2083         onMode = particlePtr->channel(i).onMode();
2084         if (onMode == 1 || onMode == 2) {
2085           gamSum += colf * ef2;
2086           intSum += colf * efvf;
2087           resSum += colf * vf2af2;
2088         }
2089
2090       // End loop over fermions.
2091       }
2092     }
2093   }
2094
2095   // Done. Return values in gamSum, intSum and resSum.
2096
2097 }
2098
2099 //--------------------------------------------------------------------------
2100
2101 // Calculate common parts of gamma/interference/Z0 propagator terms.
2102
2103 void Sigma2ffbargmZggm::propTerm() {
2104
2105   // Calculate prefactors for gamma/interference/Z0 cross section terms.
2106   gamProp = 4. * alpEM / (3. * M_PI * s3);
2107   intProp = gamProp * 2. * thetaWRat * s3 * (s3 - m2Res)
2108           / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2109   resProp = gamProp * pow2(thetaWRat * s3)
2110           / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2111
2112   // Optionally only keep gamma* or Z0 term.
2113   if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
2114   if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
2115
2116 }
2117
2118 //--------------------------------------------------------------------------
2119
2120 // Evaluate weight for gamma*/Z0 decay angle.
2121
2122 double Sigma2ffbargmZggm::weightDecay( Event& process, int iResBeg,
2123   int iResEnd) {
2124
2125   // Z should sit in entry 5 and one more parton in entry 6.
2126   if (iResBeg != 5 || iResEnd != 6) return 1.;
2127
2128   // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2129   // where f' fbar' come from gamma*/Z0 decay.
2130   int i1, i2;
2131   int i3 = (process[7].id() > 0) ? 7 : 8;
2132   int i4 = 15 - i3;
2133
2134   // Order so that fbar(1) f(2) -> gamma*/Z0 g/gamma.
2135   if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2136     i1 = (process[3].id() < 0) ? 3 : 4;
2137     i2 = 7 - i1;
2138
2139   // Order so that f(2)/fbar(1)  g/gamma -> f(1)/fbar(2) f'(3) gamma*/Z0.
2140   } else if (process[3].idAbs() < 20) {
2141     i1 = (process[3].id() < 0) ? 3 : 6;
2142     i2 = 9 - i1;
2143   } else {
2144     i1 = (process[4].id() < 0) ? 4 : 6;
2145     i2 = 10 - i1;
2146   }
2147
2148   // Charge/2, left- and righthanded couplings for in- and out-fermion.
2149   int id1Abs   = process[i1].idAbs();
2150   double ei    = 0.5 * couplingsPtr->ef(id1Abs);
2151   double li    =       couplingsPtr->lf(id1Abs);
2152   double ri    =       couplingsPtr->rf(id1Abs);
2153   int id3Abs   = process[i3].idAbs();
2154   double ef    = 0.5 * couplingsPtr->ef(id3Abs);
2155   double lf    =       couplingsPtr->lf(id3Abs);
2156   double rf    =       couplingsPtr->rf(id3Abs);
2157
2158   // Combinations of left/right for in/out, gamma*/interference/Z0.
2159   double clilf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*lf
2160                + li*li * resProp * lf*lf;
2161   double clirf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*rf
2162                + li*li * resProp * rf*rf;
2163   double crilf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*lf
2164                + ri*ri * resProp * lf*lf;
2165   double crirf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*rf
2166                + ri*ri * resProp * rf*rf;
2167
2168   // Evaluate four-vector products.
2169   double p13   = process[i1].p() * process[i3].p();
2170   double p14   = process[i1].p() * process[i4].p();
2171   double p23   = process[i2].p() * process[i3].p();
2172   double p24   = process[i2].p() * process[i4].p();
2173
2174   // Calculate weight and its maximum.
2175   double wt    = (clilf + crirf) * (p13*p13 + p24*p24)
2176                + (clirf + crilf) * (p14*p14 + p23*p23) ;
2177   double wtMax = (clilf + clirf + crilf + crirf)
2178                * (pow2(p13 + p14) + pow2(p23 + p24));
2179
2180   // Done.
2181   return (wt / wtMax);
2182
2183 }
2184
2185 //==========================================================================
2186
2187 // Sigma2qqbar2gmZg class.
2188 // Cross section for q qbar -> gamma*/Z0 g.
2189
2190 //--------------------------------------------------------------------------
2191
2192 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2193
2194 void Sigma2qqbar2gmZg::sigmaKin() {
2195
2196   // Cross section part common for all incoming flavours.
2197   sigma0 = (M_PI / sH2) * (alpEM * alpS)
2198     * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2199
2200   // Calculate flavour sums for final state.
2201   flavSum();
2202
2203   // Calculate prefactors for gamma/interference/Z0 cross section terms.
2204   propTerm();
2205
2206 }
2207
2208 //--------------------------------------------------------------------------
2209
2210 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2211
2212 double Sigma2qqbar2gmZg::sigmaHat() {
2213
2214   // Combine gamma, interference and Z0 parts.
2215   int idAbs    = abs(id1);
2216   double sigma = sigma0
2217                * ( couplingsPtr->ef2(idAbs)    * gamProp * gamSum
2218                  + couplingsPtr->efvf(idAbs)   * intProp * intSum
2219                  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2220
2221   // Correct for the running-width Z0 propagater weight in PhaseSpace.
2222   sigma       /= runBW3;
2223
2224   // Answer.
2225   return sigma;
2226
2227 }
2228
2229 //--------------------------------------------------------------------------
2230
2231 // Select identity, colour and anticolour.
2232
2233 void Sigma2qqbar2gmZg::setIdColAcol() {
2234
2235   // Flavours trivial.
2236   setId( id1, id2, 23, 21);
2237
2238   // Colour flow topologies. Swap when antiquarks.
2239   setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2240   if (id1 < 0) swapColAcol();
2241
2242 }
2243
2244 //==========================================================================
2245
2246 // Sigma2qg2gmZq class.
2247 // Cross section for q g -> gamma*/Z0 q.
2248
2249 //--------------------------------------------------------------------------
2250
2251 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2252
2253 void Sigma2qg2gmZq::sigmaKin() {
2254
2255   // Cross section part common for all incoming flavours.
2256   sigma0 = (M_PI / sH2) * (alpEM * alpS)
2257     * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2258
2259   // Calculate flavour sums for final state.
2260   flavSum();
2261
2262   // Calculate prefactors for gamma/interference/Z0 cross section terms.
2263   propTerm();
2264
2265 }
2266
2267 //--------------------------------------------------------------------------
2268
2269 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2270
2271 double Sigma2qg2gmZq::sigmaHat() {
2272
2273   // Combine gamma, interference and Z0 parts.
2274   int idAbs    = (id2 == 21) ? abs(id1) : abs(id2);
2275   double sigma = sigma0
2276                * ( couplingsPtr->ef2(idAbs)    * gamProp * gamSum
2277                  + couplingsPtr->efvf(idAbs)   * intProp * intSum
2278                  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2279
2280   // Correct for the running-width Z0 propagater weight in PhaseSpace.
2281   sigma       /= runBW3;
2282
2283   // Answer.
2284   return sigma;
2285
2286 }
2287
2288 //--------------------------------------------------------------------------
2289
2290 // Select identity, colour and anticolour.
2291
2292 void Sigma2qg2gmZq::setIdColAcol() {
2293
2294   // Flavour set up for q g -> gamma*/Z0 q.
2295   int idq = (id2 == 21) ? id1 : id2;
2296   setId( id1, id2, 23, idq);
2297
2298   // tH defined between f and f': must swap tHat <-> uHat if q g in.
2299   swapTU = (id2 == 21);
2300
2301   // Colour flow topologies. Swap when antiquarks.
2302   if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2303   else           setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2304   if (idq < 0) swapColAcol();
2305
2306 }
2307
2308 //==========================================================================
2309
2310 // Sigma2ffbar2gmZgm class.
2311 // Cross section for f fbar -> gamma*/Z0 gamma.
2312
2313 //--------------------------------------------------------------------------
2314
2315 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2316
2317 void Sigma2ffbar2gmZgm::sigmaKin() {
2318
2319   // Cross section part common for all incoming flavours.
2320   sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2321     * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2322
2323   // Calculate flavour sums for final state.
2324   flavSum();
2325
2326   // Calculate prefactors for gamma/interference/Z0 cross section terms.
2327   propTerm();
2328
2329
2330 }
2331
2332 //--------------------------------------------------------------------------
2333
2334 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2335
2336 double Sigma2ffbar2gmZgm::sigmaHat() {
2337
2338   // Combine gamma, interference and Z0 parts.
2339   int idAbs    = abs(id1);
2340   double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2341                * ( couplingsPtr->ef2(idAbs)    * gamProp * gamSum
2342                  + couplingsPtr->efvf(idAbs)   * intProp * intSum
2343                  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2344
2345   // Correct for the running-width Z0 propagater weight in PhaseSpace.
2346   sigma       /= runBW3;
2347
2348   // Colour factor. Answer.
2349   if (idAbs < 9) sigma /= 3.;
2350   return sigma;
2351
2352 }
2353
2354 //--------------------------------------------------------------------------
2355
2356 // Select identity, colour and anticolour.
2357
2358 void Sigma2ffbar2gmZgm::setIdColAcol() {
2359
2360   // Flavours trivial.
2361   setId( id1, id2, 23, 22);
2362
2363   // Colour flow topologies. Swap when antiquarks.
2364   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2365   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2366   if (id1 < 0) swapColAcol();
2367
2368 }
2369
2370 //==========================================================================
2371
2372 // Sigma2fgm2gmZf class.
2373 // Cross section for f gamma -> gamma*/Z0 f'.
2374
2375 //--------------------------------------------------------------------------
2376
2377 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2378
2379 void Sigma2fgm2gmZf::sigmaKin() {
2380
2381   // Cross section part common for all incoming flavours.
2382   sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2383     * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (- sH * uH);
2384
2385   // Calculate flavour sums for final state.
2386   flavSum();
2387
2388   // Calculate prefactors for gamma/interference/Z0 cross section terms.
2389   propTerm();
2390
2391 }
2392
2393 //--------------------------------------------------------------------------
2394
2395 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2396
2397 double Sigma2fgm2gmZf::sigmaHat() {
2398
2399   // Combine gamma, interference and Z0 parts.
2400   int idAbs    = (id2 == 22) ? abs(id1) : abs(id2);
2401   double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2402                * ( couplingsPtr->ef2(idAbs)    * gamProp * gamSum
2403                  + couplingsPtr->efvf(idAbs)   * intProp * intSum
2404                  + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2405
2406   // Correct for the running-width Z0 propagater weight in PhaseSpace.
2407   sigma         /= runBW3;
2408
2409   // Answer.
2410   return sigma;
2411
2412 }
2413
2414 //--------------------------------------------------------------------------
2415
2416 // Select identity, colour and anticolour.
2417
2418 void Sigma2fgm2gmZf::setIdColAcol() {
2419
2420   // Flavour set up for q gamma -> gamma*/Z0 q.
2421   int idq = (id2 == 22) ? id1 : id2;
2422   setId( id1, id2, 23, idq);
2423
2424   // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2425   swapTU = (id2 == 22);
2426
2427   // Colour flow topologies. Swap when antiquarks.
2428   if      (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2429   else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2430   else                   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2431   if (idq < 0) swapColAcol();
2432
2433 }
2434
2435 //==========================================================================
2436
2437 // Sigma2ffbarWggm class.
2438 // Collects common methods for f fbar -> W+- g/gamma and permutations.
2439
2440 //--------------------------------------------------------------------------
2441
2442 // Evaluate weight for W+- decay angle.
2443
2444 double Sigma2ffbarWggm::weightDecay( Event& process, int iResBeg,
2445   int iResEnd) {
2446
2447   // W should sit in entry 5 and one more parton in entry 6.
2448   if (iResBeg != 5 || iResEnd != 6) return 1.;
2449
2450   // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2451   // where f' fbar' come from W+- decay.
2452   int i1, i2;
2453   int i3 = (process[7].id() > 0) ? 7 : 8;
2454   int i4 = 15 - i3;
2455
2456   // Order so that fbar(1) f(2) -> W+- g/gamma.
2457   if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2458     i1 = (process[3].id() < 0) ? 3 : 4;
2459     i2 = 7 - i1;
2460
2461   // Order so that f(2)/fbar(1)  g/gamma -> f(1)/fbar(2) f'(3) W+-.
2462   } else if (process[3].idAbs() < 20) {
2463     i1 = (process[3].id() < 0) ? 3 : 6;
2464     i2 = 9 - i1;
2465   } else {
2466     i1 = (process[4].id() < 0) ? 4 : 6;
2467     i2 = 10 - i1;
2468   }
2469
2470   // Evaluate four-vector products.
2471   double p13 = process[i1].p() * process[i3].p();
2472   double p14 = process[i1].p() * process[i4].p();
2473   double p23 = process[i2].p() * process[i3].p();
2474   double p24 = process[i2].p() * process[i4].p();
2475
2476   // Calculate weight and its maximum.
2477   double wt    = pow2(p13) + pow2(p24);
2478   double wtMax = pow2(p13 + p14) + pow2(p23 + p24);
2479
2480   // Done.
2481   return (wt / wtMax);
2482
2483 }
2484
2485 //==========================================================================
2486
2487 // Sigma2qqbar2Wg class.
2488 // Cross section for q qbar' -> W+- g.
2489
2490 //--------------------------------------------------------------------------
2491
2492 // Initialize process.
2493
2494 void Sigma2qqbar2Wg::initProc() {
2495
2496   // Secondary open width fractions, relevant for top (or heavier).
2497   openFracPos = particleDataPtr->resOpenFrac(24);
2498   openFracNeg = particleDataPtr->resOpenFrac(-24);
2499
2500 }
2501
2502 //--------------------------------------------------------------------------
2503
2504 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2505
2506 void Sigma2qqbar2Wg::sigmaKin() {
2507
2508   // Cross section part common for all incoming flavours.
2509   sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2510     * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2511
2512 }
2513
2514 //--------------------------------------------------------------------------
2515
2516 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2517
2518 double Sigma2qqbar2Wg::sigmaHat() {
2519
2520   // CKM factor. Secondary width for W+ or W-.
2521   double sigma = sigma0 * couplingsPtr->V2CKMid(abs(id1), abs(id2));
2522   int idUp     = (abs(id1)%2 == 0) ? id1 : id2;
2523   sigma       *= (idUp > 0) ? openFracPos : openFracNeg;
2524
2525   // Answer.
2526   return sigma;
2527
2528 }
2529
2530 //--------------------------------------------------------------------------
2531
2532 // Select identity, colour and anticolour.
2533
2534 void Sigma2qqbar2Wg::setIdColAcol() {
2535
2536   // Sign of outgoing W.
2537   int sign = 1 - 2 * (abs(id1)%2);
2538   if (id1 < 0) sign = -sign;
2539   setId( id1, id2, 24 * sign, 21);
2540
2541   // Colour flow topologies. Swap when antiquarks.
2542   setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2543   if (id1 < 0) swapColAcol();
2544
2545 }
2546
2547 //==========================================================================
2548
2549 // Sigma2qg2Wq class.
2550 // Cross section for q g -> W+- q'.
2551
2552 //--------------------------------------------------------------------------
2553
2554 // Initialize process.
2555
2556 void Sigma2qg2Wq::initProc() {
2557
2558   // Secondary open width fractions, relevant for top (or heavier).
2559   openFracPos = particleDataPtr->resOpenFrac(24);
2560   openFracNeg = particleDataPtr->resOpenFrac(-24);
2561
2562 }
2563
2564 //--------------------------------------------------------------------------
2565
2566 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2567
2568 void Sigma2qg2Wq::sigmaKin() {
2569
2570   // Cross section part common for all incoming flavours.
2571   sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2572     * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2573
2574 }
2575
2576 //--------------------------------------------------------------------------
2577
2578 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2579
2580 double Sigma2qg2Wq::sigmaHat() {
2581
2582   // CKM factor. Secondary width for W+ or W-.
2583   int idAbs    = (id2 == 21) ? abs(id1) : abs(id2);
2584   double sigma = sigma0 * couplingsPtr->V2CKMsum(idAbs);
2585   int idUp     = (id2 == 21) ? id1 : id2;
2586   if (idAbs%2 == 1) idUp = -idUp;
2587   sigma       *= (idUp > 0) ? openFracPos : openFracNeg;
2588
2589   // Answer.
2590   return sigma;
2591
2592 }
2593
2594 //--------------------------------------------------------------------------
2595
2596 // Select identity, colour and anticolour.
2597
2598 void Sigma2qg2Wq::setIdColAcol() {
2599
2600   // Sign of outgoing W. Flavour of outgoing quark.
2601   int idq           = (id2 == 21) ? id1 : id2;
2602   int sign          = 1 - 2 * (abs(idq)%2);
2603   if (idq < 0) sign = -sign;
2604   id4 = couplingsPtr->V2CKMpick(idq);
2605
2606   // Flavour set up for q g -> W q.
2607   setId( id1, id2, 24 * sign, id4);
2608
2609   // tH defined between f and f': must swap tHat <-> uHat if q g in.
2610   swapTU = (id2 == 21);
2611
2612   // Colour flow topologies. Swap when antiquarks.
2613   if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2614   else           setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2615   if (idq < 0) swapColAcol();
2616
2617 }
2618
2619 //==========================================================================
2620
2621 // Sigma2ffbar2Wgm class.
2622 // Cross section for f fbar' -> W+- gamma.
2623
2624 //--------------------------------------------------------------------------
2625
2626 // Initialize process.
2627
2628 void Sigma2ffbar2Wgm::initProc() {
2629
2630   // Secondary open width fractions, relevant for top (or heavier).
2631   openFracPos = particleDataPtr->resOpenFrac(24);
2632   openFracNeg = particleDataPtr->resOpenFrac(-24);
2633
2634 }
2635
2636 //--------------------------------------------------------------------------
2637
2638 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2639
2640 void Sigma2ffbar2Wgm::sigmaKin() {
2641
2642   // Cross section part common for all incoming flavours.
2643   sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2644     * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2645 }
2646
2647 //--------------------------------------------------------------------------
2648
2649 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2650
2651 double Sigma2ffbar2Wgm::sigmaHat() {
2652
2653   // Extrafactor different for e nu and q qbar' instate.
2654   int id1Abs = abs(id1);
2655   int id2Abs = abs(id2);
2656   double chgUp = (id1Abs > 10) ? 0. : 2./3.;
2657   double sigma = sigma0 * pow2( chgUp - tH / (tH + uH) );
2658
2659   // CKM and colour factors. Secondary width for W+ or W-.
2660   if (id1Abs < 9) sigma *= couplingsPtr->V2CKMid(id1Abs, id2Abs) / 3.;
2661   int idUp     = (abs(id1)%2 == 0) ? id1 : id2;
2662   sigma       *= (idUp > 0) ? openFracPos : openFracNeg;
2663
2664   // Answer.
2665   return sigma;
2666
2667 }
2668
2669 //--------------------------------------------------------------------------
2670
2671 // Select identity, colour and anticolour.
2672
2673 void Sigma2ffbar2Wgm::setIdColAcol() {
2674
2675   // Sign of outgoing W.
2676   int sign          = 1 - 2 * (abs(id1)%2);
2677   if (id1 < 0) sign = -sign;
2678   setId( id1, id2, 24 * sign, 22);
2679
2680   // tH defined between (f,W-) or (fbar',W+).
2681   swapTU = (sign * id1 > 0);
2682
2683   // Colour flow topologies. Swap when antiquarks.
2684   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2685   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2686   if (id1 < 0) swapColAcol();
2687
2688 }
2689
2690 //==========================================================================
2691
2692 // Sigma2fgm2Wf class.
2693 // Cross section for f gamma -> W+- f'.
2694
2695 //--------------------------------------------------------------------------
2696
2697 // Initialize process.
2698
2699 void Sigma2fgm2Wf::initProc() {
2700
2701   // Secondary open width fractions, relevant for top (or heavier).
2702   openFracPos = particleDataPtr->resOpenFrac(24);
2703   openFracNeg = particleDataPtr->resOpenFrac(-24);
2704
2705 }
2706
2707 //--------------------------------------------------------------------------
2708
2709 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2710
2711 void Sigma2fgm2Wf::sigmaKin() {
2712
2713   // Cross section part common for all incoming flavours.
2714   sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2715     * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (pT2 * s3 - sH * uH);
2716
2717 }
2718
2719 //--------------------------------------------------------------------------
2720
2721 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2722
2723 double Sigma2fgm2Wf::sigmaHat() {
2724
2725   // Extrafactor dependent on charge of incoming fermion.
2726   int idAbs     = (id2 == 22) ? abs(id1) : abs(id2);
2727   double charge = (idAbs > 10) ? 1. : ( (idAbs%2 == 1) ? 1./3. : 2./3. );
2728   double sigma  = sigma0 * pow2( charge  - sH / (sH + uH) );
2729
2730   // CKM factor. Secondary width for W+ or W-.
2731   sigma        *= couplingsPtr->V2CKMsum(idAbs);
2732   int idUp      = (id2 == 22) ? id1 : id2;
2733   if (idAbs%2 == 1) idUp = -idUp;
2734   sigma        *= (idUp > 0) ? openFracPos : openFracNeg;
2735
2736   // Answer.
2737   return sigma;
2738
2739 }
2740
2741 //--------------------------------------------------------------------------
2742
2743 // Select identity, colour and anticolour.
2744
2745 void Sigma2fgm2Wf::setIdColAcol() {
2746
2747   // Sign of outgoing W. Flavour of outgoing fermion.
2748   int idq           = (id2 == 22) ? id1 : id2;
2749   int sign          = 1 - 2 * (abs(idq)%2);
2750   if (idq < 0) sign = -sign;
2751   id4 = couplingsPtr->V2CKMpick(idq);
2752
2753   // Flavour set up for q gamma -> W q.
2754   setId( id1, id2, 24 * sign, id4);
2755
2756   // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2757   swapTU = (id2 == 22);
2758
2759   // Colour flow topologies. Swap when antiquarks.
2760   if      (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2761   else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2762   else                   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2763   if (idq < 0) swapColAcol();
2764
2765 }
2766
2767 //==========================================================================
2768
2769 // Sigma2gmgm2ffbar class.
2770 // Cross section for gamma gamma -> l lbar.
2771
2772 //--------------------------------------------------------------------------
2773
2774 // Initialize process wrt flavour.
2775
2776 void Sigma2gmgm2ffbar::initProc() {
2777
2778   // Process name.
2779   nameSave = "gamma gamma -> f fbar";
2780   if (idNew ==  1) nameSave = "gamma gamma -> q qbar (uds)";
2781   if (idNew ==  4) nameSave = "gamma gamma -> c cbar";
2782   if (idNew ==  5) nameSave = "gamma gamma -> b bbar";
2783   if (idNew ==  6) nameSave = "gamma gamma -> t tbar";
2784   if (idNew == 11) nameSave = "gamma gamma -> e+ e-";
2785   if (idNew == 13) nameSave = "gamma gamma -> mu+ mu-";
2786   if (idNew == 15) nameSave = "gamma gamma -> tau+ tau-";
2787
2788   // Generate massive phase space, except for u+d+s.
2789   idMass = 0;
2790   if (idNew > 3) idMass = idNew;
2791
2792   // Charge and colour factor.
2793   ef4 = 1.;
2794   if (idNew == 1) ef4 = 3. * (pow4(2./3.) + 2. * pow4(1./3.));
2795   if (idNew == 4 || idNew == 6) ef4 = 3. * pow4(2./3.);
2796   if (idNew == 5) ef4 = 3. * pow4(1./3.);
2797
2798   // Secondary open width fraction.
2799   openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
2800
2801 }
2802
2803 //--------------------------------------------------------------------------
2804
2805 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2806
2807 void Sigma2gmgm2ffbar::sigmaKin() {
2808
2809   // Pick current flavour for u+d+s mix by e_q^4 weights.
2810   if (idNew == 1) {
2811     double rId = 18. * rndmPtr->flat();
2812     idNow = 1;
2813     if (rId > 1.)  idNow = 2;
2814     if (rId > 17.) idNow = 3;
2815     s34Avg = pow2(particleDataPtr->m0(idNow));
2816   } else {
2817     idNow = idNew;
2818     s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
2819   }
2820
2821   // Modified Mandelstam variables for massive kinematics with m3 = m4.
2822   double tHQ    = -0.5 * (sH - tH + uH);
2823   double uHQ    = -0.5 * (sH + tH - uH);
2824   double tHQ2   = tHQ * tHQ;
2825   double uHQ2   = uHQ * uHQ;
2826
2827   // Calculate kinematics dependence.
2828   if (sH < 4. * s34Avg) sigTU = 0.;
2829   else sigTU = 2. * (tHQ * uHQ - s34Avg * sH)
2830     * (tHQ2 + uHQ2 + 2. * s34Avg * sH) / (tHQ2 * uHQ2);
2831
2832   // Answer.
2833   sigma = (M_PI / sH2) * pow2(alpEM) * ef4 * sigTU * openFracPair;
2834
2835 }
2836
2837 //--------------------------------------------------------------------------
2838
2839 // Select identity, colour and anticolour.
2840
2841 void Sigma2gmgm2ffbar::setIdColAcol() {
2842
2843   // Flavours are trivial.
2844   setId( id1, id2, idNow, -idNow);
2845
2846   // Colour flow in singlet state.
2847   if (idNow < 10) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
2848   else            setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2849
2850 }
2851
2852 //==========================================================================
2853
2854 } // end namespace Pythia8