PYTHIA8/pythia8145/src/SigmaCompositeness.cxx

   1 // SigmaCompositeness.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 // compositeness simulation classes.
   8
   9 #include "SigmaCompositeness.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // Sigma1qg2qStar class.
  16 // Cross section for q g -> q^* (excited quark state).
  17 // Note: for simplicity decay is assumed isotropic.
  18
  19 //--------------------------------------------------------------------------
  20
  21 // Initialize process.
  22
  23 void Sigma1qg2qStar::initProc() {
  24
  25   // Set up process properties from the chosen quark flavour.
  26   idRes         = 4000000 + idq;
  27   codeSave      = 4000 + idq;
  28   if      (idq == 1) nameSave = "d g -> d^*";
  29   else if (idq == 2) nameSave = "u g -> u^*";
  30   else if (idq == 3) nameSave = "s g -> s^*";
  31   else if (idq == 4) nameSave = "c g -> c^*";
  32   else               nameSave = "b g -> b^*";
  33
  34   // Store q* mass and width for propagator.
  35   mRes          = particleDataPtr->m0(idRes);
  36   GammaRes      = particleDataPtr->mWidth(idRes);
  37   m2Res         = mRes*mRes;
  38   GamMRat       = GammaRes / mRes;
  39
  40   // Locally stored properties and couplings.
  41   Lambda        = settingsPtr->parm("ExcitedFermion:Lambda");
  42   coupFcol      = settingsPtr->parm("ExcitedFermion:coupFcol");
  43
  44   // Set pointer to particle properties and decay table.
  45   qStarPtr      = particleDataPtr->particleDataEntryPtr(idRes);
  46
  47 }
  48
  49 //--------------------------------------------------------------------------
  50
  51 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
  52
  53 void Sigma1qg2qStar::sigmaKin() {
  54
  55   // Incoming width for correct quark.
  56   widthIn  = pow3(mH) * alpS * pow2(coupFcol) / (3. * pow2(Lambda));
  57
  58   // Set up Breit-Wigner.
  59   sigBW    = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
  60
  61 }
  62
  63 //--------------------------------------------------------------------------
  64
  65 // Evaluate sigmaHat(sHat) for specific incoming flavours.
  66
  67 double Sigma1qg2qStar::sigmaHat() {
  68
  69   // Identify whether correct incoming flavours.
  70   int idqNow = (id2 == 21) ? id1 : id2;
  71   if (abs(idqNow) != idq) return 0.;
  72
  73   // Outgoing width and total sigma. Done.
  74   return widthIn * sigBW * qStarPtr->resWidthOpen(idqNow, mH);
  75
  76 }
  77
  78 //--------------------------------------------------------------------------
  79
  80 // Select identity, colour and anticolour.
  81
  82 void Sigma1qg2qStar::setIdColAcol() {
  83
  84   // Flavours.
  85   int idqNow = (id2 == 21) ? id1 : id2;
  86   int idqStar = (idqNow > 0) ? idRes : -idRes;
  87   setId( id1, id2, idqStar);
  88
  89   // Colour flow topology.
  90   if (id1 == idqNow) setColAcol( 1, 0, 2, 1, 2, 0);
  91   else               setColAcol( 2, 1, 1, 0, 2, 0);
  92   if (idqNow < 0) swapColAcol();
  93
  94 }
  95
  96 //--------------------------------------------------------------------------
  97
  98 // Evaluate weight for q* decay angle.
  99
 100 double Sigma1qg2qStar::weightDecay( Event& process, int iResBeg,
 101   int iResEnd) {
 102
 103   // q* should sit in entry 5. Sequential Z/W decay assumed isotropic.
 104   if (iResBeg != 5 || iResEnd != 5) return 1.;
 105
 106   // Sign of asymmetry.
 107   int sideIn    = (process[3].idAbs() < 20) ? 1 : 2;
 108   int sideOut   = (process[6].idAbs() < 20) ? 1 : 2;
 109   double eps    = (sideIn == sideOut) ? 1. : -1.;
 110
 111   // Phase space factors.
 112   double mr1    = pow2(process[6].m()) / sH;
 113   double mr2    = pow2(process[7].m()) / sH;
 114   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 115
 116   // Reconstruct decay angle. Default isotropic decay.
 117   double cosThe = (process[3].p() - process[4].p())
 118     * (process[7].p() - process[6].p()) / (sH * betaf);
 119   double wt     = 1.;
 120   double wtMax  = 1.;
 121
 122   // Decay q* -> q (g/gamma) or q (Z^0/W^+-).
 123   int idBoson   = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs();
 124   if (idBoson == 21 || idBoson == 22) {
 125     wt          = 1. + eps * cosThe;
 126     wtMax       = 2.;
 127   } else if (idBoson == 23 || idBoson == 24) {
 128     double mrB  = (sideOut == 1) ? mr2 : mr1;
 129     double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB);
 130     wt          = 1. + eps * cosThe * ratB;
 131     wtMax       = 1. + ratB;
 132   }
 133
 134   // Done.
 135   return (wt / wtMax);
 136
 137 }
 138
 139 //==========================================================================
 140
 141 // Sigma1lgm2lStar class.
 142 // Cross section for l gamma -> l^* (excited lepton state).
 143 // Note: for simplicity decay is assumed isotropic.
 144
 145 //--------------------------------------------------------------------------
 146
 147 // Initialize process.
 148
 149 void Sigma1lgm2lStar::initProc() {
 150
 151   // Set up process properties from the chosen lepton flavour.
 152   idRes         = 4000000 + idl;
 153   codeSave      = 4000 + idl;
 154   if      (idl == 11) nameSave = "e gamma -> e^*";
 155   else if (idl == 13) nameSave = "mu gamma -> mu^*";
 156   else                nameSave = "tau gamma -> tau^*";
 157
 158   // Store l* mass and width for propagator.
 159   mRes          = particleDataPtr->m0(idRes);
 160   GammaRes      = particleDataPtr->mWidth(idRes);
 161   m2Res         = mRes*mRes;
 162   GamMRat       = GammaRes / mRes;
 163
 164   // Locally stored properties and couplings.
 165   Lambda        = settingsPtr->parm("ExcitedFermion:Lambda");
 166   double coupF  = settingsPtr->parm("ExcitedFermion:coupF");
 167   double coupFp = settingsPtr->parm("ExcitedFermion:coupFprime");
 168   coupChg       = -0.5 * coupF - 0.5 * coupFp;
 169
 170   // Set pointer to particle properties and decay table.
 171   qStarPtr      = particleDataPtr->particleDataEntryPtr(idRes);
 172
 173 }
 174
 175 //--------------------------------------------------------------------------
 176
 177 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 178
 179 void Sigma1lgm2lStar::sigmaKin() {
 180
 181   // Incoming width for correct lepton.
 182   widthIn  = pow3(mH) * alpEM * pow2(coupChg) / pow2(Lambda);
 183
 184   // Set up Breit-Wigner.
 185   sigBW    = M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 186
 187 }
 188
 189 //--------------------------------------------------------------------------
 190
 191 // Evaluate sigmaHat(sHat) for specific incoming flavours.
 192
 193 double Sigma1lgm2lStar::sigmaHat() {
 194
 195   // Identify whether correct incoming flavours.
 196   int idlNow = (id2 == 22) ? id1 : id2;
 197   if (abs(idlNow) != idl) return 0.;
 198
 199   // Outgoing width and total sigma. Done.
 200   return widthIn * sigBW * qStarPtr->resWidthOpen(idlNow, mH);
 201
 202 }
 203
 204 //--------------------------------------------------------------------------
 205
 206 // Select identity, colour and anticolour.
 207
 208 void Sigma1lgm2lStar::setIdColAcol() {
 209
 210   // Flavours.
 211   int idlNow = (id2 == 22) ? id1 : id2;
 212   int idlStar = (idlNow > 0) ? idRes : -idRes;
 213   setId( id1, id2, idlStar);
 214
 215   // No colour flow.
 216   setColAcol( 0, 0, 0, 0, 0, 0);
 217
 218 }
 219
 220 //--------------------------------------------------------------------------
 221
 222 // Evaluate weight for l* decay angle.
 223
 224 double Sigma1lgm2lStar::weightDecay( Event& process, int iResBeg,
 225   int iResEnd) {
 226
 227   // l* should sit in entry 5. Sequential Z/W decay assumed isotropic.
 228   if (iResBeg != 5 || iResEnd != 5) return 1.;
 229
 230   // Sign of asymmetry.
 231   int sideIn    = (process[3].idAbs() < 20) ? 1 : 2;
 232   int sideOut   = (process[6].idAbs() < 20) ? 1 : 2;
 233   double eps    = (sideIn == sideOut) ? 1. : -1.;
 234
 235   // Phase space factors.
 236   double mr1    = pow2(process[6].m()) / sH;
 237   double mr2    = pow2(process[7].m()) / sH;
 238   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 239
 240   // Reconstruct decay angle. Default isotropic decay.
 241   double cosThe = (process[3].p() - process[4].p())
 242     * (process[7].p() - process[6].p()) / (sH * betaf);
 243   double wt     = 1.;
 244   double wtMax  = 1.;
 245
 246   // Decay l* -> l gamma or l (Z^0/W^+-).
 247   int idBoson   = (sideOut == 1) ? process[7].idAbs() : process[6].idAbs();
 248   if (idBoson == 22) {
 249     wt          = 1. + eps * cosThe;
 250     wtMax       = 2.;
 251   } else if (idBoson == 23 || idBoson == 24) {
 252     double mrB  = (sideOut == 1) ? mr2 : mr1;
 253     double ratB = (1. - 0.5 * mrB) / (1 + 0.5 * mrB);
 254     wt          = 1. + eps * cosThe * ratB;
 255     wtMax       = 1. + ratB;
 256   }
 257
 258   // Done.
 259   return (wt / wtMax);
 260
 261 }
 262
 263 //==========================================================================
 264
 265 // Sigma2qq2qStarq class.
 266 // Cross section for q q' -> q^* q' (excited quark state).
 267 // Note: for simplicity decay is assumed isotropic.
 268
 269 //--------------------------------------------------------------------------
 270
 271 // Initialize process.
 272
 273 void Sigma2qq2qStarq::initProc() {
 274
 275   // Set up process properties from the chosen quark flavour.
 276   idRes         = 4000000 + idq;
 277   codeSave      = 4020 + idq;
 278   if      (idq == 1) nameSave = "q q -> d^* q";
 279   else if (idq == 2) nameSave = "q q -> u^* q";
 280   else if (idq == 3) nameSave = "q q -> s^* q";
 281   else if (idq == 4) nameSave = "q q -> c^* q";
 282   else               nameSave = "q q -> b^* q";
 283
 284   // Locally stored properties and couplings.
 285   Lambda        = settingsPtr->parm("ExcitedFermion:Lambda");
 286   preFac        = M_PI / pow4(Lambda);
 287
 288   // Secondary open width fractions.
 289   openFracPos = particleDataPtr->resOpenFrac( idRes);
 290   openFracNeg = particleDataPtr->resOpenFrac(-idRes);
 291
 292 }
 293
 294 //--------------------------------------------------------------------------
 295
 296 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 297
 298 void Sigma2qq2qStarq::sigmaKin() {
 299
 300   // Two possible expressions, for like or unlike sign.
 301   sigmaA = preFac * (1. - s3 / sH);
 302   sigmaB = preFac * (-uH) * (sH + tH) / sH2;
 303
 304 }
 305
 306 //--------------------------------------------------------------------------
 307
 308 // Evaluate sigmaHat(sHat) for specific incoming flavours.
 309
 310 double Sigma2qq2qStarq::sigmaHat() {
 311
 312   // Identify different allowed incoming flavour combinations.
 313   int id1Abs   = abs(id1);
 314   int id2Abs   = abs(id2);
 315   double open1 = (id1 > 0) ? openFracPos : openFracNeg;
 316   double open2 = (id2 > 0) ? openFracPos : openFracNeg;
 317   double sigma = 0.;
 318   if (id1 * id2 > 0) {
 319     if (id1Abs == idq) sigma += (4./3.) * sigmaA * open1;
 320     if (id2Abs == idq) sigma += (4./3.) * sigmaA * open2;
 321   } else if (id1Abs == idq && id2 == -id1)
 322     sigma = (8./3.) * sigmaB * (open1 + open2);
 323   else if (id2 == -id1) sigma = sigmaB * (open1 + open2);
 324   else if (id1Abs == idq) sigma = sigmaB * open1;
 325   else if (id2Abs == idq) sigma = sigmaB * open2;
 326
 327   // Done.
 328   return sigma;
 329
 330 }
 331
 332 //--------------------------------------------------------------------------
 333
 334 // Select identity, colour and anticolour.
 335
 336 void Sigma2qq2qStarq::setIdColAcol() {
 337
 338   // Flavours: either side may have been excited.
 339   double open1 = 0.;
 340   double open2 = 0.;
 341   if (abs(id1) == idq) open1 = (id1 > 0) ? openFracPos : openFracNeg;
 342   if (abs(id2) == idq) open2 = (id2 > 0) ? openFracPos : openFracNeg;
 343   if (open1 == 0. && open2 == 0.) {
 344     open1  = (id1 > 0) ? openFracPos : openFracNeg;
 345     open2  = (id2 > 0) ? openFracPos : openFracNeg;
 346   }
 347   bool excite1 = (open1 > 0.);
 348   if (open1 > 0. && open2 > 0.) excite1
 349     = (rndmPtr->flat() * (open1 + open2) < open1);
 350
 351   // Always excited quark in slot 3 so colour flow flipped or not.
 352   if (excite1) {
 353     id3    = (id1 > 0) ? idq : -idq;
 354     id4    = id2;
 355     if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
 356     else               setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 357     if (id1 < 0) swapColAcol();
 358   } else {
 359     id3    = (id2 > 0) ? idq : -idq;
 360     id4    = id1;
 361     swapTU = true;
 362     if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
 363     else               setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
 364     if (id1 < 0) swapColAcol();
 365   }
 366   setId( id1, id2, id3, id4);
 367
 368 }
 369
 370 //==========================================================================
 371
 372 // Sigma2qqbar2lStarlbar class.
 373 // Cross section for q qbar -> l^* lbar (excited lepton state).
 374 // Note: for simplicity decay is assumed isotropic.
 375
 376 //--------------------------------------------------------------------------
 377
 378 // Initialize process.
 379
 380 void Sigma2qqbar2lStarlbar::initProc() {
 381
 382   // Set up process properties from the chosen lepton flavour.
 383   idRes         = 4000000 + idl;
 384   codeSave      = 4020 + idl;
 385   if      (idl == 11) nameSave = "q qbar -> e^*+- e^-+";
 386   else if (idl == 12) nameSave = "q qbar -> nu_e^* nu_ebar";
 387   else if (idl == 13) nameSave = "q qbar -> mu^*+- mu^-+";
 388   else if (idl == 14) nameSave = "q qbar -> nu_mu^* nu_mubar";
 389   else if (idl == 15) nameSave = "q qbar -> tau^*+- tau^-+";
 390   else                nameSave = "q qbar -> nu_tau^* nu_taubar";
 391
 392   // Secondary open width fractions.
 393   openFracPos = particleDataPtr->resOpenFrac( idRes);
 394   openFracNeg = particleDataPtr->resOpenFrac(-idRes);
 395
 396   // Locally stored properties and couplings.
 397   Lambda        = settingsPtr->parm("ExcitedFermion:Lambda");
 398   preFac        = (M_PI / pow4(Lambda)) * (openFracPos + openFracNeg) / 3.;
 399
 400 }
 401
 402 //--------------------------------------------------------------------------
 403
 404 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 405
 406 void Sigma2qqbar2lStarlbar::sigmaKin() {
 407
 408   // Only one possible expressions
 409   sigma = preFac * (-uH) * (sH + tH) / sH2;
 410
 411 }
 412
 413 //--------------------------------------------------------------------------
 414
 415 // Select identity, colour and anticolour.
 416
 417 void Sigma2qqbar2lStarlbar::setIdColAcol() {
 418
 419   // Flavours: either lepton or antilepton may be excited.
 420   if (rndmPtr->flat() * (openFracPos + openFracNeg) < openFracPos) {
 421     setId( id1, id2, idRes, -idl);
 422     if (id1 < 0) swapTU = true;
 423   } else {
 424     setId( id1, id2, -idRes, idl);
 425     if (id1 > 0) swapTU = true;
 426   }
 427
 428   // Colour flow trivial.
 429   if (id1 > 0) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
 430   else         setColAcol( 0, 1, 1, 0, 0, 0, 0, 0);
 431
 432 }
 433
 434 //==========================================================================
 435
 436 // Sigma2QCqq2qq class.
 437 // Cross section for q q -> q q (quark contact interactions).
 438
 439 //--------------------------------------------------------------------------
 440
 441 // Initialize process.
 442
 443 void Sigma2QCqq2qq::initProc() {
 444
 445   m_Lambda2  = settingsPtr->parm("ContactInteractions:Lambda");
 446   m_etaLL    = settingsPtr->mode("ContactInteractions:etaLL");;
 447   m_etaRR    = settingsPtr->mode("ContactInteractions:etaRR");;
 448   m_etaLR    = settingsPtr->mode("ContactInteractions:etaLR");;
 449   m_Lambda2 *= m_Lambda2;
 450
 451 }
 452
 453 //--------------------------------------------------------------------------
 454
 455 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 456
 457 void Sigma2QCqq2qq::sigmaKin() {
 458
 459   // Calculate kinematics dependence for different terms.
 460   sigT   = (4./9.) * (sH2 + uH2) / tH2;
 461   sigU   = (4./9.) * (sH2 + tH2) / uH2;
 462   sigTU  = - (8./27.) * sH2 / (tH * uH);
 463   sigST  = - (8./27.) * uH2 / (sH * tH);
 464
 465   sigQCSTU = sH2 * (1 / tH + 1 / uH);
 466   sigQCUTS = uH2 * (1 / tH + 1 / sH);
 467
 468 }
 469
 470 //--------------------------------------------------------------------------
 471
 472 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 473
 474 double Sigma2QCqq2qq::sigmaHat() {
 475
 476   // Terms from QC contact interactions.
 477   double sigQCLL = 0;
 478   double sigQCRR = 0;
 479   double sigQCLR = 0;
 480
 481   // Combine cross section terms; factor 1/2 when identical quarks.
 482   // q q -> q q
 483   if (id2 ==  id1) {
 484
 485     sigSum = 0.5 * (sigT + sigU + sigTU); // SM terms.
 486
 487     // Contact terms.
 488     sigQCLL = (8./9.) * alpS * (m_etaLL/m_Lambda2) * sigQCSTU
 489             + (8./3.) * pow2(m_etaLL/m_Lambda2) * sH2;
 490     sigQCRR = (8./9.) * alpS * (m_etaRR/m_Lambda2) * sigQCSTU
 491             + (8./3.) * pow2(m_etaRR/m_Lambda2) * sH2;
 492     sigQCLR = 2. * (uH2 + tH2) * pow2(m_etaLR/m_Lambda2);
 493
 494     sigQCLL /= 2;
 495     sigQCRR /= 2;
 496     sigQCLR /= 2;
 497
 498   // q qbar -> q qbar, without pure s-channel term.
 499   } else if (id2 == -id1) {
 500
 501     sigSum = sigT + sigST; // SM terms.
 502
 503     // Contact terms, minus the terms included in qqbar2qqbar.
 504     sigQCLL = (8./9.) * alpS * (m_etaLL/m_Lambda2) * sigQCUTS
 505             + (5./3.) * pow2(m_etaLL/m_Lambda2) * uH2;
 506     sigQCRR = (8./9.) * alpS * (m_etaRR/m_Lambda2) * sigQCUTS
 507             + (5./3.) * pow2(m_etaRR/m_Lambda2) * uH2;
 508     sigQCLR = 2. * sH2 * pow2(m_etaLR/m_Lambda2);
 509
 510   // q q' -> q q' or q qbar' -> q qbar'
 511   } else {
 512
 513     sigSum = sigT; // SM terms.
 514
 515     // Contact terms.
 516     if (id1 * id2 > 0) {
 517       sigQCLL = pow2(m_etaLL/m_Lambda2) * sH2;
 518       sigQCRR = pow2(m_etaRR/m_Lambda2) * sH2;
 519       sigQCLR = 2 * pow2(m_etaLR/m_Lambda2) * uH2;
 520     } else {
 521       sigQCLL = pow2(m_etaLL/m_Lambda2) * uH2;
 522       sigQCRR = pow2(m_etaRR/m_Lambda2) * uH2;
 523       sigQCLR = 2 * pow2(m_etaLR/m_Lambda2) * sH2;
 524     }
 525   }
 526
 527   // Answer.
 528   double sigma = (M_PI/sH2) * ( pow2(alpS) * sigSum
 529                + sigQCLL + sigQCRR + sigQCLR );
 530   return sigma;
 531
 532 }
 533
 534 //--------------------------------------------------------------------------
 535
 536 // Select identity, colour and anticolour.
 537
 538 void Sigma2QCqq2qq::setIdColAcol() {
 539
 540   // Outgoing = incoming flavours.
 541   setId( id1, id2, id1, id2);
 542
 543   // Colour flow topologies. Swap when antiquarks.
 544   if (id1 * id2 > 0)  setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
 545   else                setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
 546   if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT)
 547                       setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
 548   if (id1 < 0) swapColAcol();
 549
 550 }
 551
 552 //==========================================================================
 553
 554 // Sigma2QCqqbar2qqbar class.
 555 // Cross section for q qbar -> q' qbar' (quark contact interactions).
 556
 557 //--------------------------------------------------------------------------
 558
 559 // Initialize process.
 560
 561 void Sigma2QCqqbar2qqbar::initProc() {
 562
 563   m_nQuarkNew = settingsPtr->mode("ContactInteractions:nQuarkNew");
 564   m_Lambda2   = settingsPtr->parm("ContactInteractions:Lambda");
 565   m_etaLL     = settingsPtr->mode("ContactInteractions:etaLL");
 566   m_etaRR     = settingsPtr->mode("ContactInteractions:etaRR");
 567   m_etaLR     = settingsPtr->mode("ContactInteractions:etaLR");
 568   m_Lambda2  *= m_Lambda2;
 569
 570 }
 571
 572 //--------------------------------------------------------------------------
 573
 574 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
 575
 576 void Sigma2QCqqbar2qqbar::sigmaKin() {
 577
 578   // Pick new flavour.
 579   idNew = 1 + int( m_nQuarkNew * rndmPtr->flat() );
 580   mNew  = particleDataPtr->m0(idNew);
 581   m2New = mNew*mNew;
 582
 583   // Calculate kinematics dependence.
 584   double sigQC              = 0.;
 585   sigS                      = 0.;
 586   if (sH > 4. * m2New) {
 587     sigS = (4./9.) * (tH2 + uH2) / sH2;
 588     sigQC = pow2(m_etaLL/m_Lambda2) * uH2
 589           + pow2(m_etaRR/m_Lambda2) * uH2
 590           + 2 * pow2(m_etaLR/m_Lambda2) * tH2;
 591   }
 592
 593   // Answer is proportional to number of outgoing flavours.
 594   sigma = (M_PI / sH2) * m_nQuarkNew * ( pow2(alpS) * sigS + sigQC);
 595
 596 }
 597
 598 //--------------------------------------------------------------------------
 599
 600 // Select identity, colour and anticolour.
 601
 602 void Sigma2QCqqbar2qqbar::setIdColAcol() {
 603
 604   // Set outgoing flavours ones.
 605   id3 = (id1 > 0) ? idNew : -idNew;
 606   setId( id1, id2, id3, -id3);
 607
 608   // Colour flow topologies. Swap when antiquarks.
 609   setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
 610   if (id1 < 0) swapColAcol();
 611
 612 }
 613
 614 //==========================================================================
 615
 616 } // end namespace Pythia8