StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaQCD.cc
1 // SigmaQCD.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 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 // 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 = 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. * 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 = 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
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5