PYTHIA8/pythia8130/src/SigmaQCD.cxx

   1 // SigmaQCD.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 // 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 * Rndm::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 (Rndm::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       = Settings::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 * Rndm::flat() );
 156   mNew  = ParticleDataTable::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 * Rndm::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 * Rndm::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) * Rndm::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 * Rndm::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       = Settings::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 * Rndm::flat() );
 345   mNew  = ParticleDataTable::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 = ParticleDataTable::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 * Rndm::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 = ParticleDataTable::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 } // end namespace Pythia8