PYTHIA8/pythia8140/src/SigmaQCD.cxx

   1 // SigmaQCD.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 // QCD simulation classes.
   8
   9 #include "SigmaQCD.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // Sigma0AB2AB class.
  16 // Cross section for elastic scattering A B -> A B.
  17
  18 //--------------------------------------------------------------------------
  19
  20 // Select identity, colour and anticolour.
  21
  22 void Sigma0AB2AB::setIdColAcol() {
  23
  24   // Flavours and colours are trivial.
  25   setId( idA, idB, idA, idB);
  26   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
  27 }
  28
  29 //==========================================================================
  30
  31 // Sigma0AB2XB class.
  32 // Cross section for single diffractive scattering A B -> X B.
  33
  34 //--------------------------------------------------------------------------
  35
  36 // Select identity, colour and anticolour.
  37
  38 void Sigma0AB2XB::setIdColAcol() {
  39
  40   // Flavours and colours are trivial.
  41   int idX          = 10* (abs(idA) / 10) + 9900000;
  42   if (idA < 0) idX = -idX;
  43   setId( idA, idB, idX, idB);
  44   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
  45
  46 }
  47
  48 //==========================================================================
  49
  50 // Sigma0AB2AX class.
  51 // Cross section for single diffractive scattering A B -> A X.
  52
  53 //--------------------------------------------------------------------------
  54
  55 // Select identity, colour and anticolour.
  56
  57 void Sigma0AB2AX::setIdColAcol() {
  58
  59   // Flavours and colours are trivial.
  60   int idX          = 10* (abs(idB) / 10) + 9900000;
  61   if (idB < 0) idX = -idX;
  62   setId( idA, idB, idA, idX);
  63   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
  64
  65 }
  66
  67 //==========================================================================
  68
  69 // Sigma0AB2XX class.
  70 // Cross section for double diffractive scattering A B -> X X.
  71
  72 //--------------------------------------------------------------------------
  73
  74 // Select identity, colour and anticolour.
  75
  76 void Sigma0AB2XX::setIdColAcol() {
  77
  78   // Flavours and colours are trivial.
  79   int          idX1 = 10* (abs(idA) / 10) + 9900000;
  80   if (idA < 0) idX1 = -idX1;
  81   int          idX2 = 10* (abs(idB) / 10) + 9900000;
  82   if (idB < 0) idX2 = -idX2;
  83   setId( idA, idB, idX1, idX2);
  84   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
  85
  86 }
  87
  88 //==========================================================================
  89
  90 // Sigma2gg2gg class.
  91 // Cross section for g g -> g g.
  92
  93 //--------------------------------------------------------------------------
  94
  95 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
  96
  97 void Sigma2gg2gg::sigmaKin() {
  98
  99   // Calculate kinematics dependence.
 100   sigTS  = (9./4.) * (tH2 / sH2 + 2. * tH / sH + 3. + 2. * sH / tH
 101            + sH2 / tH2);
 102   sigUS  = (9./4.) * (uH2 / sH2 + 2. * uH / sH + 3. + 2. * sH / uH
 103            + sH2 / uH2);
 104   sigTU  = (9./4.) * (tH2 / uH2 + 2. * tH / uH + 3. + 2. * uH / tH
 105            + uH2 / tH2);
 106   sigSum = sigTS + sigUS + sigTU;
 107
 108   // Answer contains factor 1/2 from identical gluons.
 109   sigma  = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
 110
 111 }
 112
 113 //--------------------------------------------------------------------------
 114
 115 // Select identity, colour and anticolour.
 116
 117 void Sigma2gg2gg::setIdColAcol() {
 118
 119   // Flavours are trivial.
 120   setId( id1, id2, 21, 21);
 121
 122   // Three colour flow topologies, each with two orientations.
 123   double sigRand = sigSum * rndmPtr->flat();
 124   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
 125   else if (sigRand < sigTS + sigUS)
 126                        setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
 127   else                 setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
 128   if (rndmPtr->flat() > 0.5) swapColAcol();
 129
 130 }
 131
 132 //==========================================================================
 133
 134 // Sigma2gg2qqbar class.
 135 // Cross section for g g -> q qbar (q = u, d, s, i.e. almost massless).
 136
 137 //--------------------------------------------------------------------------
 138
 139 // Initialize process.
 140
 141 void Sigma2gg2qqbar::initProc() {
 142
 143   // Read number of quarks to be considered in massless approximation.
 144   nQuarkNew       = settingsPtr->mode("HardQCD:nQuarkNew");
 145
 146 }
 147
 148 //--------------------------------------------------------------------------
 149
 150 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 151
 152 void Sigma2gg2qqbar::sigmaKin() {
 153
 154   // Pick new flavour.
 155   idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
 156   mNew  = particleDataPtr->m0(idNew);
 157   m2New = mNew*mNew;
 158
 159   // Calculate kinematics dependence.
 160   sigTS = 0.;
 161   sigUS = 0.;
 162   if (sH > 4. * m2New) {
 163     sigTS = (1./6.) * uH / tH - (3./8.) * uH2 / sH2;
 164     sigUS = (1./6.) * tH / uH - (3./8.) * tH2 / sH2;
 165   }
 166   sigSum = sigTS + sigUS;
 167
 168   // Answer is proportional to number of outgoing flavours.
 169   sigma  = (M_PI / sH2) * pow2(alpS) * nQuarkNew * sigSum;
 170
 171 }
 172
 173 //--------------------------------------------------------------------------
 174
 175 // Select identity, colour and anticolour.
 176
 177 void Sigma2gg2qqbar::setIdColAcol() {
 178
 179   // Flavours are trivial.
 180   setId( id1, id2, idNew, -idNew);
 181
 182   // Two colour flow topologies.
 183   double sigRand = sigSum * rndmPtr->flat();
 184   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
 185   else                 setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
 186
 187 }
 188
 189 //==========================================================================
 190
 191 // Sigma2qg2qg class.
 192 // Cross section for q g -> q g.
 193
 194 //--------------------------------------------------------------------------
 195
 196 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 197
 198 void Sigma2qg2qg::sigmaKin() {
 199
 200   // Calculate kinematics dependence.
 201   sigTS  = uH2 / tH2 - (4./9.) * uH / sH;
 202   sigTU  = sH2 / tH2 - (4./9.) * sH / uH;
 203   sigSum = sigTS + sigTU;
 204
 205   // Answer.
 206   sigma  = (M_PI / sH2) * pow2(alpS) * sigSum;
 207
 208 }
 209
 210 //--------------------------------------------------------------------------
 211
 212 // Select identity, colour and anticolour.
 213
 214 void Sigma2qg2qg::setIdColAcol() {
 215
 216   // Outgoing = incoming flavours.
 217   setId( id1, id2, id1, id2);
 218
 219   // Two colour flow topologies. Swap if first is gluon, or when antiquark.
 220   double sigRand = sigSum * rndmPtr->flat();
 221   if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 3, 0, 2, 3);
 222   else                 setColAcol( 1, 0, 2, 3, 2, 0, 1, 3);
 223   if (id1 == 21) swapCol1234();
 224   if (id1 < 0 || id2 < 0) swapColAcol();
 225
 226 }
 227
 228 //==========================================================================
 229
 230 // Sigma2qq2qq class.
 231 // Cross section for q qbar' -> q qbar' or q q' -> q q'
 232 // (qbar qbar' -> qbar qbar'), q' may be same as q.
 233
 234 //--------------------------------------------------------------------------
 235
 236 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 237
 238 void Sigma2qq2qq::sigmaKin() {
 239
 240   // Calculate kinematics dependence for different terms.
 241   sigT   = (4./9.) * (sH2 + uH2) / tH2;
 242   sigU   = (4./9.) * (sH2 + tH2) / uH2;
 243   sigTU  = - (8./27.) * sH2 / (tH * uH);
 244   sigST  = - (8./27.) * uH2 / (sH * tH);
 245
 246 }
 247
 248 //--------------------------------------------------------------------------
 249
 250
 251 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 252
 253 double Sigma2qq2qq::sigmaHat() {
 254
 255   // Combine cross section terms; factor 1/2 when identical quarks.
 256   if      (id2 ==  id1) sigSum = 0.5 * (sigT + sigU + sigTU);
 257   else if (id2 == -id1) sigSum = sigT + sigST;
 258   else                      sigSum = sigT;
 259
 260   // Answer.
 261   return (M_PI/sH2) * pow2(alpS) * sigSum;
 262
 263 }
 264
 265 //--------------------------------------------------------------------------
 266
 267 // Select identity, colour and anticolour.
 268
 269 void Sigma2qq2qq::setIdColAcol() {
 270
 271   // Outgoing = incoming flavours.
 272   setId( id1, id2, id1, id2);
 273
 274   // Colour flow topologies. Swap when antiquarks.
 275   if (id1 * id2 > 0)  setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
 276   else                setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
 277   if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT)
 278                       setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
 279   if (id1 < 0) swapColAcol();
 280
 281 }
 282
 283 //==========================================================================
 284
 285 // Sigma2qqbar2gg class.
 286 // Cross section for q qbar -> g g.
 287
 288 //--------------------------------------------------------------------------
 289
 290 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 291
 292 void Sigma2qqbar2gg::sigmaKin() {
 293
 294   // Calculate kinematics dependence.
 295   sigTS  = (32./27.) * uH / tH - (8./3.) * uH2 / sH2;
 296   sigUS  = (32./27.) * tH / uH - (8./3.) * tH2 / sH2;
 297   sigSum = sigTS + sigUS;
 298
 299   // Answer contains factor 1/2 from identical gluons.
 300   sigma  = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
 301
 302 }
 303
 304 //--------------------------------------------------------------------------
 305
 306 // Select identity, colour and anticolour.
 307
 308 void Sigma2qqbar2gg::setIdColAcol() {
 309
 310   // Outgoing flavours trivial.
 311   setId( id1, id2, 21, 21);
 312
 313   // Two colour flow topologies. Swap if first is antiquark.
 314   double sigRand = sigSum * rndmPtr->flat();
 315   if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
 316   else                 setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
 317   if (id1 < 0) swapColAcol();
 318
 319 }
 320
 321 //==========================================================================
 322
 323 // Sigma2qqbar2qqbarNew class.
 324 // Cross section q qbar -> q' qbar'.
 325
 326 //--------------------------------------------------------------------------
 327
 328 // Initialize process.
 329
 330 void Sigma2qqbar2qqbarNew::initProc() {
 331
 332   // Read number of quarks to be considered in massless approximation.
 333   nQuarkNew       = settingsPtr->mode("HardQCD:nQuarkNew");
 334
 335 }
 336
 337 //--------------------------------------------------------------------------
 338
 339 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 340
 341 void Sigma2qqbar2qqbarNew::sigmaKin() {
 342
 343   // Pick new flavour.
 344   idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
 345   mNew  = particleDataPtr->m0(idNew);
 346   m2New = mNew*mNew;
 347
 348   // Calculate kinematics dependence.
 349   sigS                      = 0.;
 350   if (sH > 4. * m2New) sigS = (4./9.) * (tH2 + uH2) / sH2;
 351
 352   // Answer is proportional to number of outgoing flavours.
 353   sigma = (M_PI / sH2) * pow2(alpS) * nQuarkNew * sigS;
 354
 355 }
 356
 357 //--------------------------------------------------------------------------
 358
 359 // Select identity, colour and anticolour.
 360
 361 void Sigma2qqbar2qqbarNew::setIdColAcol() {
 362
 363   // Set outgoing flavours ones.
 364   id3 = (id1 > 0) ? idNew : -idNew;
 365   setId( id1, id2, id3, -id3);
 366
 367   // Colour flow topologies. Swap when antiquarks.
 368   setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 369   if (id1 < 0) swapColAcol();
 370
 371 }
 372
 373 //==========================================================================
 374
 375 // Sigma2gg2QQbar class.
 376 // Cross section g g -> Q Qbar (Q = c, b or t).
 377 // Only provided for fixed m3 = m4 so do some gymnastics:
 378 // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
 379 // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
 380 //     but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
 381 //     tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
 382
 383 //--------------------------------------------------------------------------
 384
 385 // Initialize process.
 386
 387 void Sigma2gg2QQbar::initProc() {
 388
 389   // Process name.
 390   nameSave                 = "g g -> Q Qbar";
 391   if (idNew == 4) nameSave = "g g -> c cbar";
 392   if (idNew == 5) nameSave = "g g -> b bbar";
 393   if (idNew == 6) nameSave = "g g -> t tbar";
 394   if (idNew == 7) nameSave = "g g -> b' b'bar";
 395   if (idNew == 8) nameSave = "g g -> t' t'bar";
 396
 397   // Secondary open width fraction.
 398   openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
 399
 400 }
 401
 402 //--------------------------------------------------------------------------
 403
 404 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 405
 406 void Sigma2gg2QQbar::sigmaKin() {
 407
 408   // Modified Mandelstam variables for massive kinematics with m3 = m4.
 409   double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
 410   double tHQ    = -0.5 * (sH - tH + uH);
 411   double uHQ    = -0.5 * (sH + tH - uH);
 412   double tHQ2   = tHQ * tHQ;
 413   double uHQ2   = uHQ * uHQ;
 414
 415   // Calculate kinematics dependence.
 416   double tumHQ = tHQ * uHQ - s34Avg * sH;
 417   sigTS = ( uHQ / tHQ - 2.25 * uHQ2 / sH2 + 4.5 * s34Avg * tumHQ
 418     / ( sH * tHQ2) + 0.5 * s34Avg * (tHQ + s34Avg) / tHQ2
 419     - s34Avg*s34Avg / (sH * tHQ) ) / 6.;
 420   sigUS = ( tHQ / uHQ - 2.25 * tHQ2 / sH2 + 4.5 * s34Avg * tumHQ
 421     / ( sH * uHQ2) + 0.5 * s34Avg * (uHQ + s34Avg) / uHQ2
 422     - s34Avg*s34Avg / (sH * uHQ) ) / 6.;
 423   sigSum = sigTS + sigUS;
 424
 425   // Answer.
 426   sigma = (M_PI / sH2) * pow2(alpS) * sigSum * openFracPair;
 427
 428 }
 429
 430 //--------------------------------------------------------------------------
 431
 432 // Select identity, colour and anticolour.
 433
 434 void Sigma2gg2QQbar::setIdColAcol() {
 435
 436   // Flavours are trivial.
 437   setId( id1, id2, idNew, -idNew);
 438
 439   // Two colour flow topologies.
 440   double sigRand = sigSum * rndmPtr->flat();
 441   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
 442   else                 setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
 443
 444 }
 445
 446 //--------------------------------------------------------------------------
 447
 448 // Evaluate weight for decay angles of W in top decay.
 449
 450 double Sigma2gg2QQbar::weightDecay( Event& process, int iResBeg,
 451   int iResEnd) {
 452
 453   // For top decay hand over to standard routine, else done.
 454   if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
 455        return weightTopDecay( process, iResBeg, iResEnd);
 456   else return 1.;
 457
 458 }
 459
 460 //==========================================================================
 461
 462 // Sigma2qqbar2QQbar class.
 463 // Cross section q qbar -> Q Qbar (Q = c, b or t).
 464 // Only provided for fixed m3 = m4 so do some gymnastics:
 465 // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
 466 // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
 467 //     but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
 468 //     tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
 469
 470 //--------------------------------------------------------------------------
 471
 472 // Initialize process, especially parton-flux object.
 473
 474 void Sigma2qqbar2QQbar::initProc() {
 475
 476   // Process name.
 477   nameSave                 = "q qbar -> Q Qbar";
 478   if (idNew == 4) nameSave = "q qbar -> c cbar";
 479   if (idNew == 5) nameSave = "q qbar -> b bbar";
 480   if (idNew == 6) nameSave = "q qbar -> t tbar";
 481   if (idNew == 7) nameSave = "q qbar -> b' b'bar";
 482   if (idNew == 8) nameSave = "q qbar -> t' t'bar";
 483
 484   // Secondary open width fraction.
 485   openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
 486
 487 }
 488
 489 //--------------------------------------------------------------------------
 490
 491 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 492
 493 void Sigma2qqbar2QQbar::sigmaKin() {
 494
 495   // Modified Mandelstam variables for massive kinematics with m3 = m4.
 496   double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
 497   double tHQ    = -0.5 * (sH - tH + uH);
 498   double uHQ    = -0.5 * (sH + tH - uH);
 499   double tHQ2   = tHQ * tHQ;
 500   double uHQ2   = uHQ * uHQ;
 501
 502   // Calculate kinematics dependence.
 503   double sigS = (4./9.) * ((tHQ2 + uHQ2) / sH2 + 2. * s34Avg / sH);
 504
 505   // Answer.
 506   sigma = (M_PI / sH2) * pow2(alpS) * sigS * openFracPair;
 507
 508 }
 509
 510 //--------------------------------------------------------------------------
 511
 512 // Select identity, colour and anticolour.
 513
 514 void Sigma2qqbar2QQbar::setIdColAcol() {
 515
 516   // Set outgoing flavours.
 517   id3 = (id1 > 0) ? idNew : -idNew;
 518   setId( id1, id2, id3, -id3);
 519
 520   // Colour flow topologies. Swap when antiquarks.
 521   setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 522   if (id1 < 0) swapColAcol();
 523
 524 }
 525
 526 //--------------------------------------------------------------------------
 527
 528 // Evaluate weight for decay angles of W in top decay.
 529
 530 double Sigma2qqbar2QQbar::weightDecay( Event& process, int iResBeg,
 531   int iResEnd) {
 532
 533   // For top decay hand over to standard routine, else done.
 534   if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
 535        return weightTopDecay( process, iResBeg, iResEnd);
 536   else return 1.;
 537
 538 }
 539
 540
 541 //==========================================================================
 542
 543 // Sigma3gg2ggg class.
 544 // Cross section for g g -> g g g.
 545
 546 //--------------------------------------------------------------------------
 547
 548 // Evaluate |M|^2 - no incoming flavour dependence.
 549
 550 void Sigma3gg2ggg::sigmaKin() {
 551
 552   // Calculate all four-vector products.
 553   Vec4 p1cm( 0., 0.,  0.5 * mH, 0.5 * mH);
 554   Vec4 p2cm( 0., 0., -0.5 * mH, 0.5 * mH);
 555   pp[1][2] = p1cm * p2cm;
 556   pp[1][3] = p1cm * p3cm;
 557   pp[1][4] = p1cm * p4cm;
 558   pp[1][5] = p1cm * p5cm;
 559   pp[2][3] = p2cm * p3cm;
 560   pp[2][4] = p2cm * p4cm;
 561   pp[2][5] = p2cm * p5cm;
 562   pp[3][4] = p3cm * p4cm;
 563   pp[3][5] = p3cm * p5cm;
 564   pp[4][5] = p4cm * p5cm;
 565   for (int i = 1; i < 5; ++i)
 566     for (int j = i + 1; j < 6; ++j) pp[j][i] = pp[i][j];
 567
 568   // Cross section, in three main sections.
 569   double num1 = cycle(1,2,3,4,5) + cycle(1,2,3,5,4) + cycle(1,2,4,3,5)
 570               + cycle(1,2,4,5,3) + cycle(1,2,5,3,4) + cycle(1,2,5,4,3)
 571               + cycle(1,3,2,4,5) + cycle(1,3,2,5,4) + cycle(1,3,4,2,5)
 572               + cycle(1,3,5,2,4) + cycle(1,4,2,3,5) + cycle(1,4,3,2,5);
 573   double num2 = pow4(pp[1][2]) + pow4(pp[1][3]) + pow4(pp[1][4])
 574               + pow4(pp[1][5]) + pow4(pp[2][3]) + pow4(pp[2][4])
 575               + pow4(pp[2][5]) + pow4(pp[3][4]) + pow4(pp[3][5])
 576               + pow4(pp[4][5]);
 577   double den  = pp[1][2] * pp[1][3] * pp[1][4] * pp[1][5] * pp[2][3]
 578               * pp[2][4] * pp[2][5] * pp[3][4] * pp[3][5] * pp[4][5];
 579
 580   // Answer has a factor 6 due to identical gluons
 581   // This is cancelled by phase space factor (1 / 6)
 582   sigma = pow3(4. * M_PI * alpS) * (27./16.) * num1 * num2 / den;
 583
 584 }
 585
 586 //--------------------------------------------------------------------------
 587
 588 // Select identity, colour and anticolour.
 589
 590 void Sigma3gg2ggg::setIdColAcol() {
 591
 592   // Flavours are trivial.
 593   setId( id1, id2, 21, 21, 21);
 594
 595   // Three colour flow topologies, each with two orientations.
 596   /*
 597   double sigRand = sigSum * rndmPtr->flat();
 598   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
 599   else if (sigRand < sigTS + sigUS)
 600                        setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
 601   else                 setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
 602   if (rndmPtr->flat() > 0.5) swapColAcol();
 603   */
 604
 605   // Temporary solution.
 606   setColAcol( 1, 2, 2, 3, 1, 4, 4, 5, 5, 3);
 607 }
 608
 609
 610 //==========================================================================
 611
 612 // Sigma3qqbar2ggg class.
 613 // Cross section for q qbar -> g g g.
 614
 615 //--------------------------------------------------------------------------
 616
 617 // Evaluate |M|^2 - no incoming flavour dependence.
 618 void Sigma3qqbar2ggg::sigmaKin() {
 619
 620   // Setup four-vectors
 621   pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
 622   pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
 623   pCM[2] = p3cm;
 624   pCM[3] = p4cm;
 625   pCM[4] = p5cm;
 626
 627   // Calculate |M|^2
 628   // Answer has a factor 6 due to identical gluons,
 629   // which is cancelled by phase space factor (1 / 6)
 630   sigma = m2Calc();
 631
 632 }
 633
 634 //--------------------------------------------------------------------------
 635
 636 // |M|^2
 637
 638 inline double Sigma3qqbar2ggg::m2Calc() {
 639
 640   // Calculate four-products
 641   double sHnow  = (pCM[0] + pCM[1]).m2Calc();
 642   double sHhalf = sH / 2.;
 643
 644   // qbar (p+) + q(p-) -> g(k1) g(k2) g(k3)
 645   // a_i = (p+ . ki), i = 1, 2, 3
 646   // b_i = (p- . ki), i = 1, 2, 3
 647   a[0] = pCM[0] * pCM[2];
 648   a[1] = pCM[0] * pCM[3];
 649   a[2] = pCM[0] * pCM[4];
 650   b[0] = pCM[1] * pCM[2];
 651   b[1] = pCM[1] * pCM[3];
 652   b[2] = pCM[1] * pCM[4];
 653
 654   pp[0][1] = pCM[2] * pCM[3];
 655   pp[1][2] = pCM[3] * pCM[4];
 656   pp[2][0] = pCM[4] * pCM[2];
 657
 658   // ab[i][j] = a_i * b_j + a_j * b_i
 659   ab[0][1] = a[0] * b[1] + a[1] * b[0];
 660   ab[1][2] = a[1] * b[2] + a[2] * b[1];
 661   ab[2][0] = a[2] * b[0] + a[0] * b[2];
 662
 663   // Cross section
 664   double num1 = a[0] * b[0] * (a[0] * a[0] + b[0] * b[0]) +
 665                 a[1] * b[1] * (a[1] * a[1] + b[1] * b[1]) +
 666                 a[2] * b[2] * (a[2] * a[2] + b[2] * b[2]);
 667   double den1 = a[0] * a[1] * a[2] * b[0] * b[1] * b[2];
 668   double num2 = - ( ab[0][1] / pp[0][1] )
 669                 - ( ab[1][2] / pp[1][2] )
 670                 - ( ab[2][0] / pp[2][0] );
 671   double num3 = a[2] * b[2] * ab[0][1] / (pp[1][2] * pp[2][0] ) +
 672                 a[0] * b[0] * ab[1][2] / (pp[2][0] * pp[0][1] ) +
 673                 a[1] * b[1] * ab[2][0] / (pp[0][1] * pp[1][2] );
 674
 675   // Final answer
 676   return  pow3(4. * M_PI * alpS) * (8. / 324.) * (num1 / den1) *
 677           ( sHhalf + 9. * (sHhalf + num2) + (2. * 81. / sHnow) * num3 );
 678
 679 }
 680
 681 //--------------------------------------------------------------------------
 682
 683 // Select identity, colour and anticolour.
 684
 685 void Sigma3qqbar2ggg::setIdColAcol(){
 686
 687   // Flavours are trivial.
 688   setId( id1, id2, 21, 21, 21);
 689
 690   // Temporary solution.
 691   setColAcol( 1, 0, 0, 2, 1, 3, 3, 4, 4, 2);
 692   if (id1 < 0) swapColAcol();
 693 }
 694
 695 //--------------------------------------------------------------------------
 696
 697 // Map a final state configuration
 698
 699 inline void Sigma3qqbar2ggg::mapFinal() {
 700   switch (config) {
 701   case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
 702   case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
 703   case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
 704   case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
 705   case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
 706   case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
 707   }
 708 }
 709
 710 //==========================================================================
 711
 712 // Sigma3qg2qgg class.
 713 // Cross section for q g -> q g g.
 714 // Crossed relation from q qbar -> g g g:
 715 //   qbar(p+) q(p-) -> g(k1) g(k2) g(k3)
 716
 717 //--------------------------------------------------------------------------
 718
 719 // Evaluate |M|^2 - no incoming flavour dependence
 720 // Note: two different contributions from gq and qg incoming
 721
 722 void Sigma3qg2qgg::sigmaKin() {
 723
 724   // Pick a final state configuration
 725   pickFinal();
 726
 727   // gq and qg incoming
 728   for (int i = 0; i < 2; i++) {
 729
 730     // Map incoming four-vectors to p+, p-, k1, k2, k3
 731     pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
 732     pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
 733     mapFinal();
 734
 735     // Crossing
 736     swap(pCM[i], pCM[2]);
 737
 738     // |M|^2
 739     // XXX - Extra factor of (3) from selecting a final state
 740     // configuration (already a factor of 2 in the original
 741     // answer due to two identical final state gluons)???
 742     // Extra factor of (3 / 8) as average over incoming gluon
 743     sigma[i] = 3. * (3. / 8.) * m2Calc();
 744
 745   } // for (int i = 0; i < 2; i++)
 746
 747 }
 748
 749 //--------------------------------------------------------------------------
 750
 751 // Evaluate |M|^2 - incoming flavour dependence
 752 // Pick from two configurations calculated previously
 753
 754 double Sigma3qg2qgg::sigmaHat() {
 755   // gq or qg incoming
 756   return (id1 == 21) ? sigma[0] : sigma[1];
 757 }
 758
 759 //--------------------------------------------------------------------------
 760
 761 // Select identity, colour and anticolour.
 762
 763 void Sigma3qg2qgg::setIdColAcol(){
 764   // Outgoing flavours; only need to know where the quark is
 765   int qIdx    = config / 2;
 766   int idTmp[3] = { 21, 21, 21 };
 767   idTmp[qIdx]  = (id1 == 21) ? id2 : id1;
 768   setId( id1, id2, idTmp[0], idTmp[1], idTmp[2]);
 769
 770   // Temporary solution
 771   if      (qIdx == 0) setColAcol(1, 0, 2, 1, 4, 0, 3, 4, 2, 3);
 772   else if (qIdx == 1) setColAcol(1, 0, 2, 1, 3, 4, 4, 0, 2, 3);
 773   else                setColAcol(1, 0, 2, 1, 3, 4, 2, 3, 4, 0);
 774   // gq or qg incoming
 775   if (id1 == 21) {
 776     swap( colSave[1], colSave[2]);
 777     swap(acolSave[1], acolSave[2]);
 778   }
 779   // qbar rather than q incoming
 780   if (id1 < 0 || id2 < 0) swapColAcol();
 781
 782 }
 783
 784 //==========================================================================
 785
 786 // Sigma3gg2qqbarg class.
 787 // Cross section for g g -> q qbar g
 788 // Crossed relation from q qbar -> g g g
 789
 790 //--------------------------------------------------------------------------
 791
 792 // Initialize process.
 793
 794 void Sigma3gg2qqbarg::initProc() {
 795
 796   // Read number of quarks to be considered in massless approximation.
 797   nQuarkNew       = double(settingsPtr->mode("HardQCD:nQuarkNew"));
 798
 799 }
 800
 801 //--------------------------------------------------------------------------
 802
 803 // Evaluate |M|^2 - no incoming flavour dependence.
 804
 805 void Sigma3gg2qqbarg::sigmaKin() {
 806
 807   // Incoming four-vectors
 808   pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
 809   pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
 810
 811   // Pick and map a final state configuration
 812   pickFinal();
 813   mapFinal();
 814
 815   // Crossing
 816   swap(pCM[0], pCM[2]);
 817   swap(pCM[1], pCM[3]);
 818
 819   // |M|^2
 820   // Extra factor of (6.) from picking a final state configuration
 821   // Extra factor of nQuarkNew
 822   // Extra factor of (3. / 8.) ^ 2 as averaging over two incoming gluons
 823   sigma = 6. * double(nQuarkNew) * (3. / 8.) * (3. / 8.) * m2Calc();
 824
 825 }
 826
 827 //--------------------------------------------------------------------------
 828
 829 // Select identity, colour and anticolour.
 830
 831 void Sigma3gg2qqbarg::setIdColAcol(){
 832
 833   // Pick new flavour
 834   int idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
 835
 836   // Outgoing flavours; easiest just to map by hand
 837   switch (config) {
 838   case 0: id3 =  idNew; id4 = -idNew; id5 =  21;    break;
 839   case 1: id3 =  idNew; id4 =  21;    id5 = -idNew; break;
 840   case 2: id3 = -idNew; id4 =  idNew; id5 =  21;    break;
 841   case 3: id3 =  21;    id4 =  idNew; id5 = -idNew; break;
 842   case 4: id3 = -idNew; id4 =  21;    id5 =  idNew; break;
 843   case 5: id3 =  21;    id4 = -idNew; id5 =  idNew; break;
 844   }
 845   setId(id1, id2, id3, id4, id5);
 846
 847   // Temporary solution
 848   switch (config) {
 849   case 0: setColAcol( 1, 2, 2, 3, 4, 0, 0, 3, 1, 4 ); break;
 850   case 1: setColAcol( 1, 2, 2, 3, 4, 0, 1, 4, 0, 3 ); break;
 851   case 2: setColAcol( 1, 2, 2, 3, 0, 3, 4, 0, 1, 4 ); break;
 852   case 3: setColAcol( 1, 2, 2, 3, 1, 4, 4, 0, 0, 3 ); break;
 853   case 4: setColAcol( 1, 2, 2, 3, 0, 3, 1, 4, 4, 0 ); break;
 854   case 5: setColAcol( 1, 2, 2, 3, 1, 4, 0, 3, 4, 0 ); break;
 855   }
 856
 857 }
 858
 859 //==========================================================================
 860
 861 // Sigma3qq2qqgDiff class.
 862 // Cross section for q q' -> q q' g, q != q'
 863
 864 //--------------------------------------------------------------------------
 865
 866 // Evaluate |M|^2 - no incoming flavour dependence
 867
 868 void Sigma3qq2qqgDiff::sigmaKin() {
 869
 870   // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
 871
 872   // Incoming four-vectors
 873   pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
 874   pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
 875   // Pick and map a final state configuration
 876   pickFinal();
 877   mapFinal();
 878
 879   // |M|^2
 880   // Extra factor of (6.) from picking a final state configuration
 881   sigma = 6. * m2Calc();
 882 }
 883
 884 //--------------------------------------------------------------------------
 885
 886 // |M|^2
 887
 888 inline double Sigma3qq2qqgDiff::m2Calc() {
 889
 890   // Four-products
 891   s  = (pCM[0] + pCM[1]).m2Calc();
 892   t  = (pCM[0] - pCM[2]).m2Calc();
 893   u  = (pCM[0] - pCM[3]).m2Calc();
 894   up = (pCM[1] - pCM[2]).m2Calc();
 895   sp = (pCM[2] + pCM[3]).m2Calc();
 896   tp = (pCM[1] - pCM[3]).m2Calc();
 897
 898   // |M|^2
 899   double num1 = (s * s + sp * sp + u * u + up * up) / (t * tp);
 900   double den1 = (pCM[0] * pCM[4]) * (pCM[1] * pCM[4]) *
 901                 (pCM[2] * pCM[4]) * (pCM[3] * pCM[4]);
 902   double num2 = (u + up) * (s * sp + t * tp - u * up) +
 903                 u * (s * t + sp * tp) + up * (s * tp + sp * t);
 904   double num3 = (s + sp) * (s * sp - t * tp - u * up) +
 905                 2. * t * tp * (u + up) + 2. * u * up * (t + tp);
 906
 907   // (N^2 - 1)^2 / 4N^3 = 16. / 27.
 908   // (N^2 - 1)   / 4N^3 =  2. / 27.
 909   return (1. / 8.) * pow3(4. * M_PI * alpS) * num1 / den1 *
 910          ( (16. / 27.) * num2 - (2. / 27.) * num3 );
 911
 912 }
 913
 914 //--------------------------------------------------------------------------
 915
 916 // Evaluate |M|^2 - incoming flavour dependence
 917
 918 double Sigma3qq2qqgDiff::sigmaHat() {
 919   // Different incoming flavours only
 920   if (abs(id1) == abs(id2)) return 0.;
 921   return sigma;
 922 }
 923
 924 //--------------------------------------------------------------------------
 925
 926 // Select identity, colour and anticolour.
 927
 928 void Sigma3qq2qqgDiff::setIdColAcol(){
 929
 930   // Outgoing flavours; easiest just to map by hand
 931   switch (config) {
 932   case 0: id3 = id1; id4 = id2; id5 = 21;  break;
 933   case 1: id3 = id1; id4 = 21;  id5 = id2; break;
 934   case 2: id3 = id2; id4 = id1; id5 = 21;  break;
 935   case 3: id3 = 21;  id4 = id1; id5 = id2; break;
 936   case 4: id3 = id2; id4 = 21;  id5 = id1; break;
 937   case 5: id3 = 21;  id4 = id2; id5 = id1; break;
 938   }
 939   setId(id1, id2, id3, id4, id5);
 940
 941   // Temporary solution; id1 and id2 can be q/qbar independently
 942   int cols[5][2];
 943   if (id1 > 0) {
 944     cols[0][0] = 1; cols[0][1] = 0;
 945     cols[2][0] = 1; cols[2][1] = 0;
 946   } else {
 947     cols[0][0] = 0; cols[0][1] = 1;
 948     cols[2][0] = 0; cols[2][1] = 1;
 949   }
 950   if (id2 > 0) {
 951     cols[1][0] = 2; cols[1][1] = 0;
 952     cols[3][0] = 3; cols[3][1] = 0;
 953     cols[4][0] = 2; cols[4][1] = 3;
 954   } else {
 955     cols[1][0] = 0; cols[1][1] = 2;
 956     cols[3][0] = 0; cols[3][1] = 3;
 957     cols[4][0] = 3; cols[4][1] = 2;
 958   }
 959   // Map correct final state configuration
 960   int i3 = 0, i4 = 0, i5 = 0;
 961   switch (config) {
 962   case 0: i3 = 2; i4 = 3; i5 = 4; break;
 963   case 1: i3 = 2; i4 = 4; i5 = 3; break;
 964   case 2: i3 = 3; i4 = 2; i5 = 4; break;
 965   case 3: i3 = 4; i4 = 2; i5 = 3; break;
 966   case 4: i3 = 3; i4 = 4; i5 = 2; break;
 967   case 5: i3 = 4; i4 = 3; i5 = 2; break;
 968   }
 969   // Put colours in place
 970   setColAcol(cols[0][0],  cols[0][1],  cols[1][0],  cols[1][1],
 971              cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
 972              cols[i5][0], cols[i5][1]);
 973
 974 }
 975
 976 //--------------------------------------------------------------------------
 977
 978 // Map a final state configuration
 979
 980 inline void Sigma3qq2qqgDiff::mapFinal() {
 981   switch (config) {
 982   case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
 983   case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
 984   case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
 985   case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
 986   case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
 987   case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
 988   }
 989 }
 990
 991
 992 //==========================================================================
 993
 994 // Sigma3qqbar2qqbargDiff
 995 // Cross section for q qbar -> q' qbar' g
 996 // Crossed relation from q q' -> q q' g, q != q'
 997
 998 //--------------------------------------------------------------------------
 999
1000 // Initialize process.
1001
1002 void Sigma3qqbar2qqbargDiff::initProc() {
1003
1004   // Read number of quarks to be considered in massless approximation.
1005   nQuarkNew       = double(settingsPtr->mode("HardQCD:nQuarkNew"));
1006
1007 }
1008
1009 //--------------------------------------------------------------------------
1010
1011 // Evaluate |M|^2 - no incoming flavour dependence.
1012
1013 void Sigma3qqbar2qqbargDiff::sigmaKin() {
1014   // Overall 6 possibilities for final state ordering
1015   // To keep symmetry between final states, always map to:
1016   //  1) q1(p+)    qbar1(p-)  -> qbar2(q+)  q2(q-)     g(k)
1017   //  2) qbar1(p+) q1(p-)     -> q2(q+)     qbar2(q-)  g(k)
1018   // Crossing p- and q+ gives:
1019   //  1) q1(p+)    q2(-q+)    -> q1(-p-)    q2(q-)     g(k)
1020   //  2) qbar1(p+) qbar2(-q+) -> qbar1(-p-) qbar2(q-)  g(k)
1021
1022   // Incoming four-vectors
1023   pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
1024   pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1025   // Pick and map a final state configuration
1026   pickFinal();
1027   mapFinal();
1028
1029   // Crossing
1030   swap(pCM[1], pCM[2]);
1031   pCM[1] = -pCM[1];
1032   pCM[2] = -pCM[2];
1033
1034   // |M|^2
1035   // Extra factor of (6.) from picking a final state configuration
1036   // Extra factor of (nQuarkNew - 1) from new q/qbar pairs
1037   // XXX - Extra factor of (2.) from second possible crossing???
1038   sigma = 6. * (nQuarkNew - 1) * 2. * m2Calc();
1039
1040 }
1041
1042 //--------------------------------------------------------------------------
1043
1044 // Select identity, colour and anticolour.
1045
1046 void Sigma3qqbar2qqbargDiff::setIdColAcol(){
1047
1048   // Pick new q qbar flavour with incoming flavour disallowed
1049   int idNew = 1 + int( (nQuarkNew - 1) * rndmPtr->flat() );
1050   if (idNew >= abs(id1)) ++idNew;
1051   // For qbar q incoming, q+ is always mapped to q2
1052   // For q qbar incoming, q+ is always mapped to qbar2
1053   if (id1 > 0) idNew = -idNew;
1054
1055   // Outgoing flavours; easiest just to map by hand
1056   switch (config) {
1057   case 0: id3 =  idNew; id4 = -idNew; id5 =  21;    break;
1058   case 1: id3 =  idNew; id4 =  21;    id5 = -idNew; break;
1059   case 2: id3 = -idNew; id4 =  idNew; id5 =  21;    break;
1060   case 3: id3 =  21;    id4 =  idNew; id5 = -idNew; break;
1061   case 4: id3 = -idNew; id4 =  21;    id5 =  idNew; break;
1062   case 5: id3 =  21;    id4 = -idNew; id5 =  idNew; break;
1063   }
1064   setId(id1, id2, id3, id4, id5);
1065
1066   // Temporary solution; start with q qbar -> qbar q g
1067   int cols[5][2];
1068   cols[0][0] = 1; cols[0][1] = 0;
1069   cols[1][0] = 0; cols[1][1] = 2;
1070   cols[2][0] = 0; cols[2][1] = 3;
1071   cols[3][0] = 1; cols[3][1] = 0;
1072   cols[4][0] = 3; cols[4][1] = 2;
1073   // Map into correct place
1074   int i3 = 0, i4 = 0, i5 = 0;
1075   switch (config) {
1076   case 0: i3 = 2; i4 = 3; i5 = 4; break;
1077   case 1: i3 = 2; i4 = 4; i5 = 3; break;
1078   case 2: i3 = 3; i4 = 2; i5 = 4; break;
1079   case 3: i3 = 4; i4 = 2; i5 = 3; break;
1080   case 4: i3 = 3; i4 = 4; i5 = 2; break;
1081   case 5: i3 = 4; i4 = 3; i5 = 2; break;
1082   }
1083   setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1084              cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1085              cols[i5][0], cols[i5][1]);
1086   // Swap for qbar q incoming
1087   if (id1 < 0) swapColAcol();
1088
1089 }
1090
1091 //==========================================================================
1092
1093 // Sigma3qg2qqqbarDiff class.
1094 // Cross section for q g -> q q' qbar'
1095 // Crossed relation from q q' -> q q' g, q != q'
1096 //   q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1097
1098 //--------------------------------------------------------------------------
1099
1100 // Initialize process.
1101
1102 void Sigma3qg2qqqbarDiff::initProc() {
1103
1104   // Read number of quarks to be considered in massless approximation.
1105   nQuarkNew       = settingsPtr->mode("HardQCD:nQuarkNew");
1106
1107 }
1108
1109 //--------------------------------------------------------------------------
1110
1111 // Evaluate |M|^2 - no incoming flavour dependence
1112
1113 void Sigma3qg2qqqbarDiff::sigmaKin() {
1114
1115   // Pick a final state configuration
1116   pickFinal();
1117
1118   // gq or qg incoming
1119   for (int i = 0; i < 2; i++) {
1120
1121     // Map incoming four-vectors
1122     pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
1123     pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1124     mapFinal();
1125
1126     // Crossing (note extra -ve sign in total sigma)
1127     swap(pCM[i], pCM[4]);
1128     pCM[i] = -pCM[i];
1129     pCM[4] = -pCM[4];
1130
1131     // |M|^2
1132     // Extra factor of (6) from picking a final state configuration
1133     // Extra factor of (3 / 8) as averaging over incoming gluon
1134     // Extra factor of (nQuarkNew - 1) due to new q/qbar pair
1135     sigma[i] = -6. * (3. / 8.) * (nQuarkNew - 1) * m2Calc();
1136
1137   }
1138
1139 }
1140
1141 //--------------------------------------------------------------------------
1142
1143 // Evaluate |M|^2 - incoming flavour dependence
1144
1145 double Sigma3qg2qqqbarDiff::sigmaHat() {
1146   // gq or qg incoming
1147   return (id1 == 21) ? sigma[0] : sigma[1];
1148 }
1149
1150 //--------------------------------------------------------------------------
1151
1152 // Select identity, colour and anticolour.
1153
1154 void Sigma3qg2qqqbarDiff::setIdColAcol(){
1155   // Pick new q qbar flavour with incoming flavour disallowed
1156   int sigmaIdx = (id1 == 21) ? 0 : 1;
1157   int idIn     = (id1 == 21) ? id2 : id1;
1158   int idNew    = 1 + int( (nQuarkNew - 1) * rndmPtr->flat() );
1159   if (idNew >= abs(idIn)) ++idNew;
1160
1161   // qbar instead of q incoming means swap outgoing q/qbar pair
1162   int id3Tmp = idIn, id4Tmp = idNew, id5Tmp = -idNew;
1163   if (idIn < 0)  swap(id4Tmp, id5Tmp);
1164   // If g q incoming rather than q g, idIn and idNew
1165   // should be exchanged (see sigmaKin)
1166   if (sigmaIdx == 0) swap(id3Tmp, id4Tmp);
1167   // Outgoing flavours; now just map as if q g incoming
1168   switch (config) {
1169   case 0: id3 = id3Tmp; id4 = id4Tmp; id5 = id5Tmp; break;
1170   case 1: id3 = id3Tmp; id4 = id5Tmp; id5 = id4Tmp; break;
1171   case 2: id3 = id4Tmp; id4 = id3Tmp; id5 = id5Tmp; break;
1172   case 3: id3 = id5Tmp; id4 = id3Tmp; id5 = id4Tmp; break;
1173   case 4: id3 = id4Tmp; id4 = id5Tmp; id5 = id3Tmp; break;
1174   case 5: id3 = id5Tmp; id4 = id4Tmp; id5 = id3Tmp; break;
1175   }
1176   setId(id1, id2, id3, id4, id5);
1177
1178   // Temporary solution; start with either
1179   // g q1    -> q1    q2 qbar2
1180   // g qbar1 -> qbar1 qbar2 q2
1181   int cols[5][2];
1182   cols[0][0] = 1; cols[0][1] = 2;
1183   if (idIn > 0) {
1184     cols[1][0] = 3; cols[1][1] = 0;
1185     cols[2][0] = 1; cols[2][1] = 0;
1186     cols[3][0] = 3; cols[3][1] = 0;
1187     cols[4][0] = 0; cols[4][1] = 2;
1188   } else {
1189     cols[1][0] = 0; cols[1][1] = 3;
1190     cols[2][0] = 0; cols[2][1] = 2;
1191     cols[3][0] = 0; cols[3][1] = 3;
1192     cols[4][0] = 1; cols[4][1] = 0;
1193   }
1194   // Swap incoming if q/qbar g instead
1195   if (id2 == 21) {
1196     swap(cols[0][0], cols[1][0]);
1197     swap(cols[0][1], cols[1][1]);
1198   }
1199   // Map final state
1200   int i3 = 0, i4 = 0, i5 = 0;
1201   if (sigmaIdx == 0) {
1202     switch (config) {
1203     case 0: i3 = 3; i4 = 2; i5 = 4; break;
1204     case 1: i3 = 3; i4 = 4; i5 = 2; break;
1205     case 2: i3 = 2; i4 = 3; i5 = 4; break;
1206     case 3: i3 = 4; i4 = 3; i5 = 2; break;
1207     case 4: i3 = 2; i4 = 4; i5 = 3; break;
1208     case 5: i3 = 4; i4 = 2; i5 = 3; break;
1209     }
1210   } else {
1211     switch (config) {
1212     case 0: i3 = 2; i4 = 3; i5 = 4; break;
1213     case 1: i3 = 2; i4 = 4; i5 = 3; break;
1214     case 2: i3 = 3; i4 = 2; i5 = 4; break;
1215     case 3: i3 = 4; i4 = 2; i5 = 3; break;
1216     case 4: i3 = 3; i4 = 4; i5 = 2; break;
1217     case 5: i3 = 4; i4 = 3; i5 = 2; break;
1218     }
1219   }
1220   setColAcol(cols[0][0],  cols[0][1],  cols[1][0],  cols[1][1],
1221              cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1222              cols[i5][0], cols[i5][1]);
1223 }
1224
1225 //==========================================================================
1226
1227 // Sigma3qq2qqgSame class.
1228 // Cross section for q q' -> q q' g, q == q'.
1229
1230 //--------------------------------------------------------------------------
1231
1232 // Evaluate |M|^2 - no incoming flavour dependence
1233
1234 void Sigma3qq2qqgSame::sigmaKin() {
1235   // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1236
1237   // Incoming four-vectors
1238   pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
1239   pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1240   // Pick/map a final state configuration
1241   pickFinal();
1242   mapFinal();
1243
1244   // |M|^2
1245   // Extra factor (3) from picking final state configuration
1246   // (original answer already has a factor 2 from identical
1247   // quarks in the final state)
1248   sigma = 3. * m2Calc();
1249
1250 }
1251
1252 //--------------------------------------------------------------------------
1253
1254 // |M|^2
1255
1256 inline double Sigma3qq2qqgSame::m2Calc() {
1257
1258   // Four-products
1259   s  = (pCM[0] + pCM[1]).m2Calc();
1260   t  = (pCM[0] - pCM[2]).m2Calc();
1261   u  = (pCM[0] - pCM[3]).m2Calc();
1262   sp = (pCM[2] + pCM[3]).m2Calc();
1263   tp = (pCM[1] - pCM[3]).m2Calc();
1264   up = (pCM[1] - pCM[2]).m2Calc();
1265
1266   // |M|^2
1267   ssp  = s * sp;
1268   ttp  = t * tp;
1269   uup  = u * up;
1270   s_sp = s + sp;
1271   t_tp = t + tp;
1272   u_up = u + up;
1273
1274   double den1 = (pCM[0] * pCM[4]) * (pCM[1] * pCM[4]) *
1275                 (pCM[2] * pCM[4]) * (pCM[3] * pCM[4]);
1276
1277   double fac1 = s * (t * u + tp * up) + sp * (t * up + tp * u);
1278   double fac2 = ssp - ttp - uup;
1279   double fac3 = 2. * (ttp * u_up + uup * t_tp);
1280
1281   double num1 = u_up * (ssp + ttp - uup) + fac1;
1282   double num2 = s_sp * fac2 + fac3;
1283   double num3 = (s * s + sp * sp + u * u + up * up) / (t * tp);
1284
1285   double num4 = t_tp * (ssp - ttp + uup) + fac1;
1286   double num5 = (s * s + sp * sp + t * t + tp * tp) / (u * up);
1287
1288   double num6 = s_sp * fac2 - fac3 - 2. * fac1;
1289   double num7 = (s * s + sp * sp) * fac2;
1290   double den7 = (ttp * uup);
1291
1292   // C1 = (N^2 - 1)^2 / 4N^3 = 16. / 27.
1293   // C2 = (N^2 - 1)   / 4N^3 =  2. / 27.
1294   // C3 = (N^4 - 1)   / 8N^4 = 10. / 81.
1295   // C4 = (N^2 - 1)^2 / 8N^4 =  8. / 81.
1296   return (1. / 8.) * pow3(4. * M_PI * alpS) *
1297          ( ( (16. / 27.) * num1 - (2. / 27.) * num2 ) * num3 +
1298          ( (16. / 27.) * num4 - (2. / 27.) * num2 ) * num5 +
1299          ( (10. / 81.) * num2 + (8. / 81.) * num6 ) *
1300          ( num7 / den7 ) ) / den1;
1301
1302 }
1303
1304 //--------------------------------------------------------------------------
1305
1306 // Evaluate |M|^2 - incoming flavour dependence
1307
1308 double Sigma3qq2qqgSame::sigmaHat() {
1309   // q q / qbar qbar incoming states only
1310   if (id1 != id2) return 0.;
1311   return sigma;
1312 }
1313
1314 //--------------------------------------------------------------------------
1315
1316 // Select identity, colour and anticolour.
1317
1318 void Sigma3qq2qqgSame::setIdColAcol(){
1319
1320   // Need to know where the gluon was mapped (pCM[4])
1321   int gIdx = 0;
1322   switch (config) {
1323   case 3: case 5: gIdx = 0; break;
1324   case 1: case 4: gIdx = 1; break;
1325   case 0: case 2: gIdx = 2; break;
1326   }
1327
1328   // Outgoing flavours
1329   int idTmp[3] = { id1, id1, id1 };
1330   idTmp[gIdx]  = 21;
1331   setId(id1, id2, idTmp[0], idTmp[1], idTmp[2]);
1332
1333   // Temporary solution; start with q q -> q q g
1334   setColAcol(1, 0, 2, 0, 1, 0, 3, 0, 2, 3);
1335   // Map gluon
1336   swap( colSave[5],  colSave[gIdx + 3]);
1337   swap(acolSave[5], acolSave[gIdx + 3]);
1338   // Swap if qbar qbar incoming
1339   if (id1 < 0) swapColAcol();
1340
1341 }
1342
1343 //--------------------------------------------------------------------------
1344
1345 // Map a final state configuration
1346 inline void Sigma3qq2qqgSame::mapFinal() {
1347   switch (config) {
1348   case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
1349   case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
1350   case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
1351   case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
1352   case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
1353   case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
1354   }
1355 }
1356
1357 //==========================================================================
1358
1359 // Sigma3qqbar2qqbargSame class.
1360 // Cross section for q qbar -> q qbar g
1361 // Crossed relation from q(bar) q(bar) -> q(bar) q(bar) g:
1362 //   q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1363 //--------------------------------------------------------------------------
1364
1365 // Evaluate |M|^2 - no incoming flavour dependence
1366
1367 void Sigma3qqbar2qqbargSame::sigmaKin() {
1368
1369   // Incoming four-vectors
1370   pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
1371   pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1372
1373   // Pick and map a final state configuration
1374   pickFinal();
1375   mapFinal();
1376
1377   // Crossing
1378   swap(pCM[1], pCM[3]);
1379   pCM[1] = -pCM[1];
1380   pCM[3] = -pCM[3];
1381
1382   // |M|^2
1383   // Extra factor of (6) from picking a final state configuration
1384   sigma = 6. * m2Calc();
1385
1386 }
1387
1388 //--------------------------------------------------------------------------
1389
1390 // Select identity, colour and anticolour.
1391
1392 void Sigma3qqbar2qqbargSame::setIdColAcol(){
1393   // Outgoing flavours; easiest to map by hand
1394   switch (config) {
1395   case 0: id3 = id1; id4 = id2; id5 = 21;  break;
1396   case 1: id3 = id1; id4 = 21;  id5 = id2; break;
1397   case 2: id3 = id2; id4 = id1; id5 = 21;  break;
1398   case 3: id3 = 21;  id4 = id1; id5 = id2; break;
1399   case 4: id3 = id2; id4 = 21;  id5 = id1; break;
1400   case 5: id3 = 21;  id4 = id2; id5 = id1; break;
1401   }
1402   setId(id1, id2, id3, id4, id5);
1403
1404   // Temporary solution; start with q qbar -> q qbar g
1405   int cols[5][2];
1406   cols[0][0] = 1; cols[0][1] = 0;
1407   cols[1][0] = 0; cols[1][1] = 2;
1408   cols[2][0] = 1; cols[2][1] = 0;
1409   cols[3][0] = 0; cols[3][1] = 3;
1410   cols[4][0] = 3; cols[4][1] = 2;
1411   // Map final state
1412   int i3 = 0, i4 = 0, i5 = 0;
1413   switch (config) {
1414   case 0: i3 = 2; i4 = 3; i5 = 4; break;
1415   case 1: i3 = 2; i4 = 4; i5 = 3; break;
1416   case 2: i3 = 3; i4 = 2; i5 = 4; break;
1417   case 3: i3 = 4; i4 = 2; i5 = 3; break;
1418   case 4: i3 = 3; i4 = 4; i5 = 2; break;
1419   case 5: i3 = 4; i4 = 3; i5 = 2; break;
1420   }
1421   setColAcol(cols[0][0],  cols[0][1],  cols[1][0],  cols[1][1],
1422              cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1423              cols[i5][0], cols[i5][1]);
1424   // Swap for qbar q incoming
1425   if (id1 < 0) swapColAcol();
1426 }
1427
1428 //==========================================================================
1429
1430 // Sigma3qg2qqqbarSame class.
1431 // Cross section for q g -> q q qbar.
1432 // Crossed relation from q(bar) q(bar) -> q(bar) q(bar) g:
1433 //   q1(p+) q1(p-) -> q1(q+) q1(q-) g(k)
1434
1435 //--------------------------------------------------------------------------
1436
1437 // Evaluate |M|^2 - no incoming flavour dependence
1438
1439 void Sigma3qg2qqqbarSame::sigmaKin() {
1440
1441   // Pick a final state configuration
1442   pickFinal();
1443
1444   // gq and qg incoming
1445   for (int i = 0; i < 2; i++) {
1446
1447     // Map incoming four-vectors
1448     pCM[0] = Vec4( 0., 0.,  0.5 * mH, 0.5 * mH);
1449     pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1450     mapFinal();
1451
1452     // Crossing (note extra -ve sign in total sigma)
1453     swap(pCM[i], pCM[4]);
1454     pCM[i] = -pCM[i];
1455     pCM[4] = -pCM[4];
1456
1457     // |M|^2
1458     // XXX - Extra factor of (3) from picking a final state configuration???
1459     // Extra factor of (3 / 8) as averaging over incoming gluon
1460     sigma[i] = -3. * (3. / 8.) * m2Calc();
1461
1462   }
1463
1464 }
1465
1466 //--------------------------------------------------------------------------
1467
1468 // Evaluate |M|^2 - incoming flavour dependence
1469
1470 double Sigma3qg2qqqbarSame::sigmaHat() {
1471   // gq or qg incoming
1472   return (id1 == 21) ? sigma[0] : sigma[1];
1473 }
1474
1475 //--------------------------------------------------------------------------
1476
1477 // Select identity, colour and anticolour.
1478
1479 void Sigma3qg2qqqbarSame::setIdColAcol(){
1480
1481   // Pick outgoing flavour configuration
1482   int idIn = (id1 == 21) ? id2 : id1;
1483
1484   // Outgoing flavours; easiest just to map by hand
1485   switch (config) {
1486   case 0: id3 =  idIn; id4 =  idIn;   id5 = -idIn; break;
1487   case 1: id3 =  idIn; id4 = -idIn;   id5 =  idIn; break;
1488   case 2: id3 =  idIn; id4 =  idIn;   id5 = -idIn; break;
1489   case 3: id3 = -idIn; id4 =  idIn;   id5 =  idIn; break;
1490   case 4: id3 =  idIn; id4 = -idIn;   id5 =  idIn; break;
1491   case 5: id3 = -idIn; id4 =  idIn;   id5 =  idIn; break;
1492   }
1493   setId(id1, id2, id3, id4, id5);
1494
1495   // Temporary solution; start with either
1496   // g q1    -> q1    q2 qbar2
1497   // g qbar1 -> qbar1 qbar2 q2
1498   int cols[5][2];
1499   cols[0][0] = 1; cols[0][1] = 2;
1500   if (idIn > 0) {
1501     cols[1][0] = 3; cols[1][1] = 0;
1502     cols[2][0] = 1; cols[2][1] = 0;
1503     cols[3][0] = 3; cols[3][1] = 0;
1504     cols[4][0] = 0; cols[4][1] = 2;
1505   } else {
1506     cols[1][0] = 0; cols[1][1] = 3;
1507     cols[2][0] = 0; cols[2][1] = 2;
1508     cols[3][0] = 0; cols[3][1] = 3;
1509     cols[4][0] = 1; cols[4][1] = 0;
1510   }
1511   // Swap incoming if q/qbar g instead
1512   if (id2 == 21) {
1513     swap(cols[0][0], cols[1][0]);
1514     swap(cols[0][1], cols[1][1]);
1515   }
1516   // Map final state
1517   int i3 = 0, i4 = 0, i5 = 0;
1518   switch (config) {
1519   case 0: i3 = 2; i4 = 3; i5 = 4; break;
1520   case 1: i3 = 2; i4 = 4; i5 = 3; break;
1521   case 2: i3 = 3; i4 = 2; i5 = 4; break;
1522   case 3: i3 = 4; i4 = 2; i5 = 3; break;
1523   case 4: i3 = 3; i4 = 4; i5 = 2; break;
1524   case 5: i3 = 4; i4 = 3; i5 = 2; break;
1525   }
1526   setColAcol(cols[0][0],  cols[0][1],  cols[1][0],  cols[1][1],
1527              cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1528              cols[i5][0], cols[i5][1]);
1529 }
1530
1531 //==========================================================================
1532
1533 } // end namespace Pythia8