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) 2014 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 "Pythia8/SigmaQCD.h"
10 
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // Sigma0AB2AB class.
16 // Cross section for elastic scattering A B -> A B.
17 
18 //--------------------------------------------------------------------------
19 
20 // Select identity, colour and anticolour.
21 
22 void Sigma0AB2AB::setIdColAcol() {
23 
24  // Flavours and colours are trivial.
25  setId( idA, idB, idA, idB);
26  setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
27 }
28 
29 //==========================================================================
30 
31 // Sigma0AB2XB class.
32 // Cross section for single diffractive scattering A B -> X B.
33 
34 //--------------------------------------------------------------------------
35 
36 // Select identity, colour and anticolour.
37 
38 void Sigma0AB2XB::setIdColAcol() {
39 
40  // Flavours and colours are trivial.
41  int idX = 10* (abs(idA) / 10) + 9900000;
42  if (idA < 0) idX = -idX;
43  setId( idA, idB, idX, idB);
44  setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
45 
46 }
47 
48 //==========================================================================
49 
50 // Sigma0AB2AX class.
51 // Cross section for single diffractive scattering A B -> A X.
52 
53 //--------------------------------------------------------------------------
54 
55 // Select identity, colour and anticolour.
56 
57 void Sigma0AB2AX::setIdColAcol() {
58 
59  // Flavours and colours are trivial.
60  int idX = 10* (abs(idB) / 10) + 9900000;
61  if (idB < 0) idX = -idX;
62  setId( idA, idB, idA, idX);
63  setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
64 
65 }
66 
67 //==========================================================================
68 
69 // Sigma0AB2XX class.
70 // Cross section for double diffractive scattering A B -> X X.
71 
72 //--------------------------------------------------------------------------
73 
74 // Select identity, colour and anticolour.
75 
76 void Sigma0AB2XX::setIdColAcol() {
77 
78  // Flavours and colours are trivial.
79  int idX1 = 10* (abs(idA) / 10) + 9900000;
80  if (idA < 0) idX1 = -idX1;
81  int idX2 = 10* (abs(idB) / 10) + 9900000;
82  if (idB < 0) idX2 = -idX2;
83  setId( idA, idB, idX1, idX2);
84  setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
85 
86 }
87 
88 //==========================================================================
89 
90 // Sigma0AB2AXB class.
91 // Cross section for central scattering A B -> A X B.
92 
93 //--------------------------------------------------------------------------
94 
95 // Select identity, colour and anticolour.
96 
97 void Sigma0AB2AXB::setIdColAcol() {
98 
99  // Central diffractive state represented by rho_diffr0. Colours trivial.
100  int idX = 9900110;
101  setId( idA, idB, idA, idB,idX);
102  setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
103 
104 }
105 
106 //==========================================================================
107 
108 // Sigma2gg2gg class.
109 // Cross section for g g -> g g.
110 
111 //--------------------------------------------------------------------------
112 
113 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
114 
115 void Sigma2gg2gg::sigmaKin() {
116 
117  // Calculate kinematics dependence.
118  sigTS = (9./4.) * (tH2 / sH2 + 2. * tH / sH + 3. + 2. * sH / tH
119  + sH2 / tH2);
120  sigUS = (9./4.) * (uH2 / sH2 + 2. * uH / sH + 3. + 2. * sH / uH
121  + sH2 / uH2);
122  sigTU = (9./4.) * (tH2 / uH2 + 2. * tH / uH + 3. + 2. * uH / tH
123  + uH2 / tH2);
124  sigSum = sigTS + sigUS + sigTU;
125 
126  // Answer contains factor 1/2 from identical gluons.
127  sigma = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
128 
129 }
130 
131 //--------------------------------------------------------------------------
132 
133 // Select identity, colour and anticolour.
134 
135 void Sigma2gg2gg::setIdColAcol() {
136 
137  // Flavours are trivial.
138  setId( id1, id2, 21, 21);
139 
140  // Three colour flow topologies, each with two orientations.
141  double sigRand = sigSum * rndmPtr->flat();
142  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
143  else if (sigRand < sigTS + sigUS)
144  setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
145  else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
146  if (rndmPtr->flat() > 0.5) swapColAcol();
147 
148 }
149 
150 //==========================================================================
151 
152 // Sigma2gg2qqbar class.
153 // Cross section for g g -> q qbar (q = u, d, s, i.e. almost massless).
154 
155 //--------------------------------------------------------------------------
156 
157 // Initialize process.
158 
159 void Sigma2gg2qqbar::initProc() {
160 
161  // Read number of quarks to be considered in massless approximation.
162  nQuarkNew = 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
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5