PYTHIA8/pythia8130/src/SigmaOnia.cxx

   1 // SigmaOnia.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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 // charmonia/bottomonia simulation classes.
   8
   9 #include "SigmaOnia.h"
  10
  11 namespace Pythia8 {
  12
  13 //**************************************************************************
  14
  15 // Sigma2gg2QQbar3S11g class.
  16 // Cross section g g -> QQbar[3S1(1)] g (Q = c or b).
  17
  18 //*********
  19
  20 // Initialize process.
  21
  22 void Sigma2gg2QQbar3S11g::initProc() {
  23
  24   // Produced state. Process name. Onium matrix element.
  25   idHad = (idNew == 4) ? 443 : 553;
  26   nameSave = (idNew == 4) ? "g g -> ccbar[3S1(1)] g"
  27     : "g g -> bbbar[3S1(1)] g";
  28   oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3S11")
  29     : Settings::parm("Bottomonium:OUpsilon3S11");
  30
  31 }
  32
  33 //*********
  34
  35 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
  36
  37 void Sigma2gg2QQbar3S11g::sigmaKin() {
  38
  39   // Calculate kinematics dependence.
  40   double stH = sH + tH;
  41   double tuH = tH + uH;
  42   double usH = uH + sH;
  43   double sig = (10. * M_PI / 81.) * m3 * ( pow2(sH * tuH)
  44     + pow2(tH * usH) + pow2(uH * stH) ) / pow2( stH * tuH * usH );
  45
  46   // Answer.
  47   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
  48
  49 }
  50
  51 //*********
  52
  53 // Select identity, colour and anticolour.
  54
  55 void Sigma2gg2QQbar3S11g::setIdColAcol() {
  56
  57   // Flavours are trivial.
  58   setId( id1, id2, idHad, 21);
  59
  60   // Two orientations of colour flow..
  61   setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
  62   if (Rndm::flat() > 0.5) swapColAcol();
  63
  64 }
  65
  66 //**************************************************************************
  67
  68 // Sigma2gg2QQbar3PJ1g class.
  69 // Cross section g g -> QQbar[3PJ(1)] g (Q = c or b, J = 0, 1 or 2).
  70
  71 //*********
  72
  73 // Initialize process.
  74
  75 void Sigma2gg2QQbar3PJ1g::initProc() {
  76
  77   // Produced state. Process name. Onium matrix element.
  78   idHad = 0;
  79   nameSave = "illegal process";
  80   if (jSave == 0) {
  81     idHad = (idNew == 4) ? 10441 : 10551;
  82     nameSave = (idNew == 4) ? "g g -> ccbar[3P0(1)] g"
  83       : "g g -> bbbar[3P0(1)] g";
  84   } else if (jSave == 1) {
  85     idHad = (idNew == 4) ? 20443 : 20553;
  86     nameSave = (idNew == 4) ? "g g -> ccbar[3P1(1)] g"
  87       : "g g -> bbbar[3P1(1)] g";
  88   } else if (jSave == 2) {
  89     idHad = (idNew == 4) ? 445 : 555;
  90     nameSave = (idNew == 4) ? "g g -> ccbar[3P2(1)] g"
  91       : "g g -> bbbar[3P2(1)] g";
  92   }
  93   oniumME = (idNew == 4) ? Settings::parm("Charmonium:Ochic03P01")
  94     : Settings::parm("Bottomonium:Ochib03P01");
  95
  96 }
  97
  98 //*********
  99
 100 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
 101
 102 void Sigma2gg2QQbar3PJ1g::sigmaKin() {
 103
 104   // Useful derived kinematics quantities.
 105   double pRat  = (sH * uH + uH * tH + tH * sH)/ sH2;
 106   double qRat  = tH * uH / sH2;
 107   double rRat  = s3 / sH;
 108   double pRat2 = pRat * pRat;
 109   double pRat3 = pRat2 * pRat;
 110   double pRat4 = pRat3 * pRat;
 111   double qRat2 = qRat * qRat;
 112   double qRat3 = qRat2 * qRat;
 113   double qRat4 = qRat3 * qRat;
 114   double rRat2 = rRat * rRat;
 115   double rRat3 = rRat2 * rRat;
 116   double rRat4 = rRat3 * rRat;
 117
 118   // Calculate kinematics dependence.
 119   double sig = 0.;
 120   if (jSave == 0) {
 121     sig = (8. * M_PI / (9. * m3 * sH))
 122       * ( 9. * rRat2 * pRat4 * (rRat4 - 2. * rRat2 * pRat + pRat2)
 123       - 6. * rRat * pRat3 * qRat * (2. * rRat4 - 5. * rRat2 * pRat
 124       + pRat2) - pRat2 * qRat2 * (rRat4 + 2. * rRat2 * pRat - pRat2)
 125       + 2. * rRat * pRat * qRat3 * (rRat2 - pRat) + 6. * rRat2 * qRat4)
 126       / (qRat * pow4(qRat - rRat * pRat));
 127   } else if (jSave == 1) {
 128     sig =  (8. * M_PI / (3.* m3 * sH)) * pRat2
 129       * (rRat * pRat2 * (rRat2 - 4. * pRat)
 130       + 2. * qRat * (-rRat4 + 5. * rRat2 * pRat + pRat2)
 131       - 15. * rRat * qRat2) / pow4(qRat - rRat * pRat);
 132   } else if (jSave == 2) {
 133     sig = (8. * M_PI / (9. * m3 * sH))
 134       * (12. * rRat2 * pRat4 * (rRat4 - 2. * rRat2 * pRat + pRat2)
 135       - 3. * rRat * pRat3 * qRat * (8. * rRat4 - rRat2 * pRat + 4. * pRat2)
 136       + 2. * pRat2 * qRat2 * (-7. * rRat4 + 43. * rRat2 * pRat + pRat2)
 137       + rRat * pRat * qRat3 * (16. * rRat2 - 61. * pRat)
 138       + 12. * rRat2 * qRat4) / (qRat * pow4(qRat-rRat * pRat));
 139   }
 140
 141   // Answer.
 142   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
 143
 144 }
 145
 146 //*********
 147
 148 // Select identity, colour and anticolour.
 149
 150 void Sigma2gg2QQbar3PJ1g::setIdColAcol() {
 151
 152   // Flavours are trivial.
 153   setId( id1, id2, idHad, 21);
 154
 155   // Two orientations of colour flow.
 156   setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
 157   if (Rndm::flat() > 0.5) swapColAcol();
 158
 159 }
 160
 161 //**************************************************************************
 162
 163 // Sigma2qg2QQbar3PJ1q class.
 164 // Cross section q g -> QQbar[3PJ(1)] q (Q = c or b, J = 0, 1 or 2).
 165
 166 //*********
 167
 168 // Initialize process.
 169
 170 void Sigma2qg2QQbar3PJ1q::initProc() {
 171
 172   // Produced state. Process name. Onium matrix element.
 173   idHad = 0;
 174   nameSave = "illegal process";
 175   if (jSave == 0) {
 176     idHad = (idNew == 4) ? 10441 : 10551;
 177     nameSave = (idNew == 4) ? "q g -> ccbar[3P0(1)] q"
 178       : "q g -> bbbar[3P0(1)] q";
 179   } else if (jSave == 1) {
 180     idHad = (idNew == 4) ? 20443 : 20553;
 181     nameSave = (idNew == 4) ? "q g -> ccbar[3P1(1)] q"
 182       : "q g -> bbbar[3P1(1)] q";
 183   } else if (jSave == 2) {
 184    idHad = (idNew == 4) ? 445 : 555;
 185     nameSave = (idNew == 4) ? "q g -> ccbar[3P2(1)] q"
 186       : "q g -> bbbar[3P2(1)] q";
 187   }
 188   oniumME = (idNew == 4) ? Settings::parm("Charmonium:Ochic03P01")
 189     : Settings::parm("Bottomonium:Ochib03P01");
 190
 191 }
 192
 193 //*********
 194
 195 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
 196
 197 void Sigma2qg2QQbar3PJ1q::sigmaKin() {
 198
 199   // Calculate kinematics dependence.
 200   double usH = uH + sH;
 201   double sig = 0.;
 202   if (jSave == 0) {
 203     sig = - (16. * M_PI / 81.) * pow2(tH - 3. * s3) * (sH2 + uH2)
 204       / (m3 * tH * pow4(usH));
 205   } else if (jSave == 1) {
 206     sig = - (32. * M_PI / 27.) * (4. * s3 * sH * uH + tH * (sH2 + uH2))
 207       / (m3 * pow4(usH));
 208   } else if (jSave == 2) {
 209     sig = - (32. *M_PI / 81.) * ( (6. * s3*s3 + tH2) * pow2(usH)
 210       - 2. * sH * uH * (tH2 + 6. * s3 * usH)) / (m3 * tH * pow4(usH));
 211   }
 212
 213   // Answer.
 214   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
 215
 216 }
 217
 218 //*********
 219
 220 // Select identity, colour and anticolour.
 221
 222 void Sigma2qg2QQbar3PJ1q::setIdColAcol() {
 223
 224   // Flavours are trivial.
 225   int idq = (id2 == 21) ? id1 : id2;
 226   setId( id1, id2, idHad, idq);
 227
 228   // tH defined between q_in and q_out: must swap tHat <-> uHat if q g in.
 229   swapTU = (id2 == 21);
 230
 231   // Colour flow topologies. Swap when antiquarks.
 232   if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
 233   else           setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
 234   if (idq < 0) swapColAcol();
 235
 236 }
 237
 238 //**************************************************************************
 239
 240 // Sigma2qqbar2QQbar3PJ1g class.
 241 // Cross section q qbar -> QQbar[3PJ(1)] g (Q = c or b, J = 0, 1 or 2).
 242
 243 //*********
 244
 245 // Initialize process.
 246
 247 void Sigma2qqbar2QQbar3PJ1g::initProc() {
 248
 249   // Produced state. Process name. Onium matrix element.
 250   idHad = 0;
 251   nameSave = "illegal process";
 252   if (jSave == 0) {
 253     idHad = (idNew == 4) ? 10441 : 10551;
 254     nameSave = (idNew == 4) ? "q qbar -> ccbar[3P0(1)] g"
 255       : "q qbar -> bbbar[3P0(1)] g";
 256   } else if (jSave == 1) {
 257     idHad = (idNew == 4) ? 20443 : 20553;
 258     nameSave = (idNew == 4) ? "q qbar -> ccbar[3P1(1)] g"
 259       : "q qbar -> bbbar[3P1(1)] g";
 260   } else if (jSave == 2) {
 261    idHad = (idNew == 4) ? 445 : 555;
 262     nameSave = (idNew == 4) ? "q qbar -> ccbar[3P2(1)] g"
 263       : "q qbar -> bbbar[3P2(1)] g";
 264   }
 265   oniumME = (idNew == 4) ? Settings::parm("Charmonium:Ochic03P01")
 266     : Settings::parm("Bottomonium:Ochib03P01");
 267
 268 }
 269
 270 //*********
 271
 272 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
 273
 274 void Sigma2qqbar2QQbar3PJ1g::sigmaKin() {
 275
 276   // Calculate kinematics dependence.
 277   double tuH = tH + uH;
 278   double sig = 0.;
 279   if (jSave == 0) {
 280     sig =(128. * M_PI / 243.) * pow2(sH - 3. * s3) * (tH2 + uH2)
 281       / (m3 * sH * pow4(tuH));
 282   } else if (jSave == 1) {
 283     sig = (256. * M_PI / 81.) * (4. * s3 * tH * uH + sH * (tH2 + uH2))
 284       / (m3 * pow4(tuH));
 285   } else if (jSave == 2) {
 286     sig = (256. * M_PI / 243.) * ( (6. * s3*s3 + sH2) * pow2(tuH)
 287       - 2. * tH * uH * (sH2 + 6. * s3 * tuH) )/ (m3 * sH * pow4(tuH));
 288   }
 289
 290   // Answer.
 291   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
 292
 293 }
 294
 295 //*********
 296
 297 // Select identity, colour and anticolour.
 298
 299 void Sigma2qqbar2QQbar3PJ1g::setIdColAcol() {
 300
 301   // Flavours are trivial.
 302   setId( id1, id2, idHad, 21);
 303
 304   // Colour flow topologies. Swap when antiquarks.
 305   setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
 306   if (id1 < 0) swapColAcol();
 307
 308 }
 309
 310 //**************************************************************************
 311
 312 // Sigma2gg2QQbarX8g class.
 313 // Cross section g g -> QQbar[X(8)] g (Q = c or b, X = 3S1, 1S0 or 3PJ).
 314
 315 //*********
 316
 317 // Initialize process.
 318
 319 void Sigma2gg2QQbarX8g::initProc() {
 320
 321   // Produced state. Process name. Onium matrix element.
 322   idHad = 0;
 323   nameSave = "illegal process";
 324   if (stateSave == 0) {
 325     idHad = (idNew == 4) ? 9900443 : 9900553;
 326     nameSave = (idNew == 4) ? "g g -> ccbar[3S1(8)] g"
 327       : "g g -> bbbar[3S1(8)] g";
 328     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3S18")
 329       : Settings::parm("Bottomonium:OUpsilon3S18");
 330   } else if (stateSave == 1) {
 331     idHad = (idNew == 4) ? 9900441 : 9900551;
 332     nameSave = (idNew == 4) ? "g g -> ccbar[1S0(8)] g"
 333       : "g g -> bbbar[1S0(8)] g";
 334     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi1S08")
 335       : Settings::parm("Bottomonium:OUpsilon1S08");
 336   } else if (stateSave == 2) {
 337     idHad = (idNew == 4) ? 9910441 : 9910551;
 338     nameSave = (idNew == 4) ? "g g -> ccbar[3PJ(8)] g"
 339       : "g g -> bbbar[3PJ(8)] g";
 340     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3P08")
 341       : Settings::parm("Bottomonium:OUpsilon3P08");
 342   }
 343
 344 }
 345
 346 //*********
 347
 348 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
 349
 350 void Sigma2gg2QQbarX8g::sigmaKin() {
 351
 352   // Calculate kinematics dependence.
 353   double stH = sH + tH;
 354   double tuH = tH + uH;
 355   double usH = uH + sH;
 356   double sig = 0.;
 357   if (stateSave == 0) {
 358     sig = (M_PI / 72.) * m3 * ( 27. * (pow2(stH) + pow2(tuH)
 359       + pow2(usH)) / (s3*s3) - 16. ) * ( pow2(sH * tuH)
 360       + pow2(tH * usH) + pow2(uH * stH) ) / pow2( stH * tuH * usH );
 361   } else if (stateSave == 1) {
 362     sig = (5. * M_PI / 16.) * m3 * ( pow2(uH / (tuH * usH))
 363       + pow2(sH / (stH * usH)) + pow2(tH / (stH * tuH)) ) * ( 12.
 364       + (pow4(stH) + pow4(tuH) + pow4(usH)) / (s3 * sH * tH * uH) );
 365   } else if (stateSave == 2) {
 366     double sH3 = sH2 * sH;
 367     double sH4 = sH3 * sH;
 368     double sH5 = sH4 * sH;
 369     double sH6 = sH5 * sH;
 370     double sH7 = sH6 * sH;
 371     double sH8 = sH7 * sH;
 372     double tH3 = tH2 * tH;
 373     double tH4 = tH3 * tH;
 374     double tH5 = tH4 * tH;
 375     double tH6 = tH5 * tH;
 376     double tH7 = tH6 * tH;
 377     double tH8 = tH7 * tH;
 378     double ssttH = sH * sH + sH * tH + tH * tH;
 379     sig = 5. * M_PI * (3. * sH * tH * stH * pow4(ssttH)
 380       - s3 * pow2(ssttH) * (7. * sH6 + 36. * sH5 * tH + 45. * sH4 * tH2
 381         + 28. * sH3 * tH3 + 45. * sH2 * tH4 + 36. * sH * tH5 + 7. * tH6)
 382       + pow2(s3) * stH * (35. *sH8 + 169. * sH7 * tH + 299. * sH6 * tH2
 383         + 401. * sH5 * tH3 + 418. * sH4 * tH4 + 401. * sH3 * tH5
 384         + 299. * sH2 * tH6 + 169. * sH * tH7 + 35. * tH8)
 385       - pow3(s3) * (84. *sH8+432. *sH7*tH+905. *sH6*tH2
 386         + 1287. * sH5 * tH3 + 1436. * sH4 * tH4 +1287. * sH3 * tH5
 387         + 905. * sH2 * tH6 + 432. * sH * tH7 + 84. * tH8)
 388       + pow4(s3) * stH * (126. * sH6 + 451. * sH5 * tH +677. * sH4 * tH2
 389         + 836. * sH3 * tH3 + 677. * sH2 * tH4 + 451. * sH * tH5
 390         + 126. * tH6)
 391       - pow5(s3) * 3. * (42. * sH6 + 171. * sH5 * tH + 304. * sH4 * tH2
 392         + 362. * sH3 * tH3 + 304. * sH2 * tH4 + 171. * sH * tH5
 393         + 42. * tH6)
 394       + pow3(s3 * s3) * 2. * stH * (42. * sH4 + 106. * sH3 * tH
 395         + 119. * sH2 * tH2 + 106. * sH * tH3 + 42. * tH4)
 396       - pow4(s3) * pow3(s3) * (35. * sH4 + 99. * sH3 * tH
 397         + 120. * sH2 * tH2 + 99.  * sH * tH3 + 35.  * tH4)
 398       + pow4(s3 * s3) * 7. * stH * ssttH)
 399       / (sH * tH * uH * s3 * m3 * pow3(stH * tuH * usH));
 400   }
 401
 402   // Answer.
 403   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
 404
 405 }
 406
 407 //*********
 408
 409 // Select identity, colour and anticolour.
 410
 411 void Sigma2gg2QQbarX8g::setIdColAcol() {
 412
 413   // Flavours are trivial.
 414   setId( id1, id2, idHad, 21);
 415
 416   // Split total contribution into different colour flows just like in
 417   // g g -> g g (with kinematics recalculated for massless partons).
 418   double sHr    = - (tH + uH);
 419   double sH2r   = sHr * sHr;
 420   double sigTS  = tH2/sH2r + 2.*tH/sHr + 3. + 2.*sHr/tH + sH2r/tH2;
 421   double sigUS  = uH2/sH2r + 2.*uH/sHr + 3. + 2.*sHr/uH + sH2r/uH2;
 422   double sigTU  = tH2/uH2 + 2.*tH/uH + 3. + 2.*uH/tH + uH2/tH2;
 423   double sigSum = sigTS + sigUS + sigTU;
 424
 425   // Three colour flow topologies, each with two orientations.
 426   double sigRand = sigSum * Rndm::flat();
 427   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
 428   else if (sigRand < sigTS + sigUS)
 429                        setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
 430   else                 setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
 431   if (Rndm::flat() > 0.5) swapColAcol();
 432
 433
 434 }
 435
 436 //**************************************************************************
 437
 438 // Sigma2qg2QQbarX8q class.
 439 // Cross section q g -> QQbar[X(8)] q (Q = c or b, X = 3S1, 1S0 or 3PJ).
 440
 441 //*********
 442
 443 // Initialize process.
 444
 445 void Sigma2qg2QQbarX8q::initProc() {
 446
 447   // Produced state. Process name. Onium matrix element.
 448   idHad = 0;
 449   nameSave = "illegal process";
 450   if (stateSave == 0) {
 451     idHad = (idNew == 4) ? 9900443 : 9900553;
 452     nameSave = (idNew == 4) ? "q g -> ccbar[3S1(8)] q"
 453       : "q g -> bbbar[3S1(8)] q";
 454     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3S18")
 455       : Settings::parm("Bottomonium:OUpsilon3S18");
 456   } else if (stateSave == 1) {
 457     idHad = (idNew == 4) ? 9900441 : 9900551;
 458     nameSave = (idNew == 4) ? "q g -> ccbar[1S0(8)] q"
 459       : "q g -> bbbar[1S0(8)] q";
 460     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi1S08")
 461       : Settings::parm("Bottomonium:OUpsilon1S08");
 462   } else if (stateSave == 2) {
 463     idHad = (idNew == 4) ? 9910441 : 9910551;
 464     nameSave = (idNew == 4) ? "q g -> ccbar[3PJ(8)] q"
 465       : "q g -> bbbar[3PJ(8)] q";
 466     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3P08")
 467       : Settings::parm("Bottomonium:OUpsilon3P08");
 468   }
 469
 470 }
 471
 472 //*********
 473
 474 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
 475
 476 void Sigma2qg2QQbarX8q::sigmaKin() {
 477
 478   // Calculate kinematics dependence.
 479   double stH  = sH + tH;
 480   double tuH  = tH + uH;
 481   double usH  = uH + sH;
 482   double stH2 = stH * stH;
 483   double tuH2 = tuH * tuH;
 484   double usH2 = usH * usH;
 485   double sig  = 0.;
 486   if (stateSave == 0) {
 487     sig = - (M_PI / 27.)* (4. * (sH2 + uH2) - sH * uH) * (stH2 +tuH2)
 488       / (s3 * m3 * sH * uH * usH2);
 489   } else if (stateSave == 1) {
 490     sig = - (5. * M_PI / 18.) * (sH2 + uH2) / (m3 * tH * usH2);
 491   } else if (stateSave == 2) {
 492     sig = - (10. * M_PI / 9.) * ( (7. * usH + 8. * tH) * (sH2 + uH2)
 493       + 4. * tH * (2. * pow2(s3) - stH2 - tuH2) )
 494       / (s3 * m3 * tH * usH2 * usH);
 495   }
 496
 497   // Answer.
 498   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
 499
 500 }
 501
 502 //*********
 503
 504 // Select identity, colour and anticolour.
 505
 506 void Sigma2qg2QQbarX8q::setIdColAcol() {
 507
 508   // Flavours are trivial.
 509   int idq = (id2 == 21) ? id1 : id2;
 510   setId( id1, id2, idHad, idq);
 511
 512   // tH defined between q_in and q_out: must swap tHat <-> uHat if q g in.
 513   swapTU = (id2 == 21);
 514
 515   // Split total contribution into different colour flows just like in
 516   // q g -> q g (with kinematics recalculated for massless partons).
 517   double sHr    = - (tH + uH);
 518   double sH2r   = sHr * sHr;
 519   double sigTS  = uH2/tH2 - (4./9.) * uH/sHr;
 520   double sigTU  = sH2r/tH2 - (4./9.) * sHr/uH;
 521   double sigSum = sigTS + sigTU;
 522
 523   // Two colour flow topologies. Swap if first is gluon, or when antiquark.
 524   double sigRand = sigSum * Rndm::flat();
 525   if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 2, 3, 3, 0);
 526   else                 setColAcol( 1, 0, 2, 3, 1, 3, 2, 0);
 527   if (id1 == 21) swapCol12();
 528   if (idq < 0) swapColAcol();
 529
 530 }
 531
 532 //**************************************************************************
 533
 534 // Sigma2qqbar2QQbarX8g class.
 535 // Cross section q qbar -> QQbar[X(8)] g (Q = c or b, X = 3S1, 1S0 or 3PJ).
 536
 537 //*********
 538
 539 // Initialize process.
 540
 541 void Sigma2qqbar2QQbarX8g::initProc() {
 542
 543   // Produced state. Process name. Onium matrix element.
 544   idHad = 0;
 545   nameSave = "illegal process";
 546   if (stateSave == 0) {
 547     idHad = (idNew == 4) ? 9900443 : 9900553;
 548     nameSave = (idNew == 4) ? "q qbar -> ccbar[3S1(8)] g"
 549       : "q qbar -> bbbar[3S1(8)] g";
 550     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3S18")
 551       : Settings::parm("Bottomonium:OUpsilon3S18");
 552   } else if (stateSave == 1) {
 553     idHad = (idNew == 4) ? 9900441 : 9900551;
 554     nameSave = (idNew == 4) ? "q qbar -> ccbar[1S0(8)] g"
 555       : "q qbar -> bbbar[1S0(8)] g";
 556     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi1S08")
 557       : Settings::parm("Bottomonium:OUpsilon1S08");
 558   } else if (stateSave == 2) {
 559     idHad = (idNew == 4) ? 9910441 : 9910551;
 560     nameSave = (idNew == 4) ? "q qbar -> ccbar[3PJ(8)] g"
 561       : "q qbar -> bbbar[3PJ(8)] g";
 562     oniumME = (idNew == 4) ? Settings::parm("Charmonium:OJpsi3P08")
 563       : Settings::parm("Bottomonium:OUpsilon3P08");
 564   }
 565
 566 }
 567
 568 //*********
 569
 570 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
 571
 572 void Sigma2qqbar2QQbarX8g::sigmaKin() {
 573
 574   // Calculate kinematics dependence.
 575   double stH  = sH + tH;
 576   double tuH  = tH + uH;
 577   double usH  = uH + sH;
 578   double stH2 = stH * stH;
 579   double tuH2 = tuH * tuH;
 580   double usH2 = usH * usH;
 581   double sig  = 0.;
 582   if (stateSave == 0) {
 583     sig = (8. * M_PI / 81.) * (4. * (tH2 + uH2) - tH * uH)
 584       * (stH2 + usH2) / (s3 * m3 * tH * uH * tuH2);
 585   } else if (stateSave == 1) {
 586     sig = (20. * M_PI / 27.) * (tH2 + uH2) / (m3 * sH * tuH2);
 587   } else if (stateSave == 2) {
 588     sig = (80. * M_PI / 27.) * ( (7. * tuH + 8. * sH) * (tH2 + uH2)
 589       + 4. * sH * (2. * pow2(s3) - stH2 -usH2) )
 590       / (s3 * m3 * sH * tuH2 * tuH);
 591   }
 592
 593   // Answer.
 594   sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
 595
 596 }
 597
 598 //*********
 599
 600 // Select identity, colour and anticolour.
 601
 602 void Sigma2qqbar2QQbarX8g::setIdColAcol() {
 603
 604   // Flavours are trivial.
 605   setId( id1, id2, idHad, 21);
 606
 607   // Split total contribution into different colour flows just like in
 608   // q qbar -> g g (with kinematics recalculated for massless partons).
 609   double sHr    = - (tH + uH);
 610   double sH2r   = sHr * sHr;
 611   double sigTS  = (4. / 9.) * uH / tH - uH2 / sH2r;
 612   double sigUS  = (4. / 9.) * tH / uH - tH2 / sH2r;
 613   double sigSum = sigTS + sigUS;
 614
 615   // Two colour flow topologies. Swap if first is antiquark.
 616   double sigRand = sigSum * Rndm::flat();
 617   if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
 618   else                 setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
 619   if (id1 < 0) swapColAcol();
 620
 621 }
 622
 623 //**************************************************************************
 624
 625 } // end namespace Pythia8
 626