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