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