StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaOnia.cc
1 // SigmaOnia.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 // charmonia/bottomonia simulation classes.
8 
9 #include "SigmaOnia.h"
10 
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // Sigma2gg2QQbar3S11g class.
16 // Cross section g g -> QQbar[3S1(1)] g (Q = c or b).
17 
18 //--------------------------------------------------------------------------
19 
20 // Initialize process.
21 
22 void Sigma2gg2QQbar3S11g::initProc() {
23 
24  // Produced state. Process name. Onium matrix element.
25  idHad = (idNew == 4) ? 443 : 553;
26  nameSave = (idNew == 4) ? "g g -> ccbar[3S1(1)] g"
27  : "g g -> bbbar[3S1(1)] g";
28  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3S11")
29  : settingsPtr->parm("Bottomonium:OUpsilon3S11");
30 
31 }
32 
33 //--------------------------------------------------------------------------
34 
35 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
36 
37 void Sigma2gg2QQbar3S11g::sigmaKin() {
38 
39  // Calculate kinematics dependence.
40  double stH = sH + tH;
41  double tuH = tH + uH;
42  double usH = uH + sH;
43  double sig = (10. * M_PI / 81.) * m3 * ( pow2(sH * tuH)
44  + pow2(tH * usH) + pow2(uH * stH) ) / pow2( stH * tuH * usH );
45 
46  // Answer.
47  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
48 
49 }
50 
51 //--------------------------------------------------------------------------
52 
53 // Select identity, colour and anticolour.
54 
55 void Sigma2gg2QQbar3S11g::setIdColAcol() {
56 
57  // Flavours are trivial.
58  setId( id1, id2, idHad, 21);
59 
60  // Two orientations of colour flow..
61  setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
62  if (rndmPtr->flat() > 0.5) swapColAcol();
63 
64 }
65 
66 //==========================================================================
67 
68 // Sigma2gg2QQbar3PJ1g class.
69 // Cross section g g -> QQbar[3PJ(1)] g (Q = c or b, J = 0, 1 or 2).
70 
71 //--------------------------------------------------------------------------
72 
73 // Initialize process.
74 
75 void Sigma2gg2QQbar3PJ1g::initProc() {
76 
77  // Produced state. Process name. Onium matrix element.
78  idHad = 0;
79  nameSave = "illegal process";
80  if (jSave == 0) {
81  idHad = (idNew == 4) ? 10441 : 10551;
82  nameSave = (idNew == 4) ? "g g -> ccbar[3P0(1)] g"
83  : "g g -> bbbar[3P0(1)] g";
84  } else if (jSave == 1) {
85  idHad = (idNew == 4) ? 20443 : 20553;
86  nameSave = (idNew == 4) ? "g g -> ccbar[3P1(1)] g"
87  : "g g -> bbbar[3P1(1)] g";
88  } else if (jSave == 2) {
89  idHad = (idNew == 4) ? 445 : 555;
90  nameSave = (idNew == 4) ? "g g -> ccbar[3P2(1)] g"
91  : "g g -> bbbar[3P2(1)] g";
92  }
93  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:Ochic03P01")
94  : settingsPtr->parm("Bottomonium:Ochib03P01");
95 
96 }
97 
98 //--------------------------------------------------------------------------
99 
100 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
101 
102 void Sigma2gg2QQbar3PJ1g::sigmaKin() {
103 
104  // Useful derived kinematics quantities.
105  double pRat = (sH * uH + uH * tH + tH * sH)/ sH2;
106  double qRat = tH * uH / sH2;
107  double rRat = s3 / sH;
108  double pRat2 = pRat * pRat;
109  double pRat3 = pRat2 * pRat;
110  double pRat4 = pRat3 * pRat;
111  double qRat2 = qRat * qRat;
112  double qRat3 = qRat2 * qRat;
113  double qRat4 = qRat3 * qRat;
114  double rRat2 = rRat * rRat;
115  double rRat3 = rRat2 * rRat;
116  double rRat4 = rRat3 * rRat;
117 
118  // Calculate kinematics dependence.
119  double sig = 0.;
120  if (jSave == 0) {
121  sig = (8. * M_PI / (9. * m3 * sH))
122  * ( 9. * rRat2 * pRat4 * (rRat4 - 2. * rRat2 * pRat + pRat2)
123  - 6. * rRat * pRat3 * qRat * (2. * rRat4 - 5. * rRat2 * pRat
124  + pRat2) - pRat2 * qRat2 * (rRat4 + 2. * rRat2 * pRat - pRat2)
125  + 2. * rRat * pRat * qRat3 * (rRat2 - pRat) + 6. * rRat2 * qRat4)
126  / (qRat * pow4(qRat - rRat * pRat));
127  } else if (jSave == 1) {
128  sig = (8. * M_PI / (3.* m3 * sH)) * pRat2
129  * (rRat * pRat2 * (rRat2 - 4. * pRat)
130  + 2. * qRat * (-rRat4 + 5. * rRat2 * pRat + pRat2)
131  - 15. * rRat * qRat2) / pow4(qRat - rRat * pRat);
132  } else if (jSave == 2) {
133  sig = (8. * M_PI / (9. * m3 * sH))
134  * (12. * rRat2 * pRat4 * (rRat4 - 2. * rRat2 * pRat + pRat2)
135  - 3. * rRat * pRat3 * qRat * (8. * rRat4 - rRat2 * pRat + 4. * pRat2)
136  + 2. * pRat2 * qRat2 * (-7. * rRat4 + 43. * rRat2 * pRat + pRat2)
137  + rRat * pRat * qRat3 * (16. * rRat2 - 61. * pRat)
138  + 12. * rRat2 * qRat4) / (qRat * pow4(qRat-rRat * pRat));
139  }
140 
141  // Answer.
142  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
143 
144 }
145 
146 //--------------------------------------------------------------------------
147 
148 // Select identity, colour and anticolour.
149 
150 void Sigma2gg2QQbar3PJ1g::setIdColAcol() {
151 
152  // Flavours are trivial.
153  setId( id1, id2, idHad, 21);
154 
155  // Two orientations of colour flow.
156  setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
157  if (rndmPtr->flat() > 0.5) swapColAcol();
158 
159 }
160 
161 //==========================================================================
162 
163 // Sigma2qg2QQbar3PJ1q class.
164 // Cross section q g -> QQbar[3PJ(1)] q (Q = c or b, J = 0, 1 or 2).
165 
166 //--------------------------------------------------------------------------
167 
168 // Initialize process.
169 
170 void Sigma2qg2QQbar3PJ1q::initProc() {
171 
172  // Produced state. Process name. Onium matrix element.
173  idHad = 0;
174  nameSave = "illegal process";
175  if (jSave == 0) {
176  idHad = (idNew == 4) ? 10441 : 10551;
177  nameSave = (idNew == 4) ? "q g -> ccbar[3P0(1)] q"
178  : "q g -> bbbar[3P0(1)] q";
179  } else if (jSave == 1) {
180  idHad = (idNew == 4) ? 20443 : 20553;
181  nameSave = (idNew == 4) ? "q g -> ccbar[3P1(1)] q"
182  : "q g -> bbbar[3P1(1)] q";
183  } else if (jSave == 2) {
184  idHad = (idNew == 4) ? 445 : 555;
185  nameSave = (idNew == 4) ? "q g -> ccbar[3P2(1)] q"
186  : "q g -> bbbar[3P2(1)] q";
187  }
188  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:Ochic03P01")
189  : settingsPtr->parm("Bottomonium:Ochib03P01");
190 
191 }
192 
193 //--------------------------------------------------------------------------
194 
195 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
196 
197 void Sigma2qg2QQbar3PJ1q::sigmaKin() {
198 
199  // Calculate kinematics dependence.
200  double usH = uH + sH;
201  double sig = 0.;
202  if (jSave == 0) {
203  sig = - (16. * M_PI / 81.) * pow2(tH - 3. * s3) * (sH2 + uH2)
204  / (m3 * tH * pow4(usH));
205  } else if (jSave == 1) {
206  sig = - (32. * M_PI / 27.) * (4. * s3 * sH * uH + tH * (sH2 + uH2))
207  / (m3 * pow4(usH));
208  } else if (jSave == 2) {
209  sig = - (32. *M_PI / 81.) * ( (6. * s3*s3 + tH2) * pow2(usH)
210  - 2. * sH * uH * (tH2 + 6. * s3 * usH)) / (m3 * tH * pow4(usH));
211  }
212 
213  // Answer.
214  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
215 
216 }
217 
218 //--------------------------------------------------------------------------
219 
220 // Select identity, colour and anticolour.
221 
222 void Sigma2qg2QQbar3PJ1q::setIdColAcol() {
223 
224  // Flavours are trivial.
225  int idq = (id2 == 21) ? id1 : id2;
226  setId( id1, id2, idHad, idq);
227 
228  // tH defined between q_in and q_out: must swap tHat <-> uHat if q g in.
229  swapTU = (id2 == 21);
230 
231  // Colour flow topologies. Swap when antiquarks.
232  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
233  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
234  if (idq < 0) swapColAcol();
235 
236 }
237 
238 //==========================================================================
239 
240 // Sigma2qqbar2QQbar3PJ1g class.
241 // Cross section q qbar -> QQbar[3PJ(1)] g (Q = c or b, J = 0, 1 or 2).
242 
243 //--------------------------------------------------------------------------
244 
245 // Initialize process.
246 
247 void Sigma2qqbar2QQbar3PJ1g::initProc() {
248 
249  // Produced state. Process name. Onium matrix element.
250  idHad = 0;
251  nameSave = "illegal process";
252  if (jSave == 0) {
253  idHad = (idNew == 4) ? 10441 : 10551;
254  nameSave = (idNew == 4) ? "q qbar -> ccbar[3P0(1)] g"
255  : "q qbar -> bbbar[3P0(1)] g";
256  } else if (jSave == 1) {
257  idHad = (idNew == 4) ? 20443 : 20553;
258  nameSave = (idNew == 4) ? "q qbar -> ccbar[3P1(1)] g"
259  : "q qbar -> bbbar[3P1(1)] g";
260  } else if (jSave == 2) {
261  idHad = (idNew == 4) ? 445 : 555;
262  nameSave = (idNew == 4) ? "q qbar -> ccbar[3P2(1)] g"
263  : "q qbar -> bbbar[3P2(1)] g";
264  }
265  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:Ochic03P01")
266  : settingsPtr->parm("Bottomonium:Ochib03P01");
267 
268 }
269 
270 //--------------------------------------------------------------------------
271 
272 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
273 
274 void Sigma2qqbar2QQbar3PJ1g::sigmaKin() {
275 
276  // Calculate kinematics dependence.
277  double tuH = tH + uH;
278  double sig = 0.;
279  if (jSave == 0) {
280  sig =(128. * M_PI / 243.) * pow2(sH - 3. * s3) * (tH2 + uH2)
281  / (m3 * sH * pow4(tuH));
282  } else if (jSave == 1) {
283  sig = (256. * M_PI / 81.) * (4. * s3 * tH * uH + sH * (tH2 + uH2))
284  / (m3 * pow4(tuH));
285  } else if (jSave == 2) {
286  sig = (256. * M_PI / 243.) * ( (6. * s3*s3 + sH2) * pow2(tuH)
287  - 2. * tH * uH * (sH2 + 6. * s3 * tuH) )/ (m3 * sH * pow4(tuH));
288  }
289 
290  // Answer.
291  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
292 
293 }
294 
295 //--------------------------------------------------------------------------
296 
297 // Select identity, colour and anticolour.
298 
299 void Sigma2qqbar2QQbar3PJ1g::setIdColAcol() {
300 
301  // Flavours are trivial.
302  setId( id1, id2, idHad, 21);
303 
304  // Colour flow topologies. Swap when antiquarks.
305  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
306  if (id1 < 0) swapColAcol();
307 
308 }
309 
310 //==========================================================================
311 
312 // Sigma2gg2QQbarX8g class.
313 // Cross section g g -> QQbar[X(8)] g (Q = c or b, X = 3S1, 1S0 or 3PJ).
314 
315 //--------------------------------------------------------------------------
316 
317 // Initialize process.
318 
319 void Sigma2gg2QQbarX8g::initProc() {
320 
321  // Produced state. Process name. Onium matrix element.
322  idHad = 0;
323  nameSave = "illegal process";
324  if (stateSave == 0) {
325  idHad = (idNew == 4) ? 9900443 : 9900553;
326  nameSave = (idNew == 4) ? "g g -> ccbar[3S1(8)] g"
327  : "g g -> bbbar[3S1(8)] g";
328  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3S18")
329  : settingsPtr->parm("Bottomonium:OUpsilon3S18");
330  } else if (stateSave == 1) {
331  idHad = (idNew == 4) ? 9900441 : 9900551;
332  nameSave = (idNew == 4) ? "g g -> ccbar[1S0(8)] g"
333  : "g g -> bbbar[1S0(8)] g";
334  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi1S08")
335  : settingsPtr->parm("Bottomonium:OUpsilon1S08");
336  } else if (stateSave == 2) {
337  idHad = (idNew == 4) ? 9910441 : 9910551;
338  nameSave = (idNew == 4) ? "g g -> ccbar[3PJ(8)] g"
339  : "g g -> bbbar[3PJ(8)] g";
340  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3P08")
341  : settingsPtr->parm("Bottomonium:OUpsilon3P08");
342  }
343 
344 }
345 
346 //--------------------------------------------------------------------------
347 
348 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
349 
350 void Sigma2gg2QQbarX8g::sigmaKin() {
351 
352  // Calculate kinematics dependence.
353  double stH = sH + tH;
354  double tuH = tH + uH;
355  double usH = uH + sH;
356  double sig = 0.;
357  if (stateSave == 0) {
358  sig = (M_PI / 72.) * m3 * ( 27. * (pow2(stH) + pow2(tuH)
359  + pow2(usH)) / (s3*s3) - 16. ) * ( pow2(sH * tuH)
360  + pow2(tH * usH) + pow2(uH * stH) ) / pow2( stH * tuH * usH );
361  } else if (stateSave == 1) {
362  sig = (5. * M_PI / 16.) * m3 * ( pow2(uH / (tuH * usH))
363  + pow2(sH / (stH * usH)) + pow2(tH / (stH * tuH)) ) * ( 12.
364  + (pow4(stH) + pow4(tuH) + pow4(usH)) / (s3 * sH * tH * uH) );
365  } else if (stateSave == 2) {
366  double sH3 = sH2 * sH;
367  double sH4 = sH3 * sH;
368  double sH5 = sH4 * sH;
369  double sH6 = sH5 * sH;
370  double sH7 = sH6 * sH;
371  double sH8 = sH7 * sH;
372  double tH3 = tH2 * tH;
373  double tH4 = tH3 * tH;
374  double tH5 = tH4 * tH;
375  double tH6 = tH5 * tH;
376  double tH7 = tH6 * tH;
377  double tH8 = tH7 * tH;
378  double ssttH = sH * sH + sH * tH + tH * tH;
379  sig = 5. * M_PI * (3. * sH * tH * stH * pow4(ssttH)
380  - s3 * pow2(ssttH) * (7. * sH6 + 36. * sH5 * tH + 45. * sH4 * tH2
381  + 28. * sH3 * tH3 + 45. * sH2 * tH4 + 36. * sH * tH5 + 7. * tH6)
382  + pow2(s3) * stH * (35. *sH8 + 169. * sH7 * tH + 299. * sH6 * tH2
383  + 401. * sH5 * tH3 + 418. * sH4 * tH4 + 401. * sH3 * tH5
384  + 299. * sH2 * tH6 + 169. * sH * tH7 + 35. * tH8)
385  - pow3(s3) * (84. *sH8+432. *sH7*tH+905. *sH6*tH2
386  + 1287. * sH5 * tH3 + 1436. * sH4 * tH4 +1287. * sH3 * tH5
387  + 905. * sH2 * tH6 + 432. * sH * tH7 + 84. * tH8)
388  + pow4(s3) * stH * (126. * sH6 + 451. * sH5 * tH +677. * sH4 * tH2
389  + 836. * sH3 * tH3 + 677. * sH2 * tH4 + 451. * sH * tH5
390  + 126. * tH6)
391  - pow5(s3) * 3. * (42. * sH6 + 171. * sH5 * tH + 304. * sH4 * tH2
392  + 362. * sH3 * tH3 + 304. * sH2 * tH4 + 171. * sH * tH5
393  + 42. * tH6)
394  + pow3(s3 * s3) * 2. * stH * (42. * sH4 + 106. * sH3 * tH
395  + 119. * sH2 * tH2 + 106. * sH * tH3 + 42. * tH4)
396  - pow4(s3) * pow3(s3) * (35. * sH4 + 99. * sH3 * tH
397  + 120. * sH2 * tH2 + 99. * sH * tH3 + 35. * tH4)
398  + pow4(s3 * s3) * 7. * stH * ssttH)
399  / (sH * tH * uH * s3 * m3 * pow3(stH * tuH * usH));
400  }
401 
402  // Answer.
403  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
404 
405 }
406 
407 //--------------------------------------------------------------------------
408 
409 // Select identity, colour and anticolour.
410 
411 void Sigma2gg2QQbarX8g::setIdColAcol() {
412 
413  // Flavours are trivial.
414  setId( id1, id2, idHad, 21);
415 
416  // Split total contribution into different colour flows just like in
417  // g g -> g g (with kinematics recalculated for massless partons).
418  double sHr = - (tH + uH);
419  double sH2r = sHr * sHr;
420  double sigTS = tH2/sH2r + 2.*tH/sHr + 3. + 2.*sHr/tH + sH2r/tH2;
421  double sigUS = uH2/sH2r + 2.*uH/sHr + 3. + 2.*sHr/uH + sH2r/uH2;
422  double sigTU = tH2/uH2 + 2.*tH/uH + 3. + 2.*uH/tH + uH2/tH2;
423  double sigSum = sigTS + sigUS + sigTU;
424 
425  // Three colour flow topologies, each with two orientations.
426  double sigRand = sigSum * rndmPtr->flat();
427  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
428  else if (sigRand < sigTS + sigUS)
429  setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
430  else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
431  if (rndmPtr->flat() > 0.5) swapColAcol();
432 
433 
434 }
435 
436 //==========================================================================
437 
438 // Sigma2qg2QQbarX8q class.
439 // Cross section q g -> QQbar[X(8)] q (Q = c or b, X = 3S1, 1S0 or 3PJ).
440 
441 //--------------------------------------------------------------------------
442 
443 // Initialize process.
444 
445 void Sigma2qg2QQbarX8q::initProc() {
446 
447  // Produced state. Process name. Onium matrix element.
448  idHad = 0;
449  nameSave = "illegal process";
450  if (stateSave == 0) {
451  idHad = (idNew == 4) ? 9900443 : 9900553;
452  nameSave = (idNew == 4) ? "q g -> ccbar[3S1(8)] q"
453  : "q g -> bbbar[3S1(8)] q";
454  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3S18")
455  : settingsPtr->parm("Bottomonium:OUpsilon3S18");
456  } else if (stateSave == 1) {
457  idHad = (idNew == 4) ? 9900441 : 9900551;
458  nameSave = (idNew == 4) ? "q g -> ccbar[1S0(8)] q"
459  : "q g -> bbbar[1S0(8)] q";
460  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi1S08")
461  : settingsPtr->parm("Bottomonium:OUpsilon1S08");
462  } else if (stateSave == 2) {
463  idHad = (idNew == 4) ? 9910441 : 9910551;
464  nameSave = (idNew == 4) ? "q g -> ccbar[3PJ(8)] q"
465  : "q g -> bbbar[3PJ(8)] q";
466  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3P08")
467  : settingsPtr->parm("Bottomonium:OUpsilon3P08");
468  }
469 
470 }
471 
472 //--------------------------------------------------------------------------
473 
474 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
475 
476 void Sigma2qg2QQbarX8q::sigmaKin() {
477 
478  // Calculate kinematics dependence.
479  double stH = sH + tH;
480  double tuH = tH + uH;
481  double usH = uH + sH;
482  double stH2 = stH * stH;
483  double tuH2 = tuH * tuH;
484  double usH2 = usH * usH;
485  double sig = 0.;
486  if (stateSave == 0) {
487  sig = - (M_PI / 27.)* (4. * (sH2 + uH2) - sH * uH) * (stH2 +tuH2)
488  / (s3 * m3 * sH * uH * usH2);
489  } else if (stateSave == 1) {
490  sig = - (5. * M_PI / 18.) * (sH2 + uH2) / (m3 * tH * usH2);
491  } else if (stateSave == 2) {
492  sig = - (10. * M_PI / 9.) * ( (7. * usH + 8. * tH) * (sH2 + uH2)
493  + 4. * tH * (2. * pow2(s3) - stH2 - tuH2) )
494  / (s3 * m3 * tH * usH2 * usH);
495  }
496 
497  // Answer.
498  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
499 
500 }
501 
502 //--------------------------------------------------------------------------
503 
504 // Select identity, colour and anticolour.
505 
506 void Sigma2qg2QQbarX8q::setIdColAcol() {
507 
508  // Flavours are trivial.
509  int idq = (id2 == 21) ? id1 : id2;
510  setId( id1, id2, idHad, idq);
511 
512  // tH defined between q_in and q_out: must swap tHat <-> uHat if q g in.
513  swapTU = (id2 == 21);
514 
515  // Split total contribution into different colour flows just like in
516  // q g -> q g (with kinematics recalculated for massless partons).
517  double sHr = - (tH + uH);
518  double sH2r = sHr * sHr;
519  double sigTS = uH2/tH2 - (4./9.) * uH/sHr;
520  double sigTU = sH2r/tH2 - (4./9.) * sHr/uH;
521  double sigSum = sigTS + sigTU;
522 
523  // Two colour flow topologies. Swap if first is gluon, or when antiquark.
524  double sigRand = sigSum * rndmPtr->flat();
525  if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 2, 3, 3, 0);
526  else setColAcol( 1, 0, 2, 3, 1, 3, 2, 0);
527  if (id1 == 21) swapCol12();
528  if (idq < 0) swapColAcol();
529 
530 }
531 
532 //==========================================================================
533 
534 // Sigma2qqbar2QQbarX8g class.
535 // Cross section q qbar -> QQbar[X(8)] g (Q = c or b, X = 3S1, 1S0 or 3PJ).
536 
537 //--------------------------------------------------------------------------
538 
539 // Initialize process.
540 
541 void Sigma2qqbar2QQbarX8g::initProc() {
542 
543  // Produced state. Process name. Onium matrix element.
544  idHad = 0;
545  nameSave = "illegal process";
546  if (stateSave == 0) {
547  idHad = (idNew == 4) ? 9900443 : 9900553;
548  nameSave = (idNew == 4) ? "q qbar -> ccbar[3S1(8)] g"
549  : "q qbar -> bbbar[3S1(8)] g";
550  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3S18")
551  : settingsPtr->parm("Bottomonium:OUpsilon3S18");
552  } else if (stateSave == 1) {
553  idHad = (idNew == 4) ? 9900441 : 9900551;
554  nameSave = (idNew == 4) ? "q qbar -> ccbar[1S0(8)] g"
555  : "q qbar -> bbbar[1S0(8)] g";
556  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi1S08")
557  : settingsPtr->parm("Bottomonium:OUpsilon1S08");
558  } else if (stateSave == 2) {
559  idHad = (idNew == 4) ? 9910441 : 9910551;
560  nameSave = (idNew == 4) ? "q qbar -> ccbar[3PJ(8)] g"
561  : "q qbar -> bbbar[3PJ(8)] g";
562  oniumME = (idNew == 4) ? settingsPtr->parm("Charmonium:OJpsi3P08")
563  : settingsPtr->parm("Bottomonium:OUpsilon3P08");
564  }
565 
566 }
567 
568 //--------------------------------------------------------------------------
569 
570 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
571 
572 void Sigma2qqbar2QQbarX8g::sigmaKin() {
573 
574  // Calculate kinematics dependence.
575  double stH = sH + tH;
576  double tuH = tH + uH;
577  double usH = uH + sH;
578  double stH2 = stH * stH;
579  double tuH2 = tuH * tuH;
580  double usH2 = usH * usH;
581  double sig = 0.;
582  if (stateSave == 0) {
583  sig = (8. * M_PI / 81.) * (4. * (tH2 + uH2) - tH * uH)
584  * (stH2 + usH2) / (s3 * m3 * tH * uH * tuH2);
585  } else if (stateSave == 1) {
586  sig = (20. * M_PI / 27.) * (tH2 + uH2) / (m3 * sH * tuH2);
587  } else if (stateSave == 2) {
588  sig = (80. * M_PI / 27.) * ( (7. * tuH + 8. * sH) * (tH2 + uH2)
589  + 4. * sH * (2. * pow2(s3) - stH2 -usH2) )
590  / (s3 * m3 * sH * tuH2 * tuH);
591  }
592 
593  // Answer.
594  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
595 
596 }
597 
598 //--------------------------------------------------------------------------
599 
600 // Select identity, colour and anticolour.
601 
602 void Sigma2qqbar2QQbarX8g::setIdColAcol() {
603 
604  // Flavours are trivial.
605  setId( id1, id2, idHad, 21);
606 
607  // Split total contribution into different colour flows just like in
608  // q qbar -> g g (with kinematics recalculated for massless partons).
609  double sHr = - (tH + uH);
610  double sH2r = sHr * sHr;
611  double sigTS = (4. / 9.) * uH / tH - uH2 / sH2r;
612  double sigUS = (4. / 9.) * tH / uH - tH2 / sH2r;
613  double sigSum = sigTS + sigUS;
614 
615  // Two colour flow topologies. Swap if first is antiquark.
616  double sigRand = sigSum * rndmPtr->flat();
617  if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
618  else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
619  if (id1 < 0) swapColAcol();
620 
621 }
622 
623 //==========================================================================
624 
625 } // end namespace Pythia8
626