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