1 // SigmaQCD.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2013 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.
6 // Function definitions (not found in the header) for the
7 // QCD simulation classes.
13 //==========================================================================
16 // Cross section for elastic scattering A B -> A B.
18 //--------------------------------------------------------------------------
20 // Select identity, colour and anticolour.
22 void Sigma0AB2AB::setIdColAcol() {
24 // Flavours and colours are trivial.
25 setId( idA, idB, idA, idB);
26 setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
29 //==========================================================================
32 // Cross section for single diffractive scattering A B -> X B.
34 //--------------------------------------------------------------------------
36 // Select identity, colour and anticolour.
38 void Sigma0AB2XB::setIdColAcol() {
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);
48 //==========================================================================
51 // Cross section for single diffractive scattering A B -> A X.
53 //--------------------------------------------------------------------------
55 // Select identity, colour and anticolour.
57 void Sigma0AB2AX::setIdColAcol() {
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);
67 //==========================================================================
70 // Cross section for double diffractive scattering A B -> X X.
72 //--------------------------------------------------------------------------
74 // Select identity, colour and anticolour.
76 void Sigma0AB2XX::setIdColAcol() {
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);
88 //==========================================================================
90 // Sigma0AB2AXB class.
91 // Cross section for central scattering A B -> A X B.
93 //--------------------------------------------------------------------------
95 // Select identity, colour and anticolour.
97 void Sigma0AB2AXB::setIdColAcol() {
99 // Central diffractive state represented by rho_diffr0. Colours trivial.
101 setId( idA, idB, idA, idB,idX);
102 setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
106 //==========================================================================
108 // Sigma2gg2gg class.
109 // Cross section for g g -> g g.
111 //--------------------------------------------------------------------------
113 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
115 void Sigma2gg2gg::sigmaKin() {
117 // Calculate kinematics dependence.
118 sigTS = (9./4.) * (tH2 / sH2 + 2. * tH / sH + 3. + 2. * sH / tH
120 sigUS = (9./4.) * (uH2 / sH2 + 2. * uH / sH + 3. + 2. * sH / uH
122 sigTU = (9./4.) * (tH2 / uH2 + 2. * tH / uH + 3. + 2. * uH / tH
124 sigSum = sigTS + sigUS + sigTU;
126 // Answer contains factor 1/2 from identical gluons.
127 sigma = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
131 //--------------------------------------------------------------------------
133 // Select identity, colour and anticolour.
135 void Sigma2gg2gg::setIdColAcol() {
137 // Flavours are trivial.
138 setId( id1, id2, 21, 21);
140 // Three colour flow topologies, each with two orientations.
141 double sigRand = sigSum * rndmPtr->flat();
142 if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
143 else if (sigRand < sigTS + sigUS)
144 setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
145 else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
146 if (rndmPtr->flat() > 0.5) swapColAcol();
150 //==========================================================================
152 // Sigma2gg2qqbar class.
153 // Cross section for g g -> q qbar (q = u, d, s, i.e. almost massless).
155 //--------------------------------------------------------------------------
157 // Initialize process.
159 void Sigma2gg2qqbar::initProc() {
161 // Read number of quarks to be considered in massless approximation.
162 nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
166 //--------------------------------------------------------------------------
168 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
170 void Sigma2gg2qqbar::sigmaKin() {
173 idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
174 mNew = particleDataPtr->m0(idNew);
177 // Calculate kinematics dependence.
180 if (sH > 4. * m2New) {
181 sigTS = (1./6.) * uH / tH - (3./8.) * uH2 / sH2;
182 sigUS = (1./6.) * tH / uH - (3./8.) * tH2 / sH2;
184 sigSum = sigTS + sigUS;
186 // Answer is proportional to number of outgoing flavours.
187 sigma = (M_PI / sH2) * pow2(alpS) * nQuarkNew * sigSum;
191 //--------------------------------------------------------------------------
193 // Select identity, colour and anticolour.
195 void Sigma2gg2qqbar::setIdColAcol() {
197 // Flavours are trivial.
198 setId( id1, id2, idNew, -idNew);
200 // Two colour flow topologies.
201 double sigRand = sigSum * rndmPtr->flat();
202 if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
203 else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
207 //==========================================================================
209 // Sigma2qg2qg class.
210 // Cross section for q g -> q g.
212 //--------------------------------------------------------------------------
214 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
216 void Sigma2qg2qg::sigmaKin() {
218 // Calculate kinematics dependence.
219 sigTS = uH2 / tH2 - (4./9.) * uH / sH;
220 sigTU = sH2 / tH2 - (4./9.) * sH / uH;
221 sigSum = sigTS + sigTU;
224 sigma = (M_PI / sH2) * pow2(alpS) * sigSum;
228 //--------------------------------------------------------------------------
230 // Select identity, colour and anticolour.
232 void Sigma2qg2qg::setIdColAcol() {
234 // Outgoing = incoming flavours.
235 setId( id1, id2, id1, id2);
237 // Two colour flow topologies. Swap if first is gluon, or when antiquark.
238 double sigRand = sigSum * rndmPtr->flat();
239 if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 3, 0, 2, 3);
240 else setColAcol( 1, 0, 2, 3, 2, 0, 1, 3);
241 if (id1 == 21) swapCol1234();
242 if (id1 < 0 || id2 < 0) swapColAcol();
246 //==========================================================================
248 // Sigma2qq2qq class.
249 // Cross section for q qbar' -> q qbar' or q q' -> q q'
250 // (qbar qbar' -> qbar qbar'), q' may be same as q.
252 //--------------------------------------------------------------------------
254 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
256 void Sigma2qq2qq::sigmaKin() {
258 // Calculate kinematics dependence for different terms.
259 sigT = (4./9.) * (sH2 + uH2) / tH2;
260 sigU = (4./9.) * (sH2 + tH2) / uH2;
261 sigTU = - (8./27.) * sH2 / (tH * uH);
262 sigST = - (8./27.) * uH2 / (sH * tH);
266 //--------------------------------------------------------------------------
269 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
271 double Sigma2qq2qq::sigmaHat() {
273 // Combine cross section terms; factor 1/2 when identical quarks.
274 if (id2 == id1) sigSum = 0.5 * (sigT + sigU + sigTU);
275 else if (id2 == -id1) sigSum = sigT + sigST;
279 return (M_PI/sH2) * pow2(alpS) * sigSum;
283 //--------------------------------------------------------------------------
285 // Select identity, colour and anticolour.
287 void Sigma2qq2qq::setIdColAcol() {
289 // Outgoing = incoming flavours.
290 setId( id1, id2, id1, id2);
292 // Colour flow topologies. Swap when antiquarks.
293 if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
294 else setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
295 if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT)
296 setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
297 if (id1 < 0) swapColAcol();
301 //==========================================================================
303 // Sigma2qqbar2gg class.
304 // Cross section for q qbar -> g g.
306 //--------------------------------------------------------------------------
308 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
310 void Sigma2qqbar2gg::sigmaKin() {
312 // Calculate kinematics dependence.
313 sigTS = (32./27.) * uH / tH - (8./3.) * uH2 / sH2;
314 sigUS = (32./27.) * tH / uH - (8./3.) * tH2 / sH2;
315 sigSum = sigTS + sigUS;
317 // Answer contains factor 1/2 from identical gluons.
318 sigma = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
322 //--------------------------------------------------------------------------
324 // Select identity, colour and anticolour.
326 void Sigma2qqbar2gg::setIdColAcol() {
328 // Outgoing flavours trivial.
329 setId( id1, id2, 21, 21);
331 // Two colour flow topologies. Swap if first is antiquark.
332 double sigRand = sigSum * rndmPtr->flat();
333 if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
334 else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
335 if (id1 < 0) swapColAcol();
339 //==========================================================================
341 // Sigma2qqbar2qqbarNew class.
342 // Cross section q qbar -> q' qbar'.
344 //--------------------------------------------------------------------------
346 // Initialize process.
348 void Sigma2qqbar2qqbarNew::initProc() {
350 // Read number of quarks to be considered in massless approximation.
351 nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
355 //--------------------------------------------------------------------------
357 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
359 void Sigma2qqbar2qqbarNew::sigmaKin() {
362 idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
363 mNew = particleDataPtr->m0(idNew);
366 // Calculate kinematics dependence.
368 if (sH > 4. * m2New) sigS = (4./9.) * (tH2 + uH2) / sH2;
370 // Answer is proportional to number of outgoing flavours.
371 sigma = (M_PI / sH2) * pow2(alpS) * nQuarkNew * sigS;
375 //--------------------------------------------------------------------------
377 // Select identity, colour and anticolour.
379 void Sigma2qqbar2qqbarNew::setIdColAcol() {
381 // Set outgoing flavours ones.
382 id3 = (id1 > 0) ? idNew : -idNew;
383 setId( id1, id2, id3, -id3);
385 // Colour flow topologies. Swap when antiquarks.
386 setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
387 if (id1 < 0) swapColAcol();
391 //==========================================================================
393 // Sigma2gg2QQbar class.
394 // Cross section g g -> Q Qbar (Q = c, b or t).
395 // Only provided for fixed m3 = m4 so do some gymnastics:
396 // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
397 // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
398 // but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
399 // tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
401 //--------------------------------------------------------------------------
403 // Initialize process.
405 void Sigma2gg2QQbar::initProc() {
408 nameSave = "g g -> Q Qbar";
409 if (idNew == 4) nameSave = "g g -> c cbar";
410 if (idNew == 5) nameSave = "g g -> b bbar";
411 if (idNew == 6) nameSave = "g g -> t tbar";
412 if (idNew == 7) nameSave = "g g -> b' b'bar";
413 if (idNew == 8) nameSave = "g g -> t' t'bar";
415 // Secondary open width fraction.
416 openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
420 //--------------------------------------------------------------------------
422 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
424 void Sigma2gg2QQbar::sigmaKin() {
426 // Modified Mandelstam variables for massive kinematics with m3 = m4.
427 double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
428 double tHQ = -0.5 * (sH - tH + uH);
429 double uHQ = -0.5 * (sH + tH - uH);
430 double tHQ2 = tHQ * tHQ;
431 double uHQ2 = uHQ * uHQ;
433 // Calculate kinematics dependence.
434 double tumHQ = tHQ * uHQ - s34Avg * sH;
435 sigTS = ( uHQ / tHQ - 2.25 * uHQ2 / sH2 + 4.5 * s34Avg * tumHQ
436 / ( sH * tHQ2) + 0.5 * s34Avg * (tHQ + s34Avg) / tHQ2
437 - s34Avg*s34Avg / (sH * tHQ) ) / 6.;
438 sigUS = ( tHQ / uHQ - 2.25 * tHQ2 / sH2 + 4.5 * s34Avg * tumHQ
439 / ( sH * uHQ2) + 0.5 * s34Avg * (uHQ + s34Avg) / uHQ2
440 - s34Avg*s34Avg / (sH * uHQ) ) / 6.;
441 sigSum = sigTS + sigUS;
444 sigma = (M_PI / sH2) * pow2(alpS) * sigSum * openFracPair;
448 //--------------------------------------------------------------------------
450 // Select identity, colour and anticolour.
452 void Sigma2gg2QQbar::setIdColAcol() {
454 // Flavours are trivial.
455 setId( id1, id2, idNew, -idNew);
457 // Two colour flow topologies.
458 double sigRand = sigSum * rndmPtr->flat();
459 if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
460 else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
464 //--------------------------------------------------------------------------
466 // Evaluate weight for decay angles of W in top decay.
468 double Sigma2gg2QQbar::weightDecay( Event& process, int iResBeg,
471 // For top decay hand over to standard routine, else done.
472 if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
473 return weightTopDecay( process, iResBeg, iResEnd);
478 //==========================================================================
480 // Sigma2qqbar2QQbar class.
481 // Cross section q qbar -> Q Qbar (Q = c, b or t).
482 // Only provided for fixed m3 = m4 so do some gymnastics:
483 // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
484 // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
485 // but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
486 // tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
488 //--------------------------------------------------------------------------
490 // Initialize process, especially parton-flux object.
492 void Sigma2qqbar2QQbar::initProc() {
495 nameSave = "q qbar -> Q Qbar";
496 if (idNew == 4) nameSave = "q qbar -> c cbar";
497 if (idNew == 5) nameSave = "q qbar -> b bbar";
498 if (idNew == 6) nameSave = "q qbar -> t tbar";
499 if (idNew == 7) nameSave = "q qbar -> b' b'bar";
500 if (idNew == 8) nameSave = "q qbar -> t' t'bar";
502 // Secondary open width fraction.
503 openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
507 //--------------------------------------------------------------------------
509 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
511 void Sigma2qqbar2QQbar::sigmaKin() {
513 // Modified Mandelstam variables for massive kinematics with m3 = m4.
514 double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
515 double tHQ = -0.5 * (sH - tH + uH);
516 double uHQ = -0.5 * (sH + tH - uH);
517 double tHQ2 = tHQ * tHQ;
518 double uHQ2 = uHQ * uHQ;
520 // Calculate kinematics dependence.
521 double sigS = (4./9.) * ((tHQ2 + uHQ2) / sH2 + 2. * s34Avg / sH);
524 sigma = (M_PI / sH2) * pow2(alpS) * sigS * openFracPair;
528 //--------------------------------------------------------------------------
530 // Select identity, colour and anticolour.
532 void Sigma2qqbar2QQbar::setIdColAcol() {
534 // Set outgoing flavours.
535 id3 = (id1 > 0) ? idNew : -idNew;
536 setId( id1, id2, id3, -id3);
538 // Colour flow topologies. Swap when antiquarks.
539 setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
540 if (id1 < 0) swapColAcol();
544 //--------------------------------------------------------------------------
546 // Evaluate weight for decay angles of W in top decay.
548 double Sigma2qqbar2QQbar::weightDecay( Event& process, int iResBeg,
551 // For top decay hand over to standard routine, else done.
552 if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
553 return weightTopDecay( process, iResBeg, iResEnd);
559 //==========================================================================
561 // Sigma3gg2ggg class.
562 // Cross section for g g -> g g g.
564 //--------------------------------------------------------------------------
566 // Evaluate |M|^2 - no incoming flavour dependence.
568 void Sigma3gg2ggg::sigmaKin() {
570 // Calculate all four-vector products.
571 Vec4 p1cm( 0., 0., 0.5 * mH, 0.5 * mH);
572 Vec4 p2cm( 0., 0., -0.5 * mH, 0.5 * mH);
573 pp[1][2] = p1cm * p2cm;
574 pp[1][3] = p1cm * p3cm;
575 pp[1][4] = p1cm * p4cm;
576 pp[1][5] = p1cm * p5cm;
577 pp[2][3] = p2cm * p3cm;
578 pp[2][4] = p2cm * p4cm;
579 pp[2][5] = p2cm * p5cm;
580 pp[3][4] = p3cm * p4cm;
581 pp[3][5] = p3cm * p5cm;
582 pp[4][5] = p4cm * p5cm;
583 for (int i = 1; i < 5; ++i)
584 for (int j = i + 1; j < 6; ++j) pp[j][i] = pp[i][j];
586 // Cross section, in three main sections.
587 double num1 = cycle(1,2,3,4,5) + cycle(1,2,3,5,4) + cycle(1,2,4,3,5)
588 + cycle(1,2,4,5,3) + cycle(1,2,5,3,4) + cycle(1,2,5,4,3)
589 + cycle(1,3,2,4,5) + cycle(1,3,2,5,4) + cycle(1,3,4,2,5)
590 + cycle(1,3,5,2,4) + cycle(1,4,2,3,5) + cycle(1,4,3,2,5);
591 double num2 = pow4(pp[1][2]) + pow4(pp[1][3]) + pow4(pp[1][4])
592 + pow4(pp[1][5]) + pow4(pp[2][3]) + pow4(pp[2][4])
593 + pow4(pp[2][5]) + pow4(pp[3][4]) + pow4(pp[3][5])
595 double den = pp[1][2] * pp[1][3] * pp[1][4] * pp[1][5] * pp[2][3]
596 * pp[2][4] * pp[2][5] * pp[3][4] * pp[3][5] * pp[4][5];
598 // Answer has a factor 6 due to identical gluons
599 // This is cancelled by phase space factor (1 / 6)
600 sigma = pow3(4. * M_PI * alpS) * (27./16.) * num1 * num2 / den;
604 //--------------------------------------------------------------------------
606 // Select identity, colour and anticolour.
608 void Sigma3gg2ggg::setIdColAcol() {
610 // Flavours are trivial.
611 setId( id1, id2, 21, 21, 21);
613 // Three colour flow topologies, each with two orientations.
615 double sigRand = sigSum * rndmPtr->flat();
616 if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
617 else if (sigRand < sigTS + sigUS)
618 setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
619 else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
620 if (rndmPtr->flat() > 0.5) swapColAcol();
623 // Temporary solution.
624 setColAcol( 1, 2, 2, 3, 1, 4, 4, 5, 5, 3);
628 //==========================================================================
630 // Sigma3qqbar2ggg class.
631 // Cross section for q qbar -> g g g.
633 //--------------------------------------------------------------------------
635 // Evaluate |M|^2 - no incoming flavour dependence.
636 void Sigma3qqbar2ggg::sigmaKin() {
638 // Setup four-vectors
639 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
640 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
646 // Answer has a factor 6 due to identical gluons,
647 // which is cancelled by phase space factor (1 / 6)
652 //--------------------------------------------------------------------------
656 inline double Sigma3qqbar2ggg::m2Calc() {
658 // Calculate four-products
659 double sHnow = (pCM[0] + pCM[1]).m2Calc();
660 double sHhalf = sH / 2.;
662 // qbar (p+) + q(p-) -> g(k1) g(k2) g(k3)
663 // a_i = (p+ . ki), i = 1, 2, 3
664 // b_i = (p- . ki), i = 1, 2, 3
665 a[0] = pCM[0] * pCM[2];
666 a[1] = pCM[0] * pCM[3];
667 a[2] = pCM[0] * pCM[4];
668 b[0] = pCM[1] * pCM[2];
669 b[1] = pCM[1] * pCM[3];
670 b[2] = pCM[1] * pCM[4];
672 pp[0][1] = pCM[2] * pCM[3];
673 pp[1][2] = pCM[3] * pCM[4];
674 pp[2][0] = pCM[4] * pCM[2];
676 // ab[i][j] = a_i * b_j + a_j * b_i
677 ab[0][1] = a[0] * b[1] + a[1] * b[0];
678 ab[1][2] = a[1] * b[2] + a[2] * b[1];
679 ab[2][0] = a[2] * b[0] + a[0] * b[2];
682 double num1 = a[0] * b[0] * (a[0] * a[0] + b[0] * b[0]) +
683 a[1] * b[1] * (a[1] * a[1] + b[1] * b[1]) +
684 a[2] * b[2] * (a[2] * a[2] + b[2] * b[2]);
685 double den1 = a[0] * a[1] * a[2] * b[0] * b[1] * b[2];
686 double num2 = - ( ab[0][1] / pp[0][1] )
687 - ( ab[1][2] / pp[1][2] )
688 - ( ab[2][0] / pp[2][0] );
689 double num3 = a[2] * b[2] * ab[0][1] / (pp[1][2] * pp[2][0] ) +
690 a[0] * b[0] * ab[1][2] / (pp[2][0] * pp[0][1] ) +
691 a[1] * b[1] * ab[2][0] / (pp[0][1] * pp[1][2] );
694 return pow3(4. * M_PI * alpS) * (8. / 324.) * (num1 / den1) *
695 ( sHhalf + 9. * (sHhalf + num2) + (2. * 81. / sHnow) * num3 );
699 //--------------------------------------------------------------------------
701 // Select identity, colour and anticolour.
703 void Sigma3qqbar2ggg::setIdColAcol(){
705 // Flavours are trivial.
706 setId( id1, id2, 21, 21, 21);
708 // Temporary solution.
709 setColAcol( 1, 0, 0, 2, 1, 3, 3, 4, 4, 2);
710 if (id1 < 0) swapColAcol();
713 //--------------------------------------------------------------------------
715 // Map a final state configuration
717 inline void Sigma3qqbar2ggg::mapFinal() {
719 case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
720 case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
721 case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
722 case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
723 case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
724 case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
728 //==========================================================================
730 // Sigma3qg2qgg class.
731 // Cross section for q g -> q g g.
732 // Crossed relation from q qbar -> g g g:
733 // qbar(p+) q(p-) -> g(k1) g(k2) g(k3)
735 //--------------------------------------------------------------------------
737 // Evaluate |M|^2 - no incoming flavour dependence
738 // Note: two different contributions from gq and qg incoming
740 void Sigma3qg2qgg::sigmaKin() {
742 // Pick a final state configuration
745 // gq and qg incoming
746 for (int i = 0; i < 2; i++) {
748 // Map incoming four-vectors to p+, p-, k1, k2, k3
749 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
750 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
754 swap(pCM[i], pCM[2]);
757 // XXX - Extra factor of (3) from selecting a final state
758 // configuration (already a factor of 2 in the original
759 // answer due to two identical final state gluons)???
760 // Extra factor of (3 / 8) as average over incoming gluon
761 sigma[i] = 3. * (3. / 8.) * m2Calc();
763 } // for (int i = 0; i < 2; i++)
767 //--------------------------------------------------------------------------
769 // Evaluate |M|^2 - incoming flavour dependence
770 // Pick from two configurations calculated previously
772 double Sigma3qg2qgg::sigmaHat() {
774 return (id1 == 21) ? sigma[0] : sigma[1];
777 //--------------------------------------------------------------------------
779 // Select identity, colour and anticolour.
781 void Sigma3qg2qgg::setIdColAcol(){
782 // Outgoing flavours; only need to know where the quark is
783 int qIdx = config / 2;
784 int idTmp[3] = { 21, 21, 21 };
785 idTmp[qIdx] = (id1 == 21) ? id2 : id1;
786 setId( id1, id2, idTmp[0], idTmp[1], idTmp[2]);
788 // Temporary solution
789 if (qIdx == 0) setColAcol(1, 0, 2, 1, 4, 0, 3, 4, 2, 3);
790 else if (qIdx == 1) setColAcol(1, 0, 2, 1, 3, 4, 4, 0, 2, 3);
791 else setColAcol(1, 0, 2, 1, 3, 4, 2, 3, 4, 0);
794 swap( colSave[1], colSave[2]);
795 swap(acolSave[1], acolSave[2]);
797 // qbar rather than q incoming
798 if (id1 < 0 || id2 < 0) swapColAcol();
802 //==========================================================================
804 // Sigma3gg2qqbarg class.
805 // Cross section for g g -> q qbar g
806 // Crossed relation from q qbar -> g g g
808 //--------------------------------------------------------------------------
810 // Initialize process.
812 void Sigma3gg2qqbarg::initProc() {
814 // Read number of quarks to be considered in massless approximation.
815 nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
819 //--------------------------------------------------------------------------
821 // Evaluate |M|^2 - no incoming flavour dependence.
823 void Sigma3gg2qqbarg::sigmaKin() {
825 // Incoming four-vectors
826 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
827 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
829 // Pick and map a final state configuration
834 swap(pCM[0], pCM[2]);
835 swap(pCM[1], pCM[3]);
838 // Extra factor of (6.) from picking a final state configuration
839 // Extra factor of nQuarkNew
840 // Extra factor of (3. / 8.) ^ 2 as averaging over two incoming gluons
841 sigma = 6. * nQuarkNew * (3. / 8.) * (3. / 8.) * m2Calc();
845 //--------------------------------------------------------------------------
847 // Select identity, colour and anticolour.
849 void Sigma3gg2qqbarg::setIdColAcol(){
852 int idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
854 // Outgoing flavours; easiest just to map by hand
856 case 0: id3 = idNew; id4 = -idNew; id5 = 21; break;
857 case 1: id3 = idNew; id4 = 21; id5 = -idNew; break;
858 case 2: id3 = -idNew; id4 = idNew; id5 = 21; break;
859 case 3: id3 = 21; id4 = idNew; id5 = -idNew; break;
860 case 4: id3 = -idNew; id4 = 21; id5 = idNew; break;
861 case 5: id3 = 21; id4 = -idNew; id5 = idNew; break;
863 setId(id1, id2, id3, id4, id5);
865 // Temporary solution
867 case 0: setColAcol( 1, 2, 2, 3, 4, 0, 0, 3, 1, 4 ); break;
868 case 1: setColAcol( 1, 2, 2, 3, 4, 0, 1, 4, 0, 3 ); break;
869 case 2: setColAcol( 1, 2, 2, 3, 0, 3, 4, 0, 1, 4 ); break;
870 case 3: setColAcol( 1, 2, 2, 3, 1, 4, 4, 0, 0, 3 ); break;
871 case 4: setColAcol( 1, 2, 2, 3, 0, 3, 1, 4, 4, 0 ); break;
872 case 5: setColAcol( 1, 2, 2, 3, 1, 4, 0, 3, 4, 0 ); break;
877 //==========================================================================
879 // Sigma3qq2qqgDiff class.
880 // Cross section for q q' -> q q' g, q != q'
882 //--------------------------------------------------------------------------
884 // Evaluate |M|^2 - no incoming flavour dependence
886 void Sigma3qq2qqgDiff::sigmaKin() {
888 // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
890 // Incoming four-vectors
891 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
892 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
893 // Pick and map a final state configuration
898 // Extra factor of (6.) from picking a final state configuration
899 sigma = 6. * m2Calc();
902 //--------------------------------------------------------------------------
906 inline double Sigma3qq2qqgDiff::m2Calc() {
909 s = (pCM[0] + pCM[1]).m2Calc();
910 t = (pCM[0] - pCM[2]).m2Calc();
911 u = (pCM[0] - pCM[3]).m2Calc();
912 up = (pCM[1] - pCM[2]).m2Calc();
913 sp = (pCM[2] + pCM[3]).m2Calc();
914 tp = (pCM[1] - pCM[3]).m2Calc();
917 double num1 = (s * s + sp * sp + u * u + up * up) / (t * tp);
918 double den1 = (pCM[0] * pCM[4]) * (pCM[1] * pCM[4]) *
919 (pCM[2] * pCM[4]) * (pCM[3] * pCM[4]);
920 double num2 = (u + up) * (s * sp + t * tp - u * up) +
921 u * (s * t + sp * tp) + up * (s * tp + sp * t);
922 double num3 = (s + sp) * (s * sp - t * tp - u * up) +
923 2. * t * tp * (u + up) + 2. * u * up * (t + tp);
925 // (N^2 - 1)^2 / 4N^3 = 16. / 27.
926 // (N^2 - 1) / 4N^3 = 2. / 27.
927 return (1. / 8.) * pow3(4. * M_PI * alpS) * num1 / den1 *
928 ( (16. / 27.) * num2 - (2. / 27.) * num3 );
932 //--------------------------------------------------------------------------
934 // Evaluate |M|^2 - incoming flavour dependence
936 double Sigma3qq2qqgDiff::sigmaHat() {
937 // Different incoming flavours only
938 if (abs(id1) == abs(id2)) return 0.;
942 //--------------------------------------------------------------------------
944 // Select identity, colour and anticolour.
946 void Sigma3qq2qqgDiff::setIdColAcol(){
948 // Outgoing flavours; easiest just to map by hand
950 case 0: id3 = id1; id4 = id2; id5 = 21; break;
951 case 1: id3 = id1; id4 = 21; id5 = id2; break;
952 case 2: id3 = id2; id4 = id1; id5 = 21; break;
953 case 3: id3 = 21; id4 = id1; id5 = id2; break;
954 case 4: id3 = id2; id4 = 21; id5 = id1; break;
955 case 5: id3 = 21; id4 = id2; id5 = id1; break;
957 setId(id1, id2, id3, id4, id5);
959 // Temporary solution; id1 and id2 can be q/qbar independently
962 cols[0][0] = 1; cols[0][1] = 0;
963 cols[2][0] = 1; cols[2][1] = 0;
965 cols[0][0] = 0; cols[0][1] = 1;
966 cols[2][0] = 0; cols[2][1] = 1;
969 cols[1][0] = 2; cols[1][1] = 0;
970 cols[3][0] = 3; cols[3][1] = 0;
971 cols[4][0] = 2; cols[4][1] = 3;
973 cols[1][0] = 0; cols[1][1] = 2;
974 cols[3][0] = 0; cols[3][1] = 3;
975 cols[4][0] = 3; cols[4][1] = 2;
977 // Map correct final state configuration
978 int i3 = 0, i4 = 0, i5 = 0;
980 case 0: i3 = 2; i4 = 3; i5 = 4; break;
981 case 1: i3 = 2; i4 = 4; i5 = 3; break;
982 case 2: i3 = 3; i4 = 2; i5 = 4; break;
983 case 3: i3 = 4; i4 = 2; i5 = 3; break;
984 case 4: i3 = 3; i4 = 4; i5 = 2; break;
985 case 5: i3 = 4; i4 = 3; i5 = 2; break;
987 // Put colours in place
988 setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
989 cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
990 cols[i5][0], cols[i5][1]);
994 //--------------------------------------------------------------------------
996 // Map a final state configuration
998 inline void Sigma3qq2qqgDiff::mapFinal() {
1000 case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
1001 case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
1002 case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
1003 case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
1004 case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
1005 case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
1010 //==========================================================================
1012 // Sigma3qqbar2qqbargDiff
1013 // Cross section for q qbar -> q' qbar' g
1014 // Crossed relation from q q' -> q q' g, q != q'
1016 //--------------------------------------------------------------------------
1018 // Initialize process.
1020 void Sigma3qqbar2qqbargDiff::initProc() {
1022 // Read number of quarks to be considered in massless approximation.
1023 nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
1027 //--------------------------------------------------------------------------
1029 // Evaluate |M|^2 - no incoming flavour dependence.
1031 void Sigma3qqbar2qqbargDiff::sigmaKin() {
1032 // Overall 6 possibilities for final state ordering
1033 // To keep symmetry between final states, always map to:
1034 // 1) q1(p+) qbar1(p-) -> qbar2(q+) q2(q-) g(k)
1035 // 2) qbar1(p+) q1(p-) -> q2(q+) qbar2(q-) g(k)
1036 // Crossing p- and q+ gives:
1037 // 1) q1(p+) q2(-q+) -> q1(-p-) q2(q-) g(k)
1038 // 2) qbar1(p+) qbar2(-q+) -> qbar1(-p-) qbar2(q-) g(k)
1040 // Incoming four-vectors
1041 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1042 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1043 // Pick and map a final state configuration
1048 swap(pCM[1], pCM[2]);
1053 // Extra factor of (6.) from picking a final state configuration
1054 // Extra factor of (nQuarkNew - 1) from new q/qbar pairs
1055 // XXX - Extra factor of (2.) from second possible crossing???
1056 sigma = 6. * (nQuarkNew - 1) * 2. * m2Calc();
1060 //--------------------------------------------------------------------------
1062 // Select identity, colour and anticolour.
1064 void Sigma3qqbar2qqbargDiff::setIdColAcol(){
1066 // Pick new q qbar flavour with incoming flavour disallowed
1067 int idNew = 1 + int( (nQuarkNew - 1) * rndmPtr->flat() );
1068 if (idNew >= abs(id1)) ++idNew;
1069 // For qbar q incoming, q+ is always mapped to q2
1070 // For q qbar incoming, q+ is always mapped to qbar2
1071 if (id1 > 0) idNew = -idNew;
1073 // Outgoing flavours; easiest just to map by hand
1075 case 0: id3 = idNew; id4 = -idNew; id5 = 21; break;
1076 case 1: id3 = idNew; id4 = 21; id5 = -idNew; break;
1077 case 2: id3 = -idNew; id4 = idNew; id5 = 21; break;
1078 case 3: id3 = 21; id4 = idNew; id5 = -idNew; break;
1079 case 4: id3 = -idNew; id4 = 21; id5 = idNew; break;
1080 case 5: id3 = 21; id4 = -idNew; id5 = idNew; break;
1082 setId(id1, id2, id3, id4, id5);
1084 // Temporary solution; start with q qbar -> qbar q g
1086 cols[0][0] = 1; cols[0][1] = 0;
1087 cols[1][0] = 0; cols[1][1] = 2;
1088 cols[2][0] = 0; cols[2][1] = 3;
1089 cols[3][0] = 1; cols[3][1] = 0;
1090 cols[4][0] = 3; cols[4][1] = 2;
1091 // Map into correct place
1092 int i3 = 0, i4 = 0, i5 = 0;
1094 case 0: i3 = 2; i4 = 3; i5 = 4; break;
1095 case 1: i3 = 2; i4 = 4; i5 = 3; break;
1096 case 2: i3 = 3; i4 = 2; i5 = 4; break;
1097 case 3: i3 = 4; i4 = 2; i5 = 3; break;
1098 case 4: i3 = 3; i4 = 4; i5 = 2; break;
1099 case 5: i3 = 4; i4 = 3; i5 = 2; break;
1101 setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1102 cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1103 cols[i5][0], cols[i5][1]);
1104 // Swap for qbar q incoming
1105 if (id1 < 0) swapColAcol();
1109 //==========================================================================
1111 // Sigma3qg2qqqbarDiff class.
1112 // Cross section for q g -> q q' qbar'
1113 // Crossed relation from q q' -> q q' g, q != q'
1114 // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1116 //--------------------------------------------------------------------------
1118 // Initialize process.
1120 void Sigma3qg2qqqbarDiff::initProc() {
1122 // Read number of quarks to be considered in massless approximation.
1123 nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
1127 //--------------------------------------------------------------------------
1129 // Evaluate |M|^2 - no incoming flavour dependence
1131 void Sigma3qg2qqqbarDiff::sigmaKin() {
1133 // Pick a final state configuration
1136 // gq or qg incoming
1137 for (int i = 0; i < 2; i++) {
1139 // Map incoming four-vectors
1140 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1141 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1144 // Crossing (note extra -ve sign in total sigma)
1145 swap(pCM[i], pCM[4]);
1150 // Extra factor of (6) from picking a final state configuration
1151 // Extra factor of (3 / 8) as averaging over incoming gluon
1152 // Extra factor of (nQuarkNew - 1) due to new q/qbar pair
1153 sigma[i] = -6. * (3. / 8.) * (nQuarkNew - 1) * m2Calc();
1159 //--------------------------------------------------------------------------
1161 // Evaluate |M|^2 - incoming flavour dependence
1163 double Sigma3qg2qqqbarDiff::sigmaHat() {
1164 // gq or qg incoming
1165 return (id1 == 21) ? sigma[0] : sigma[1];
1168 //--------------------------------------------------------------------------
1170 // Select identity, colour and anticolour.
1172 void Sigma3qg2qqqbarDiff::setIdColAcol(){
1173 // Pick new q qbar flavour with incoming flavour disallowed
1174 int sigmaIdx = (id1 == 21) ? 0 : 1;
1175 int idIn = (id1 == 21) ? id2 : id1;
1176 int idNew = 1 + int( (nQuarkNew - 1) * rndmPtr->flat() );
1177 if (idNew >= abs(idIn)) ++idNew;
1179 // qbar instead of q incoming means swap outgoing q/qbar pair
1180 int id3Tmp = idIn, id4Tmp = idNew, id5Tmp = -idNew;
1181 if (idIn < 0) swap(id4Tmp, id5Tmp);
1182 // If g q incoming rather than q g, idIn and idNew
1183 // should be exchanged (see sigmaKin)
1184 if (sigmaIdx == 0) swap(id3Tmp, id4Tmp);
1185 // Outgoing flavours; now just map as if q g incoming
1187 case 0: id3 = id3Tmp; id4 = id4Tmp; id5 = id5Tmp; break;
1188 case 1: id3 = id3Tmp; id4 = id5Tmp; id5 = id4Tmp; break;
1189 case 2: id3 = id4Tmp; id4 = id3Tmp; id5 = id5Tmp; break;
1190 case 3: id3 = id5Tmp; id4 = id3Tmp; id5 = id4Tmp; break;
1191 case 4: id3 = id4Tmp; id4 = id5Tmp; id5 = id3Tmp; break;
1192 case 5: id3 = id5Tmp; id4 = id4Tmp; id5 = id3Tmp; break;
1194 setId(id1, id2, id3, id4, id5);
1196 // Temporary solution; start with either
1197 // g q1 -> q1 q2 qbar2
1198 // g qbar1 -> qbar1 qbar2 q2
1200 cols[0][0] = 1; cols[0][1] = 2;
1202 cols[1][0] = 3; cols[1][1] = 0;
1203 cols[2][0] = 1; cols[2][1] = 0;
1204 cols[3][0] = 3; cols[3][1] = 0;
1205 cols[4][0] = 0; cols[4][1] = 2;
1207 cols[1][0] = 0; cols[1][1] = 3;
1208 cols[2][0] = 0; cols[2][1] = 2;
1209 cols[3][0] = 0; cols[3][1] = 3;
1210 cols[4][0] = 1; cols[4][1] = 0;
1212 // Swap incoming if q/qbar g instead
1214 swap(cols[0][0], cols[1][0]);
1215 swap(cols[0][1], cols[1][1]);
1218 int i3 = 0, i4 = 0, i5 = 0;
1219 if (sigmaIdx == 0) {
1221 case 0: i3 = 3; i4 = 2; i5 = 4; break;
1222 case 1: i3 = 3; i4 = 4; i5 = 2; break;
1223 case 2: i3 = 2; i4 = 3; i5 = 4; break;
1224 case 3: i3 = 4; i4 = 3; i5 = 2; break;
1225 case 4: i3 = 2; i4 = 4; i5 = 3; break;
1226 case 5: i3 = 4; i4 = 2; i5 = 3; break;
1230 case 0: i3 = 2; i4 = 3; i5 = 4; break;
1231 case 1: i3 = 2; i4 = 4; i5 = 3; break;
1232 case 2: i3 = 3; i4 = 2; i5 = 4; break;
1233 case 3: i3 = 4; i4 = 2; i5 = 3; break;
1234 case 4: i3 = 3; i4 = 4; i5 = 2; break;
1235 case 5: i3 = 4; i4 = 3; i5 = 2; break;
1238 setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1239 cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1240 cols[i5][0], cols[i5][1]);
1243 //==========================================================================
1245 // Sigma3qq2qqgSame class.
1246 // Cross section for q q' -> q q' g, q == q'.
1248 //--------------------------------------------------------------------------
1250 // Evaluate |M|^2 - no incoming flavour dependence
1252 void Sigma3qq2qqgSame::sigmaKin() {
1253 // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1255 // Incoming four-vectors
1256 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1257 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1258 // Pick/map a final state configuration
1263 // Extra factor (3) from picking final state configuration
1264 // (original answer already has a factor 2 from identical
1265 // quarks in the final state)
1266 sigma = 3. * m2Calc();
1270 //--------------------------------------------------------------------------
1274 inline double Sigma3qq2qqgSame::m2Calc() {
1277 s = (pCM[0] + pCM[1]).m2Calc();
1278 t = (pCM[0] - pCM[2]).m2Calc();
1279 u = (pCM[0] - pCM[3]).m2Calc();
1280 sp = (pCM[2] + pCM[3]).m2Calc();
1281 tp = (pCM[1] - pCM[3]).m2Calc();
1282 up = (pCM[1] - pCM[2]).m2Calc();
1292 double den1 = (pCM[0] * pCM[4]) * (pCM[1] * pCM[4]) *
1293 (pCM[2] * pCM[4]) * (pCM[3] * pCM[4]);
1295 double fac1 = s * (t * u + tp * up) + sp * (t * up + tp * u);
1296 double fac2 = ssp - ttp - uup;
1297 double fac3 = 2. * (ttp * u_up + uup * t_tp);
1299 double num1 = u_up * (ssp + ttp - uup) + fac1;
1300 double num2 = s_sp * fac2 + fac3;
1301 double num3 = (s * s + sp * sp + u * u + up * up) / (t * tp);
1303 double num4 = t_tp * (ssp - ttp + uup) + fac1;
1304 double num5 = (s * s + sp * sp + t * t + tp * tp) / (u * up);
1306 double num6 = s_sp * fac2 - fac3 - 2. * fac1;
1307 double num7 = (s * s + sp * sp) * fac2;
1308 double den7 = (ttp * uup);
1310 // C1 = (N^2 - 1)^2 / 4N^3 = 16. / 27.
1311 // C2 = (N^2 - 1) / 4N^3 = 2. / 27.
1312 // C3 = (N^4 - 1) / 8N^4 = 10. / 81.
1313 // C4 = (N^2 - 1)^2 / 8N^4 = 8. / 81.
1314 return (1. / 8.) * pow3(4. * M_PI * alpS) *
1315 ( ( (16. / 27.) * num1 - (2. / 27.) * num2 ) * num3 +
1316 ( (16. / 27.) * num4 - (2. / 27.) * num2 ) * num5 +
1317 ( (10. / 81.) * num2 + (8. / 81.) * num6 ) *
1318 ( num7 / den7 ) ) / den1;
1322 //--------------------------------------------------------------------------
1324 // Evaluate |M|^2 - incoming flavour dependence
1326 double Sigma3qq2qqgSame::sigmaHat() {
1327 // q q / qbar qbar incoming states only
1328 if (id1 != id2) return 0.;
1332 //--------------------------------------------------------------------------
1334 // Select identity, colour and anticolour.
1336 void Sigma3qq2qqgSame::setIdColAcol(){
1338 // Need to know where the gluon was mapped (pCM[4])
1341 case 3: case 5: gIdx = 0; break;
1342 case 1: case 4: gIdx = 1; break;
1343 case 0: case 2: gIdx = 2; break;
1346 // Outgoing flavours
1347 int idTmp[3] = { id1, id1, id1 };
1349 setId(id1, id2, idTmp[0], idTmp[1], idTmp[2]);
1351 // Temporary solution; start with q q -> q q g
1352 setColAcol(1, 0, 2, 0, 1, 0, 3, 0, 2, 3);
1354 swap( colSave[5], colSave[gIdx + 3]);
1355 swap(acolSave[5], acolSave[gIdx + 3]);
1356 // Swap if qbar qbar incoming
1357 if (id1 < 0) swapColAcol();
1361 //--------------------------------------------------------------------------
1363 // Map a final state configuration
1364 inline void Sigma3qq2qqgSame::mapFinal() {
1366 case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
1367 case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
1368 case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
1369 case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
1370 case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
1371 case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
1375 //==========================================================================
1377 // Sigma3qqbar2qqbargSame class.
1378 // Cross section for q qbar -> q qbar g
1379 // Crossed relation from q(bar) q(bar) -> q(bar) q(bar) g:
1380 // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1381 //--------------------------------------------------------------------------
1383 // Evaluate |M|^2 - no incoming flavour dependence
1385 void Sigma3qqbar2qqbargSame::sigmaKin() {
1387 // Incoming four-vectors
1388 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1389 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1391 // Pick and map a final state configuration
1396 swap(pCM[1], pCM[3]);
1401 // Extra factor of (6) from picking a final state configuration
1402 sigma = 6. * m2Calc();
1406 //--------------------------------------------------------------------------
1408 // Select identity, colour and anticolour.
1410 void Sigma3qqbar2qqbargSame::setIdColAcol(){
1411 // Outgoing flavours; easiest to map by hand
1413 case 0: id3 = id1; id4 = id2; id5 = 21; break;
1414 case 1: id3 = id1; id4 = 21; id5 = id2; break;
1415 case 2: id3 = id2; id4 = id1; id5 = 21; break;
1416 case 3: id3 = 21; id4 = id1; id5 = id2; break;
1417 case 4: id3 = id2; id4 = 21; id5 = id1; break;
1418 case 5: id3 = 21; id4 = id2; id5 = id1; break;
1420 setId(id1, id2, id3, id4, id5);
1422 // Temporary solution; start with q qbar -> q qbar g
1424 cols[0][0] = 1; cols[0][1] = 0;
1425 cols[1][0] = 0; cols[1][1] = 2;
1426 cols[2][0] = 1; cols[2][1] = 0;
1427 cols[3][0] = 0; cols[3][1] = 3;
1428 cols[4][0] = 3; cols[4][1] = 2;
1430 int i3 = 0, i4 = 0, i5 = 0;
1432 case 0: i3 = 2; i4 = 3; i5 = 4; break;
1433 case 1: i3 = 2; i4 = 4; i5 = 3; break;
1434 case 2: i3 = 3; i4 = 2; i5 = 4; break;
1435 case 3: i3 = 4; i4 = 2; i5 = 3; break;
1436 case 4: i3 = 3; i4 = 4; i5 = 2; break;
1437 case 5: i3 = 4; i4 = 3; i5 = 2; break;
1439 setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1440 cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1441 cols[i5][0], cols[i5][1]);
1442 // Swap for qbar q incoming
1443 if (id1 < 0) swapColAcol();
1446 //==========================================================================
1448 // Sigma3qg2qqqbarSame class.
1449 // Cross section for q g -> q q qbar.
1450 // Crossed relation from q(bar) q(bar) -> q(bar) q(bar) g:
1451 // q1(p+) q1(p-) -> q1(q+) q1(q-) g(k)
1453 //--------------------------------------------------------------------------
1455 // Evaluate |M|^2 - no incoming flavour dependence
1457 void Sigma3qg2qqqbarSame::sigmaKin() {
1459 // Pick a final state configuration
1462 // gq and qg incoming
1463 for (int i = 0; i < 2; i++) {
1465 // Map incoming four-vectors
1466 pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1467 pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1470 // Crossing (note extra -ve sign in total sigma)
1471 swap(pCM[i], pCM[4]);
1476 // XXX - Extra factor of (3) from picking a final state configuration???
1477 // Extra factor of (3 / 8) as averaging over incoming gluon
1478 sigma[i] = -3. * (3. / 8.) * m2Calc();
1484 //--------------------------------------------------------------------------
1486 // Evaluate |M|^2 - incoming flavour dependence
1488 double Sigma3qg2qqqbarSame::sigmaHat() {
1489 // gq or qg incoming
1490 return (id1 == 21) ? sigma[0] : sigma[1];
1493 //--------------------------------------------------------------------------
1495 // Select identity, colour and anticolour.
1497 void Sigma3qg2qqqbarSame::setIdColAcol(){
1499 // Pick outgoing flavour configuration
1500 int idIn = (id1 == 21) ? id2 : id1;
1502 // Outgoing flavours; easiest just to map by hand
1504 case 0: id3 = idIn; id4 = idIn; id5 = -idIn; break;
1505 case 1: id3 = idIn; id4 = -idIn; id5 = idIn; break;
1506 case 2: id3 = idIn; id4 = idIn; id5 = -idIn; break;
1507 case 3: id3 = -idIn; id4 = idIn; id5 = idIn; break;
1508 case 4: id3 = idIn; id4 = -idIn; id5 = idIn; break;
1509 case 5: id3 = -idIn; id4 = idIn; id5 = idIn; break;
1511 setId(id1, id2, id3, id4, id5);
1513 // Temporary solution; start with either
1514 // g q1 -> q1 q2 qbar2
1515 // g qbar1 -> qbar1 qbar2 q2
1517 cols[0][0] = 1; cols[0][1] = 2;
1519 cols[1][0] = 3; cols[1][1] = 0;
1520 cols[2][0] = 1; cols[2][1] = 0;
1521 cols[3][0] = 3; cols[3][1] = 0;
1522 cols[4][0] = 0; cols[4][1] = 2;
1524 cols[1][0] = 0; cols[1][1] = 3;
1525 cols[2][0] = 0; cols[2][1] = 2;
1526 cols[3][0] = 0; cols[3][1] = 3;
1527 cols[4][0] = 1; cols[4][1] = 0;
1529 // Swap incoming if q/qbar g instead
1531 swap(cols[0][0], cols[1][0]);
1532 swap(cols[0][1], cols[1][1]);
1535 int i3 = 0, i4 = 0, i5 = 0;
1537 case 0: i3 = 2; i4 = 3; i5 = 4; break;
1538 case 1: i3 = 2; i4 = 4; i5 = 3; break;
1539 case 2: i3 = 3; i4 = 2; i5 = 4; break;
1540 case 3: i3 = 4; i4 = 2; i5 = 3; break;
1541 case 4: i3 = 3; i4 = 4; i5 = 2; break;
1542 case 5: i3 = 4; i4 = 3; i5 = 2; break;
1544 setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1545 cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1546 cols[i5][0], cols[i5][1]);
1549 //==========================================================================
1551 } // end namespace Pythia8