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) 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 // charmonia/bottomonia simulation classes.
8 
9 #include "Pythia8/SigmaOnia.h"
10 #include <limits>
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // SigmaOniaSetup class.
16 // A helper class used to setup the SigmaOnia processes.
17 
18 //--------------------------------------------------------------------------
19 
20 // The constructor.
21 
22 SigmaOniaSetup::SigmaOniaSetup(Info* infoPtrIn, Settings* settingsPtrIn,
23  ParticleData* particleDataPtrIn, int flavourIn)
24  : valid3S1(true), valid3PJ(true), valid3DJ(true), flavour(flavourIn) {
25 
26  // Set the pointers and category/key strings and mass splitting.
27  infoPtr = infoPtrIn;
28  settingsPtr = settingsPtrIn;
29  particleDataPtr = particleDataPtrIn;
30  cat = (flavour == 4) ? "Charmonium" : "Bottomonium";
31  key = (flavour == 4) ? "ccbar" : "bbbar";
32  mSplit = settingsPtr->parm("Onia:massSplit");
33  if (!settingsPtr->flag("Onia:forceMassSplit")) mSplit = -mSplit;
34 
35  // Set the general switch settings.
36  onia = settingsPtr->flag("Onia:all");
37  onia3S1 = settingsPtr->flag("Onia:all(3S1)");
38  onia3PJ = settingsPtr->flag("Onia:all(3PJ)");
39  onia3DJ = settingsPtr->flag("Onia:all(3DJ)");
40  oniaFlavour = settingsPtr->flag(cat + ":all");
41 
42  // Set the names of the matrix element settings.
43  meNames3S1.push_back(cat + ":O(3S1)[3S1(1)]");
44  meNames3S1.push_back(cat + ":O(3S1)[3S1(8)]");
45  meNames3S1.push_back(cat + ":O(3S1)[1S0(8)]");
46  meNames3S1.push_back(cat + ":O(3S1)[3P0(8)]");
47  meNames3PJ.push_back(cat + ":O(3PJ)[3P0(1)]");
48  meNames3PJ.push_back(cat + ":O(3PJ)[3S1(8)]");
49  meNames3DJ.push_back(cat + ":O(3DJ)[3D1(1)]");
50  meNames3DJ.push_back(cat + ":O(3DJ)[3P0(8)]");
51 
52  // Set the names of the production settings.
53  ggNames3S1.push_back(cat + ":gg2" + key + "(3S1)[3S1(1)]g");
54  ggNames3S1.push_back(cat + ":gg2" + key + "(3S1)[3S1(8)]g");
55  ggNames3S1.push_back(cat + ":gg2" + key + "(3S1)[1S0(8)]g");
56  ggNames3S1.push_back(cat + ":gg2" + key + "(3S1)[3PJ(8)]g");
57  qgNames3S1.push_back(cat + ":qg2" + key + "(3S1)[3S1(8)]q");
58  qgNames3S1.push_back(cat + ":qg2" + key + "(3S1)[1S0(8)]q");
59  qgNames3S1.push_back(cat + ":qg2" + key + "(3S1)[3PJ(8)]q");
60  qqNames3S1.push_back(cat + ":qqbar2" + key + "(3S1)[3S1(8)]g");
61  qqNames3S1.push_back(cat + ":qqbar2" + key + "(3S1)[1S0(8)]g");
62  qqNames3S1.push_back(cat + ":qqbar2" + key + "(3S1)[3PJ(8)]g");
63  ggNames3PJ.push_back(cat + ":gg2" + key + "(3PJ)[3PJ(1)]g");
64  ggNames3PJ.push_back(cat + ":gg2" + key + "(3PJ)[3S1(8)]g");
65  qgNames3PJ.push_back(cat + ":qg2" + key + "(3PJ)[3PJ(1)]q");
66  qgNames3PJ.push_back(cat + ":qg2" + key + "(3PJ)[3S1(8)]q");
67  qqNames3PJ.push_back(cat + ":qqbar2" + key + "(3PJ)[3PJ(1)]g");
68  qqNames3PJ.push_back(cat + ":qqbar2" + key + "(3PJ)[3S1(8)]g");
69  ggNames3DJ.push_back(cat + ":gg2" + key + "(3DJ)[3DJ(1)]g");
70  ggNames3DJ.push_back(cat + ":gg2" + key + "(3DJ)[3PJ(8)]g");
71  qgNames3DJ.push_back(cat + ":qg2" + key + "(3DJ)[3PJ(8)]q");
72  qqNames3DJ.push_back(cat + ":qqbar2" + key + "(3DJ)[3PJ(8)]g");
73 
74  // Initialise and check all settings.
75  states3S1 = settingsPtr->mvec(cat + ":states(3S1)");
76  initStates("3S1", states3S1, spins3S1, valid3S1);
77  initSettings("3S1", states3S1.size(), meNames3S1, mes3S1, valid3S1);
78  initSettings("3S1", states3S1.size(), ggNames3S1, ggs3S1, valid3S1);
79  initSettings("3S1", states3S1.size(), qgNames3S1, qgs3S1, valid3S1);
80  initSettings("3S1", states3S1.size(), qqNames3S1, qqs3S1, valid3S1);
81  states3PJ = settingsPtr->mvec(cat + ":states(3PJ)");
82  initStates("3PJ", states3PJ, spins3PJ, valid3PJ);
83  initSettings("3PJ", states3PJ.size(), meNames3PJ, mes3PJ, valid3PJ);
84  initSettings("3PJ", states3PJ.size(), ggNames3PJ, ggs3PJ, valid3PJ);
85  initSettings("3PJ", states3PJ.size(), qgNames3PJ, qgs3PJ, valid3PJ);
86  initSettings("3PJ", states3PJ.size(), qqNames3PJ, qqs3PJ, valid3PJ);
87  states3DJ = settingsPtr->mvec(cat + ":states(3DJ)");
88  initStates("3DJ", states3DJ, spins3DJ, valid3DJ);
89  initSettings("3DJ", states3DJ.size(), meNames3DJ, mes3DJ, valid3DJ);
90  initSettings("3DJ", states3DJ.size(), ggNames3DJ, ggs3DJ, valid3DJ);
91  initSettings("3DJ", states3DJ.size(), qgNames3DJ, qgs3DJ, valid3DJ);
92  initSettings("3DJ", states3DJ.size(), qqNames3DJ, qqs3DJ, valid3DJ);
93 
94 }
95 
96 //--------------------------------------------------------------------------
97 
98 // Initialise the SigmaProcesses for g g -> X g production.
99 
100 void SigmaOniaSetup::setupSigma2gg(vector<SigmaProcess*> &procs, bool oniaIn) {
101 
102  // Initialise the 3S1 processes.
103  if (valid3S1) {
104  for (unsigned int i = 0; i < states3S1.size(); ++i) {
105  bool flag = oniaIn || onia || onia3S1 || oniaFlavour;
106  // Colour-singlet.
107  if (flag || ggs3S1[0][i])
108  procs.push_back(new Sigma2gg2QQbar3S11g
109  (states3S1[i], mes3S1[0][i], flavour*100 + 1));
110  // Colour-octet.
111  if (flag || ggs3S1[1][i])
112  procs.push_back(new Sigma2gg2QQbarX8g
113  (states3S1[i], mes3S1[1][i], 0, mSplit, flavour*100+2));
114  if (flag || ggs3S1[2][i])
115  procs.push_back(new Sigma2gg2QQbarX8g
116  (states3S1[i], mes3S1[2][i], 1, mSplit, flavour*100+5));
117  if (flag || ggs3S1[3][i])
118  procs.push_back(new Sigma2gg2QQbarX8g
119  (states3S1[i], mes3S1[3][i], 2, mSplit, flavour*100+8));
120  }
121  }
122 
123  // Initialise the 3PJ processes.
124  if (valid3PJ) {
125  for (unsigned int i = 0; i < states3PJ.size(); ++i) {
126  bool flag = oniaIn || onia || onia3PJ || oniaFlavour;
127  // Colour-singlet.
128  if (flag || ggs3PJ[0][i]) {
129  procs.push_back(new Sigma2gg2QQbar3PJ1g
130  (states3PJ[i], mes3PJ[0][i], spins3PJ[i], flavour*100 + 11));
131  }
132  // Colour-octet.
133  if (flag || ggs3PJ[1][i])
134  procs.push_back(new Sigma2gg2QQbarX8g
135  (states3PJ[i], mes3PJ[1][i], 0, mSplit, flavour*100+14));
136  }
137  }
138 
139  // Initialise the 3DJ processes.
140  if (valid3DJ) {
141  for (unsigned int i = 0; i < states3DJ.size(); ++i) {
142  bool flag = oniaIn || onia || onia3DJ || oniaFlavour;
143  // Colour-singlet.
144  if (flag || ggs3DJ[0][i]) {
145  procs.push_back(new Sigma2gg2QQbar3DJ1g
146  (states3DJ[i], mes3DJ[0][i], spins3DJ[i], flavour*100 + 17));
147  }
148  // Colour-octet.
149  if (flag || ggs3DJ[1][i]) {
150  procs.push_back(new Sigma2gg2QQbarX8g
151  (states3DJ[i], mes3DJ[1][i], 2, mSplit, flavour*100+18));
152  }
153  }
154  }
155 
156 }
157 
158 //--------------------------------------------------------------------------
159 
160 // Initialise the SigmaProcesses for q g -> X q production.
161 
162 void SigmaOniaSetup::setupSigma2qg(vector<SigmaProcess*> &procs, bool oniaIn) {
163 
164  // Initialise the 3S1 processes.
165  if (valid3S1) {
166  for (unsigned int i = 0; i < states3S1.size(); ++i) {
167  bool flag = oniaIn || onia || onia3S1 || oniaFlavour;
168  // Colour-octet.
169  if (flag || qgs3S1[0][i])
170  procs.push_back(new Sigma2qg2QQbarX8q
171  (states3S1[i], mes3S1[1][i], 0, mSplit, flavour*100+3));
172  if (flag || qgs3S1[1][i])
173  procs.push_back(new Sigma2qg2QQbarX8q
174  (states3S1[i], mes3S1[2][i], 1, mSplit, flavour*100+6));
175  if (flag || qgs3S1[2][i])
176  procs.push_back(new Sigma2qg2QQbarX8q
177  (states3S1[i], mes3S1[3][i], 2, mSplit, flavour*100+9));
178  }
179  }
180 
181  // Initialise the 3PJ processes.
182  if (valid3PJ) {
183  for (unsigned int i = 0; i < states3PJ.size(); ++i) {
184  bool flag = oniaIn || onia || onia3PJ || oniaFlavour;
185  // Colour-singlet.
186  if (flag || qgs3PJ[0][i])
187  procs.push_back(new Sigma2qg2QQbar3PJ1q
188  (states3PJ[i], mes3PJ[0][i], spins3PJ[i], flavour*100 + 12));
189  // Colour-octet.
190  if (flag || qgs3PJ[1][i])
191  procs.push_back(new Sigma2qg2QQbarX8q
192  (states3PJ[i], mes3PJ[1][i], 0, mSplit, flavour*100+15));
193  }
194  }
195 
196  // Initialise the 3DJ processes.
197  if (valid3DJ) {
198  for (unsigned int i = 0; i < states3DJ.size(); ++i) {
199  bool flag = oniaIn || onia || onia3DJ || oniaFlavour;
200  // Colour-octet.
201  if (flag || qgs3DJ[0][i])
202  procs.push_back(new Sigma2qg2QQbarX8q
203  (states3DJ[i], mes3DJ[1][i], 2, mSplit, flavour*100+19));
204  }
205  }
206 
207 }
208 
209 //--------------------------------------------------------------------------
210 
211 // Initialise the SigmaProcesses for q qbar -> X g production.
212 
213 void SigmaOniaSetup::setupSigma2qq(vector<SigmaProcess*> &procs, bool oniaIn) {
214 
215  // Initialise the 3S1 processes.
216  if (valid3S1) {
217  for (unsigned int i = 0; i < states3S1.size(); ++i) {
218  bool flag = oniaIn || onia || onia3S1 || oniaFlavour;
219  // Colour-octet.
220  if (flag || qqs3S1[0][i])
221  procs.push_back(new Sigma2qqbar2QQbarX8g
222  (states3S1[i], mes3S1[1][i], 0, mSplit, flavour*100+4));
223  if (flag || qqs3S1[1][i])
224  procs.push_back(new Sigma2qqbar2QQbarX8g
225  (states3S1[i], mes3S1[2][i], 1, mSplit, flavour*100+7));
226  if (flag || qqs3S1[2][i])
227  procs.push_back(new Sigma2qqbar2QQbarX8g
228  (states3S1[i], mes3S1[3][i], 2, mSplit, flavour*100+10));
229  }
230  }
231 
232  // Initialise the 3PJ processes.
233  if (valid3PJ) {
234  for (unsigned int i = 0; i < states3PJ.size(); ++i) {
235  bool flag = oniaIn || onia || onia3PJ || oniaFlavour;
236  // Colour-singlet.
237  if (flag || qqs3PJ[0][i])
238  procs.push_back(new Sigma2qqbar2QQbar3PJ1g
239  (states3PJ[i], mes3PJ[0][i], spins3PJ[i], flavour*100 + 13));
240  // Colour-octet.
241  if (flag || qqs3PJ[1][i])
242  procs.push_back(new Sigma2qqbar2QQbarX8g
243  (states3PJ[i], mes3PJ[1][i], 0, mSplit, flavour*100+16));
244  }
245  }
246 
247  // Initialise the 3DJ processes.
248  if (valid3DJ) {
249  for (unsigned int i = 0; i < states3DJ.size(); ++i) {
250  bool flag = oniaIn || onia || onia3DJ || oniaFlavour;
251  // Colour-octet.
252  if (flag || qqs3DJ[0][i])
253  procs.push_back(new Sigma2qqbar2QQbarX8g
254  (states3DJ[i], mes3DJ[1][i], 2, mSplit, flavour*100+20));
255  }
256  }
257 
258 }
259 
260 //--------------------------------------------------------------------------
261 
262 // Initialise and check the flavour, j-number, and validity of states.
263 
264 void SigmaOniaSetup::initStates(string wave, const vector<int> &states,
265  vector<int> &jnums, bool &valid) {
266 
267  set<int> unique;
268  unsigned int nstates(0);
269  for (unsigned int i = 0; i < states.size(); ++i) {
270 
271  // Check state is unique and remove if not.
272  stringstream state;
273  state << states[i];
274  unique.insert(states[i]);
275  if (nstates + 1 != unique.size()) {
276  infoPtr->errorMsg("Error in SigmaOniaSetup::initStates: particle "
277  + state.str() + " in mvec " + cat + ":states("
278  + wave + ")", "has duplicates");
279  valid = false;
280  } else ++nstates;
281 
282  // Determine quark composition and quantum numbers.
283  int mod1(10), mod2(1);
284  vector<int> digits;
285  while (digits.size() < 7) {
286  digits.push_back((states[i]%mod1 - states[i]%mod2) / mod2);
287  mod1 *= 10;
288  mod2 *= 10;
289  }
290  int s, l, j((digits[0] - 1)/2);
291  if (j != 0) {
292  if (digits[4] == 0) {l = j - 1; s = 1;}
293  else if (digits[4] == 1) {l = j; s = 0;}
294  else if (digits[4] == 2) {l = j; s = 1;}
295  else {l = j + 1; s = 1;}
296  } else {
297  if (digits[4] == 0) {l = 0; s = 0;}
298  else {l = 1; s = 1;}
299  }
300 
301  // Check state validity.
302  if (states[i] != 0) {
303  if (!particleDataPtr->isParticle(states[i])) {
304  infoPtr->errorMsg("Error in SigmaOniaSetup::initStates: particle "
305  + state.str() + " in mvec " + cat + ":states("
306  + wave + ")", "is unknown");
307  valid = false;
308  }
309  if (digits[3] != 0) {
310  infoPtr->errorMsg("Error in SigmaOniaSetup::initStates: particle "
311  + state.str() + " in mvec " + cat + ":states("
312  + wave + ")", " is not a meson");
313  valid = false;
314  }
315  if (digits[2] != digits[1] || digits[1] != flavour) {
316  infoPtr->errorMsg("Error in SigmaOniaSetup::initStates: particle "
317  + state.str() + " in mvec " + cat + ":states("
318  + wave + ")", "is not a " + key + " state");
319  valid = false;
320  }
321  if ((wave == "3S1" && (s != 1 || l != 0 || j != 1)) ||
322  (wave == "3PJ" && (s != 1 || l != 1 || j < 0 || j > 2)) ||
323  (wave == "3DJ" && (s != 1 || l != 2 || j < 1 || j > 3))) {
324  infoPtr->errorMsg("Error in SigmaOniaSetup::initStates: particle "
325  + state.str() + " in mvec " + cat + ":states("
326  + wave + ")", "is not a " + wave + " state");
327  valid = false;
328  }
329  } else valid = false;
330  jnums.push_back(j);
331  }
332 
333 }
334 
335 //--------------------------------------------------------------------------
336 
337 // Initialise and check a group of PVec settings.
338 
339 void SigmaOniaSetup::initSettings(string wave, unsigned int size,
340  const vector<string> &names, vector< vector<double> > &pvecs,
341  bool &valid) {
342 
343  for (unsigned int i = 0; i < names.size(); ++i) {
344  pvecs.push_back(settingsPtr->pvec(names[i]));
345  if (pvecs.back().size() != size) {
346  infoPtr->errorMsg("Error in SigmaOniaSetup::initSettings: mvec " + cat
347  + ":states(" + wave + ")", "is not the same size as"
348  " pvec " + names[i]);
349  valid = false;
350  }
351  }
352 
353 }
354 
355 //--------------------------------------------------------------------------
356 
357 // Initialise and check a group of FVec settings.
358 
359 void SigmaOniaSetup::initSettings(string wave, unsigned int size,
360  const vector<string> &names, vector< vector<bool> > &fvecs,
361  bool &valid) {
362 
363  for (unsigned int i = 0; i < names.size(); ++i) {
364  fvecs.push_back(settingsPtr->fvec(names[i]));
365  if (fvecs.back().size() != size) {
366  infoPtr->errorMsg("Error in SigmaOniaSetup::initSettings: mvec " + cat
367  + ":states(" + wave + ")", "is not the same size as"
368  " fvec " + names[i]);
369  valid = false;
370  }
371  }
372 
373 }
374 
375 //==========================================================================
376 
377 // Sigma2gg2QQbar3S11g class.
378 // Cross section g g -> QQbar[3S1(1)] g (Q = c or b).
379 
380 //--------------------------------------------------------------------------
381 
382 // Initialize process.
383 
384 void Sigma2gg2QQbar3S11g::initProc() {
385 
386  // Process name.
387  nameSave = "g g -> "
388  + string((codeSave - codeSave%100)/100 == 4 ? "ccbar" : "bbbar")
389  + "(3S1)[3S1(1)] g";
390 
391 }
392 
393 //--------------------------------------------------------------------------
394 
395 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
396 
397 void Sigma2gg2QQbar3S11g::sigmaKin() {
398 
399  // Calculate kinematics dependence.
400  double stH = sH + tH;
401  double tuH = tH + uH;
402  double usH = uH + sH;
403  double sig = (10. * M_PI / 81.) * m3 * ( pow2(sH * tuH)
404  + pow2(tH * usH) + pow2(uH * stH) ) / pow2( stH * tuH * usH );
405 
406  // Answer.
407  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
408 
409 }
410 
411 //--------------------------------------------------------------------------
412 
413 // Select identity, colour and anticolour.
414 
415 void Sigma2gg2QQbar3S11g::setIdColAcol() {
416 
417  // Flavours are trivial.
418  setId( id1, id2, idHad, 21);
419 
420  // Two orientations of colour flow.
421  setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
422  if (rndmPtr->flat() > 0.5) swapColAcol();
423 
424 }
425 
426 //==========================================================================
427 
428 // Sigma2gg2QQbar3PJ1g class.
429 // Cross section g g -> QQbar[3PJ(1)] g (Q = c or b, J = 0, 1 or 2).
430 
431 //--------------------------------------------------------------------------
432 
433 // Initialize process.
434 
435 void Sigma2gg2QQbar3PJ1g::initProc() {
436 
437  // Process name.
438  if (jSave >= 0 && jSave <= 2)
439  nameSave = namePrefix() + " -> " + nameMidfix() + "(3PJ)[3PJ(1)] "
440  + namePostfix();
441  else
442  nameSave = "illegal process";
443 
444 }
445 
446 //--------------------------------------------------------------------------
447 
448 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
449 
450 void Sigma2gg2QQbar3PJ1g::sigmaKin() {
451 
452  // Useful derived kinematics quantities.
453  double pRat = (sH * uH + uH * tH + tH * sH)/ sH2;
454  double qRat = tH * uH / sH2;
455  double rRat = s3 / sH;
456  double pRat2 = pRat * pRat;
457  double pRat3 = pRat2 * pRat;
458  double pRat4 = pRat3 * pRat;
459  double qRat2 = qRat * qRat;
460  double qRat3 = qRat2 * qRat;
461  double qRat4 = qRat3 * qRat;
462  double rRat2 = rRat * rRat;
463  double rRat3 = rRat2 * rRat;
464  double rRat4 = rRat3 * rRat;
465 
466  // Calculate kinematics dependence.
467  double sig = 0.;
468  if (jSave == 0) {
469  sig = (8. * M_PI / (9. * m3 * sH))
470  * ( 9. * rRat2 * pRat4 * (rRat4 - 2. * rRat2 * pRat + pRat2)
471  - 6. * rRat * pRat3 * qRat * (2. * rRat4 - 5. * rRat2 * pRat
472  + pRat2) - pRat2 * qRat2 * (rRat4 + 2. * rRat2 * pRat - pRat2)
473  + 2. * rRat * pRat * qRat3 * (rRat2 - pRat) + 6. * rRat2 * qRat4)
474  / (qRat * pow4(qRat - rRat * pRat));
475  } else if (jSave == 1) {
476  sig = (8. * M_PI / (3.* m3 * sH)) * pRat2
477  * (rRat * pRat2 * (rRat2 - 4. * pRat)
478  + 2. * qRat * (-rRat4 + 5. * rRat2 * pRat + pRat2)
479  - 15. * rRat * qRat2) / pow4(qRat - rRat * pRat);
480  } else if (jSave == 2) {
481  sig = (8. * M_PI / (9. * m3 * sH))
482  * (12. * rRat2 * pRat4 * (rRat4 - 2. * rRat2 * pRat + pRat2)
483  - 3. * rRat * pRat3 * qRat * (8. * rRat4 - rRat2 * pRat + 4. * pRat2)
484  + 2. * pRat2 * qRat2 * (-7. * rRat4 + 43. * rRat2 * pRat + pRat2)
485  + rRat * pRat * qRat3 * (16. * rRat2 - 61. * pRat)
486  + 12. * rRat2 * qRat4) / (qRat * pow4(qRat-rRat * pRat));
487  }
488 
489  // Answer.
490  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
491 
492 }
493 
494 //--------------------------------------------------------------------------
495 
496 // Select identity, colour and anticolour.
497 
498 void Sigma2gg2QQbar3PJ1g::setIdColAcol() {
499 
500  // Flavours are trivial.
501  setId( id1, id2, idHad, 21);
502 
503  // Two orientations of colour flow.
504  setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
505  if (rndmPtr->flat() > 0.5) swapColAcol();
506 
507 }
508 
509 //==========================================================================
510 
511 // Sigma2qg2QQbar3PJ1q class.
512 // Cross section q g -> QQbar[3PJ(1)] q (Q = c or b, J = 0, 1 or 2).
513 
514 //--------------------------------------------------------------------------
515 
516 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
517 
518 void Sigma2qg2QQbar3PJ1q::sigmaKin() {
519 
520  // Calculate kinematics dependence.
521  double usH = uH + sH;
522  double sig = 0.;
523  if (jSave == 0) {
524  sig = - (16. * M_PI / 81.) * pow2(tH - 3. * s3) * (sH2 + uH2)
525  / (m3 * tH * pow4(usH));
526  } else if (jSave == 1) {
527  sig = - (32. * M_PI / 27.) * (4. * s3 * sH * uH + tH * (sH2 + uH2))
528  / (m3 * pow4(usH));
529  } else if (jSave == 2) {
530  sig = - (32. *M_PI / 81.) * ( (6. * s3*s3 + tH2) * pow2(usH)
531  - 2. * sH * uH * (tH2 + 6. * s3 * usH)) / (m3 * tH * pow4(usH));
532  }
533 
534  // Answer.
535  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
536 
537 }
538 
539 //--------------------------------------------------------------------------
540 
541 // Select identity, colour and anticolour.
542 
543 void Sigma2qg2QQbar3PJ1q::setIdColAcol() {
544 
545  // Flavours are trivial.
546  int idq = (id2 == 21) ? id1 : id2;
547  setId( id1, id2, idHad, idq);
548 
549  // tH defined between q_in and q_out: must swap tHat <-> uHat if q g in.
550  swapTU = (id2 == 21);
551 
552  // Colour flow topologies. Swap when antiquarks.
553  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
554  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
555  if (idq < 0) swapColAcol();
556 
557 }
558 
559 //==========================================================================
560 
561 // Sigma2qqbar2QQbar3PJ1g class.
562 // Cross section q qbar -> QQbar[3PJ(1)] g (Q = c or b, J = 0, 1 or 2).
563 
564 //--------------------------------------------------------------------------
565 
566 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
567 
568 void Sigma2qqbar2QQbar3PJ1g::sigmaKin() {
569 
570  // Calculate kinematics dependence.
571  double tuH = tH + uH;
572  double sig = 0.;
573  if (jSave == 0) {
574  sig =(128. * M_PI / 243.) * pow2(sH - 3. * s3) * (tH2 + uH2)
575  / (m3 * sH * pow4(tuH));
576  } else if (jSave == 1) {
577  sig = (256. * M_PI / 81.) * (4. * s3 * tH * uH + sH * (tH2 + uH2))
578  / (m3 * pow4(tuH));
579  } else if (jSave == 2) {
580  sig = (256. * M_PI / 243.) * ( (6. * s3*s3 + sH2) * pow2(tuH)
581  - 2. * tH * uH * (sH2 + 6. * s3 * tuH) )/ (m3 * sH * pow4(tuH));
582  }
583 
584  // Answer.
585  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
586 
587 }
588 
589 //--------------------------------------------------------------------------
590 
591 // Select identity, colour and anticolour.
592 
593 void Sigma2qqbar2QQbar3PJ1g::setIdColAcol() {
594 
595  // Flavours are trivial.
596  setId( id1, id2, idHad, 21);
597 
598  // Colour flow topologies. Swap when antiquarks.
599  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
600  if (id1 < 0) swapColAcol();
601 
602 }
603 
604 //==========================================================================
605 
606 // Sigma2gg2QQbar3DJ1g class.
607 // Cross section g g -> QQbar[3DJ(1)] g (Q = c or b).
608 
609 //--------------------------------------------------------------------------
610 
611 // Initialize process.
612 
613 void Sigma2gg2QQbar3DJ1g::initProc() {
614 
615  // Process name.
616  if (jSave >= 1 && jSave <= 3)
617  nameSave = namePrefix() + " -> " + nameMidfix() + "(3DJ)[3DJ(1)] "
618  + namePostfix();
619  else
620  nameSave = "illegal process";
621 
622 }
623 
624 //--------------------------------------------------------------------------
625 
626 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
627 
628 void Sigma2gg2QQbar3DJ1g::sigmaKin() {
629 
630  // Calculate kinematics dependence.
631  double m2V[12], sHV[12], mpsV[8], mmsV[6], mmtV[6], sptV[6];
632  m2V[0] = 1;
633  sHV[0] = 1;
634  mmtV[0] = 1;
635  mpsV[0] = 1;
636  mmsV[0] = 1;
637  sptV[0] = 1;
638  for (int i = 1; i < 12; ++i) {
639  m2V[i] = m2V[i - 1] * s3;
640  sHV[i] = sHV[i - 1] * sH;
641  if (i < 8) {
642  mpsV[i] = mpsV[i - 1] * (s3 + sH);
643  if (i < 6) {
644  mmsV[i] = mmsV[i - 1] * (s3 - sH);
645  mmtV[i] = mmtV[i - 1] * (s3 - tH);
646  sptV[i] = sptV[i - 1] * (sH + tH);
647  }
648  }
649  }
650  double fac = (pow3(alpS)*pow2(M_PI));
651  double sig = 0;
652  if (jSave == 1) {
653  fac *= 16. / 81.;
654  sig = -25/(sqrt(m2V[1])*mmsV[5]) + (49*sqrt(m2V[3]))/(mmsV[5]*sHV[2])
655  + (48*sqrt(m2V[3])*sHV[2]*(m2V[2] + sHV[2]))/(mmsV[3]*mmtV[5]*mpsV[3])
656  - (67*sqrt(m2V[1]))/(mmsV[5]*sHV[1]) - (5*sHV[1])/(sqrt(m2V[3])*mmsV[5])
657  + (4*sqrt(m2V[1])*(m2V[6] + 97*m2V[4]*sHV[2] - 48*m2V[3]*sHV[3]
658  + 105*m2V[2]*sHV[4] + 33*sHV[6] -
659  24*m2V[5]*sHV[1]))/(mmsV[4]*mmtV[4]*mpsV[4]) - (4*(m2V[9] +
660  197*m2V[7]*sHV[2] - 50*m2V[6]*sHV[3] + 509*m2V[5]*sHV[4] -
661  416*m2V[4]*sHV[5] + 237*m2V[3]*sHV[6] - 400*m2V[2]*sHV[7] -
662  10*sHV[9] -
663  164*m2V[8]*sHV[1]))/(sqrt(m2V[1])*mmsV[5]*mmtV[3]*mpsV[5]*sHV[1]) +
664  (224*m2V[10] + 1825*m2V[8]*sHV[2] - 3980*m2V[7]*sHV[3] +
665  3996*m2V[6]*sHV[4] - 4766*m2V[5]*sHV[5] + 10022*m2V[4]*sHV[6] -
666  5212*m2V[3]*sHV[7] + 6124*m2V[2]*sHV[8] - 869*m2V[1]*sHV[9] +
667  145*sHV[10] -
668  597*m2V[9]*sHV[1])/(sqrt(m2V[1])*mmsV[5]*mmtV[1]*mpsV[7]*sHV[2]) +
669  (102*m2V[11] + 331*m2V[9]*sHV[2] - 2021*m2V[8]*sHV[3] +
670  3616*m2V[7]*sHV[4] - 968*m2V[6]*sHV[5] + 3386*m2V[5]*sHV[6] -
671  6150*m2V[4]*sHV[7] + 666*m2V[3]*sHV[8] - 1134*m2V[2]*sHV[9] -
672  5*m2V[1]*sHV[10] - 5*sHV[11] -
673  506*m2V[10]*sHV[1])/(sqrt(m2V[3])*mmsV[5]*mmtV[2]*mpsV[6]*sHV[2]) +
674  (48*sqrt(m2V[3])*sHV[2]*(m2V[2] + sHV[2]))/(mmsV[3]*mpsV[3]*sptV[5])
675  + (4*sqrt(m2V[1])*(m2V[6] + 97*m2V[4]*sHV[2] - 48*m2V[3]*sHV[3] +
676  105*m2V[2]*sHV[4] + 33*sHV[6] -
677  24*m2V[5]*sHV[1]))/(mmsV[4]*mpsV[4]*sptV[4]) - (4*(m2V[9] +
678  197*m2V[7]*sHV[2] - 50*m2V[6]*sHV[3] + 509*m2V[5]*sHV[4] -
679  416*m2V[4]*sHV[5] + 237*m2V[3]*sHV[6] - 400*m2V[2]*sHV[7] -
680  10*sHV[9] -
681  164*m2V[8]*sHV[1]))/(sqrt(m2V[1])*mmsV[5]*mpsV[5]*sHV[1]*sptV[3]) +
682  (102*m2V[11] + 331*m2V[9]*sHV[2] - 2021*m2V[8]*sHV[3] +
683  3616*m2V[7]*sHV[4] - 968*m2V[6]*sHV[5] + 3386*m2V[5]*sHV[6] -
684  6150*m2V[4]*sHV[7] + 666*m2V[3]*sHV[8] - 1134*m2V[2]*sHV[9] -
685  5*m2V[1]*sHV[10] - 5*sHV[11] -
686  506*m2V[10]*sHV[1])/(sqrt(m2V[3])*mmsV[5]*mpsV[6]*sHV[2]*sptV[2]) +
687  (224*m2V[10] + 1825*m2V[8]*sHV[2] - 3980*m2V[7]*sHV[3] +
688  3996*m2V[6]*sHV[4] - 4766*m2V[5]*sHV[5] + 10022*m2V[4]*sHV[6] -
689  5212*m2V[3]*sHV[7] + 6124*m2V[2]*sHV[8] - 869*m2V[1]*sHV[9] +
690  145*sHV[10] -
691  597*m2V[9]*sHV[1])/(sqrt(m2V[1])*mmsV[5]*mpsV[7]*sHV[2]*sptV[1]);
692  } else if (jSave == 2) {
693  fac *= 32. / 27.;
694  sig = 16/(sqrt(m2V[1])*mmsV[5]) +
695  (2*sqrt(m2V[3]))/(mmsV[5]*sHV[2]) - (8*sqrt(m2V[3])*sHV[2]*(m2V[2] +
696  sHV[2]))/(mmsV[3]*mmtV[5]*mpsV[3]) +
697  (6*sqrt(m2V[1]))/(mmsV[5]*sHV[1]) -
698  (16*sHV[1])/(sqrt(m2V[3])*mmsV[5]) - (2*sqrt(m2V[1])*(3*m2V[6] -
699  25*m2V[4]*sHV[2] - 16*m2V[3]*sHV[3] - 33*m2V[2]*sHV[4] - 5*sHV[6] -
700  8*m2V[5]*sHV[1]))/(mmsV[4]*mmtV[4]*mpsV[4]) + (2*(3*m2V[9] -
701  41*m2V[7]*sHV[2] - 37*m2V[6]*sHV[3] - 149*m2V[5]*sHV[4] +
702  55*m2V[4]*sHV[5] - 53*m2V[3]*sHV[6] + 167*m2V[2]*sHV[7] + 16*sHV[9]
703  + 7*m2V[8]*sHV[1]))/(sqrt(m2V[1])*mmsV[5]*mmtV[3]*mpsV[5]*sHV[1]) +
704  (2*(m2V[10] + 34*m2V[8]*sHV[2] - 198*m2V[7]*sHV[3] -
705  140*m2V[6]*sHV[4] - 746*m2V[5]*sHV[5] + 226*m2V[4]*sHV[6] -
706  486*m2V[3]*sHV[7] + 679*m2V[2]*sHV[8] - 50*m2V[1]*sHV[9] +
707  112*sHV[10] -
708  8*m2V[9]*sHV[1]))/(sqrt(m2V[1])*mmsV[5]*mmtV[1]*mpsV[7]*sHV[2]) +
709  (m2V[11] + 19*m2V[9]*sHV[2] - m2V[8]*sHV[3] + 597*m2V[7]*sHV[4] +
710  321*m2V[6]*sHV[5] + 797*m2V[5]*sHV[6] - 791*m2V[4]*sHV[7] +
711  26*m2V[3]*sHV[8] - 468*m2V[2]*sHV[9] - 16*m2V[1]*sHV[10] -
712  16*sHV[11] -
713  21*m2V[10]*sHV[1])/(sqrt(m2V[3])*mmsV[5]*mmtV[2]*mpsV[6]*sHV[2]) -
714  (8*sqrt(m2V[3])*sHV[2]*(m2V[2] + sHV[2]))/(mmsV[3]*mpsV[3]*sptV[5])
715  - (2*sqrt(m2V[1])*(3*m2V[6] - 25*m2V[4]*sHV[2] - 16*m2V[3]*sHV[3] -
716  33*m2V[2]*sHV[4] - 5*sHV[6] -
717  8*m2V[5]*sHV[1]))/(mmsV[4]*mpsV[4]*sptV[4]) + (2*(3*m2V[9] -
718  41*m2V[7]*sHV[2] - 37*m2V[6]*sHV[3] - 149*m2V[5]*sHV[4] +
719  55*m2V[4]*sHV[5] - 53*m2V[3]*sHV[6] + 167*m2V[2]*sHV[7] + 16*sHV[9]
720  + 7*m2V[8]*sHV[1]))/(sqrt(m2V[1])*mmsV[5]*mpsV[5]*sHV[1]*sptV[3]) +
721  (m2V[11] + 19*m2V[9]*sHV[2] - m2V[8]*sHV[3] + 597*m2V[7]*sHV[4] +
722  321*m2V[6]*sHV[5] + 797*m2V[5]*sHV[6] - 791*m2V[4]*sHV[7] +
723  26*m2V[3]*sHV[8] - 468*m2V[2]*sHV[9] - 16*m2V[1]*sHV[10] -
724  16*sHV[11] -
725  21*m2V[10]*sHV[1])/(sqrt(m2V[3])*mmsV[5]*mpsV[6]*sHV[2]*sptV[2]) +
726  (2*(m2V[10] + 34*m2V[8]*sHV[2] - 198*m2V[7]*sHV[3] -
727  140*m2V[6]*sHV[4] - 746*m2V[5]*sHV[5] + 226*m2V[4]*sHV[6] -
728  486*m2V[3]*sHV[7] + 679*m2V[2]*sHV[8] - 50*m2V[1]*sHV[9] +
729  112*sHV[10] -
730  8*m2V[9]*sHV[1]))/(sqrt(m2V[1])*mmsV[5]*mpsV[7]*sHV[2]*sptV[1]);
731  } else if (jSave == 3) {
732  fac *= 256. / 189.;
733  sig = 5/(sqrt(m2V[1])*mmsV[5]) + sqrt(m2V[3])/(mmsV[5]*sHV[2]) +
734  (2*sqrt(m2V[3])*sHV[2]*(m2V[2] + sHV[2]))/(mmsV[3]*mmtV[5]*mpsV[3])
735  - (3*sqrt(m2V[1]))/(mmsV[5]*sHV[1]) -
736  (5*sHV[1])/(sqrt(m2V[3])*mmsV[5]) + (sqrt(m2V[1])*(6*m2V[6] +
737  67*m2V[4]*sHV[2] - 8*m2V[3]*sHV[3] + 45*m2V[2]*sHV[4] + 8*sHV[6] -
738  4*m2V[5]*sHV[1]))/(mmsV[4]*mmtV[4]*mpsV[4]) + (-6*m2V[9] -
739  152*m2V[7]*sHV[2] + 80*m2V[6]*sHV[3] - 269*m2V[5]*sHV[4] +
740  211*m2V[4]*sHV[5] - 77*m2V[3]*sHV[6] + 155*m2V[2]*sHV[7] + 10*sHV[9]
741  + 64*m2V[8]*sHV[1])/(sqrt(m2V[1])*mmsV[5]*mmtV[3]*mpsV[5]*sHV[1]) +
742  (16*m2V[10] + 295*m2V[8]*sHV[2] - 555*m2V[7]*sHV[3] +
743  769*m2V[6]*sHV[4] - 1079*m2V[5]*sHV[5] + 913*m2V[4]*sHV[6] -
744  603*m2V[3]*sHV[7] + 601*m2V[2]*sHV[8] - 56*m2V[1]*sHV[9] +
745  70*sHV[10] -
746  83*m2V[9]*sHV[1])/(sqrt(m2V[1])*mmsV[5]*mmtV[1]*mpsV[7]*sHV[2]) +
747  (8*m2V[11] + 104*m2V[9]*sHV[2] - 284*m2V[8]*sHV[3] +
748  549*m2V[7]*sHV[4] - 282*m2V[6]*sHV[5] + 514*m2V[5]*sHV[6] -
749  520*m2V[4]*sHV[7] + 34*m2V[3]*sHV[8] - 171*m2V[2]*sHV[9] -
750  5*m2V[1]*sHV[10] - 5*sHV[11] -
751  54*m2V[10]*sHV[1])/(sqrt(m2V[3])*mmsV[5]*mmtV[2]*mpsV[6]*sHV[2]) +
752  (2*sqrt(m2V[3])*sHV[2]*(m2V[2] + sHV[2]))/(mmsV[3]*mpsV[3]*sptV[5])
753  + (sqrt(m2V[1])*(6*m2V[6] + 67*m2V[4]*sHV[2] - 8*m2V[3]*sHV[3] +
754  45*m2V[2]*sHV[4] + 8*sHV[6] -
755  4*m2V[5]*sHV[1]))/(mmsV[4]*mpsV[4]*sptV[4]) + (-6*m2V[9] -
756  152*m2V[7]*sHV[2] + 80*m2V[6]*sHV[3] - 269*m2V[5]*sHV[4] +
757  211*m2V[4]*sHV[5] - 77*m2V[3]*sHV[6] + 155*m2V[2]*sHV[7] + 10*sHV[9]
758  + 64*m2V[8]*sHV[1])/(sqrt(m2V[1])*mmsV[5]*mpsV[5]*sHV[1]*sptV[3]) +
759  (8*m2V[11] + 104*m2V[9]*sHV[2] - 284*m2V[8]*sHV[3] +
760  549*m2V[7]*sHV[4] - 282*m2V[6]*sHV[5] + 514*m2V[5]*sHV[6] -
761  520*m2V[4]*sHV[7] + 34*m2V[3]*sHV[8] - 171*m2V[2]*sHV[9] -
762  5*m2V[1]*sHV[10] - 5*sHV[11] -
763  54*m2V[10]*sHV[1])/(sqrt(m2V[3])*mmsV[5]*mpsV[6]*sHV[2]*sptV[2]) +
764  (16*m2V[10] + 295*m2V[8]*sHV[2] - 555*m2V[7]*sHV[3] +
765  769*m2V[6]*sHV[4] - 1079*m2V[5]*sHV[5] + 913*m2V[4]*sHV[6] -
766  603*m2V[3]*sHV[7] + 601*m2V[2]*sHV[8] - 56*m2V[1]*sHV[9] +
767  70*sHV[10] -
768  83*m2V[9]*sHV[1])/(sqrt(m2V[1])*mmsV[5]*mpsV[7]*sHV[2]*sptV[1]);
769  }
770 
771  // Answer.
772  sigma = ((2.*jSave + 1.) / 3.) * oniumME * fac * sig;
773 
774 }
775 
776 //==========================================================================
777 
778 // Sigma2gg2QQbarX8g class.
779 // Cross section g g -> QQbar[X(8)] g (Q = c or b, X = 3S1, 1S0 or 3PJ).
780 
781 //--------------------------------------------------------------------------
782 
783 // Initialize process.
784 
785 void Sigma2gg2QQbarX8g::initProc() {
786 
787  // Return for illegal process.
788  if (stateSave < 0 || stateSave > 2) {
789  idHad = 0;
790  nameSave = "illegal process";
791  return;
792  }
793 
794  // Determine quark composition and quantum numbers.
795  int mod1(10), mod2(1);
796  vector<int> digits;
797  while (digits.size() < 7) {
798  digits.push_back((idHad%mod1 - idHad%mod2) / mod2);
799  mod1 *= 10;
800  mod2 *= 10;
801  }
802  int s, l, j((digits[0] - 1)/2);
803  if (j != 0) {
804  if (digits[4] == 0) {l = j - 1; s = 1;}
805  else if (digits[4] == 1) {l = j; s = 0;}
806  else if (digits[4] == 2) {l = j; s = 1;}
807  else {l = j + 1; s = 1;}
808  } else {
809  if (digits[4] == 0) {l = 0; s = 0;}
810  else {l = 1; s = 1;}
811  }
812 
813  // Set the process name.
814  stringstream sName, jName;
815  string lName, stateName;
816  sName << 2*s + 1;
817  if (l == 0) jName << j;
818  else jName << "J";
819  if (l == 0) lName = "S";
820  else if (l == 1) lName = "P";
821  else if (l == 2) lName = "D";
822  if (stateSave == 0) stateName = "[3S1(8)]";
823  else if (stateSave == 1) stateName = "[1S0(8)]";
824  else if (stateSave == 2) stateName = "[3PJ(8)]";
825  nameSave = namePrefix() + " -> " + (digits[1] == 4 ? "ccbar" : "bbbar")
826  + "(" + sName.str() + lName + jName.str() + ")" + stateName
827  + " " + namePostfix();
828 
829  // Ensure the dummy particle for the colour-octet state is valid.
830  int idOct = 9900000 + digits[1]*10000 + stateSave*1000 + digits[5]*100
831  + digits[4]*10 + digits[0];
832  double m0 = particleDataPtr->m0(idHad) + abs(mSplit);
833  double mWidth = 0.0;
834  if (!particleDataPtr->isParticle(idOct)) {
835  string nameOct = particleDataPtr->name(idHad) + stateName;
836  int spinType = stateSave == 1 ? 1 : 3;
837  int chargeType = particleDataPtr->chargeType(idHad);
838  int colType = 2;
839  particleDataPtr->addParticle(idOct, nameOct, spinType, chargeType, colType,
840  m0, mWidth, m0, m0);
841  ParticleDataEntry* entry = particleDataPtr->particleDataEntryPtr(idOct);
842  if (entry) entry->addChannel(1, 1.0, 0, idHad, 21);
843  } else if (mSplit > 0 && abs(particleDataPtr->m0(idOct) - m0) > 1E-15) {
844  particleDataPtr->m0(idOct, m0);
845  particleDataPtr->mWidth(idOct, mWidth);
846  particleDataPtr->mMin(idOct, m0);
847  particleDataPtr->mMax(idOct, m0);
848  } else if (particleDataPtr->m0(idOct) <= particleDataPtr->m0(idHad)) {
849  infoPtr->errorMsg("Warning in Sigma2gg2QQbarX8g::initProc: mass of "
850  "intermediate colour-octet state"
851  "increased to be greater than the physical state");
852  particleDataPtr->m0(idOct, m0);
853  particleDataPtr->mWidth(idOct, mWidth);
854  particleDataPtr->mMin(idOct, m0);
855  particleDataPtr->mMax(idOct, m0);
856  }
857  idHad = idOct;
858 
859 }
860 
861 //--------------------------------------------------------------------------
862 
863 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
864 
865 void Sigma2gg2QQbarX8g::sigmaKin() {
866 
867  // Calculate kinematics dependence.
868  double stH = sH + tH;
869  double tuH = tH + uH;
870  double usH = uH + sH;
871  double sig = 0.;
872  if (stateSave == 0) {
873  sig = (M_PI / 72.) * m3 * ( 27. * (pow2(stH) + pow2(tuH)
874  + pow2(usH)) / (s3*s3) - 16. ) * ( pow2(sH * tuH)
875  + pow2(tH * usH) + pow2(uH * stH) ) / pow2( stH * tuH * usH );
876  } else if (stateSave == 1) {
877  sig = (5. * M_PI / 16.) * m3 * ( pow2(uH / (tuH * usH))
878  + pow2(sH / (stH * usH)) + pow2(tH / (stH * tuH)) ) * ( 12.
879  + (pow4(stH) + pow4(tuH) + pow4(usH)) / (s3 * sH * tH * uH) );
880  } else if (stateSave == 2) {
881  double sH3 = sH2 * sH;
882  double sH4 = sH3 * sH;
883  double sH5 = sH4 * sH;
884  double sH6 = sH5 * sH;
885  double sH7 = sH6 * sH;
886  double sH8 = sH7 * sH;
887  double tH3 = tH2 * tH;
888  double tH4 = tH3 * tH;
889  double tH5 = tH4 * tH;
890  double tH6 = tH5 * tH;
891  double tH7 = tH6 * tH;
892  double tH8 = tH7 * tH;
893  double ssttH = sH * sH + sH * tH + tH * tH;
894  sig = 5. * M_PI * (3. * sH * tH * stH * pow4(ssttH)
895  - s3 * pow2(ssttH) * (7. * sH6 + 36. * sH5 * tH + 45. * sH4 * tH2
896  + 28. * sH3 * tH3 + 45. * sH2 * tH4 + 36. * sH * tH5 + 7. * tH6)
897  + pow2(s3) * stH * (35. *sH8 + 169. * sH7 * tH + 299. * sH6 * tH2
898  + 401. * sH5 * tH3 + 418. * sH4 * tH4 + 401. * sH3 * tH5
899  + 299. * sH2 * tH6 + 169. * sH * tH7 + 35. * tH8)
900  - pow3(s3) * (84. *sH8+432. *sH7*tH+905. *sH6*tH2
901  + 1287. * sH5 * tH3 + 1436. * sH4 * tH4 +1287. * sH3 * tH5
902  + 905. * sH2 * tH6 + 432. * sH * tH7 + 84. * tH8)
903  + pow4(s3) * stH * (126. * sH6 + 451. * sH5 * tH +677. * sH4 * tH2
904  + 836. * sH3 * tH3 + 677. * sH2 * tH4 + 451. * sH * tH5
905  + 126. * tH6)
906  - pow5(s3) * 3. * (42. * sH6 + 171. * sH5 * tH + 304. * sH4 * tH2
907  + 362. * sH3 * tH3 + 304. * sH2 * tH4 + 171. * sH * tH5
908  + 42. * tH6)
909  + pow3(s3 * s3) * 2. * stH * (42. * sH4 + 106. * sH3 * tH
910  + 119. * sH2 * tH2 + 106. * sH * tH3 + 42. * tH4)
911  - pow4(s3) * pow3(s3) * (35. * sH4 + 99. * sH3 * tH
912  + 120. * sH2 * tH2 + 99. * sH * tH3 + 35. * tH4)
913  + pow4(s3 * s3) * 7. * stH * ssttH)
914  / (sH * tH * uH * s3 * m3 * pow3(stH * tuH * usH));
915  }
916 
917  // Answer.
918  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
919 
920 }
921 
922 //--------------------------------------------------------------------------
923 
924 // Select identity, colour and anticolour.
925 
926 void Sigma2gg2QQbarX8g::setIdColAcol() {
927 
928  // Flavours are trivial.
929  setId( id1, id2, idHad, 21);
930 
931  // Split total contribution into different colour flows just like in
932  // g g -> g g (with kinematics recalculated for massless partons).
933  double sHr = - (tH + uH);
934  double sH2r = sHr * sHr;
935  double sigTS = tH2/sH2r + 2.*tH/sHr + 3. + 2.*sHr/tH + sH2r/tH2;
936  double sigUS = uH2/sH2r + 2.*uH/sHr + 3. + 2.*sHr/uH + sH2r/uH2;
937  double sigTU = tH2/uH2 + 2.*tH/uH + 3. + 2.*uH/tH + uH2/tH2;
938  double sigSum = sigTS + sigUS + sigTU;
939 
940  // Three colour flow topologies, each with two orientations.
941  double sigRand = sigSum * rndmPtr->flat();
942  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
943  else if (sigRand < sigTS + sigUS)
944  setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
945  else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
946  if (rndmPtr->flat() > 0.5) swapColAcol();
947 
948 }
949 
950 //==========================================================================
951 
952 // Sigma2qg2QQbarX8q class.
953 // Cross section q g -> QQbar[X(8)] q (Q = c or b, X = 3S1, 1S0 or 3PJ).
954 
955 //--------------------------------------------------------------------------
956 
957 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
958 
959 void Sigma2qg2QQbarX8q::sigmaKin() {
960 
961  // Calculate kinematics dependence.
962  double stH = sH + tH;
963  double tuH = tH + uH;
964  double usH = uH + sH;
965  double stH2 = stH * stH;
966  double tuH2 = tuH * tuH;
967  double usH2 = usH * usH;
968  double sig = 0.;
969  if (stateSave == 0) {
970  sig = - (M_PI / 27.)* (4. * (sH2 + uH2) - sH * uH) * (stH2 +tuH2)
971  / (s3 * m3 * sH * uH * usH2);
972  } else if (stateSave == 1) {
973  sig = - (5. * M_PI / 18.) * (sH2 + uH2) / (m3 * tH * usH2);
974  } else if (stateSave == 2) {
975  sig = - (10. * M_PI / 9.) * ( (7. * usH + 8. * tH) * (sH2 + uH2)
976  + 4. * tH * (2. * pow2(s3) - stH2 - tuH2) )
977  / (s3 * m3 * tH * usH2 * usH);
978  }
979 
980  // Answer.
981  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
982 
983 }
984 
985 //--------------------------------------------------------------------------
986 
987 // Select identity, colour and anticolour.
988 
989 void Sigma2qg2QQbarX8q::setIdColAcol() {
990 
991  // Flavours are trivial.
992  int idq = (id2 == 21) ? id1 : id2;
993  setId( id1, id2, idHad, idq);
994 
995  // tH defined between q_in and q_out: must swap tHat <-> uHat if q g in.
996  swapTU = (id2 == 21);
997 
998  // Split total contribution into different colour flows just like in
999  // q g -> q g (with kinematics recalculated for massless partons).
1000  double sHr = - (tH + uH);
1001  double sH2r = sHr * sHr;
1002  double sigTS = uH2/tH2 - (4./9.) * uH/sHr;
1003  double sigTU = sH2r/tH2 - (4./9.) * sHr/uH;
1004  double sigSum = sigTS + sigTU;
1005 
1006  // Two colour flow topologies. Swap if first is gluon, or when antiquark.
1007  double sigRand = sigSum * rndmPtr->flat();
1008  if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 2, 3, 3, 0);
1009  else setColAcol( 1, 0, 2, 3, 1, 3, 2, 0);
1010  if (id1 == 21) swapCol12();
1011  if (idq < 0) swapColAcol();
1012 
1013 }
1014 
1015 //==========================================================================
1016 
1017 // Sigma2qqbar2QQbarX8g class.
1018 // Cross section q qbar -> QQbar[X(8)] g (Q = c or b, X = 3S1, 1S0 or 3PJ).
1019 
1020 //--------------------------------------------------------------------------
1021 
1022 // Evaluate d(sigmaHat)/d(tHat); no explicit flavour dependence.
1023 
1024 void Sigma2qqbar2QQbarX8g::sigmaKin() {
1025 
1026  // Calculate kinematics dependence.
1027  double stH = sH + tH;
1028  double tuH = tH + uH;
1029  double usH = uH + sH;
1030  double stH2 = stH * stH;
1031  double tuH2 = tuH * tuH;
1032  double usH2 = usH * usH;
1033  double sig = 0.;
1034  if (stateSave == 0) {
1035  sig = (8. * M_PI / 81.) * (4. * (tH2 + uH2) - tH * uH)
1036  * (stH2 + usH2) / (s3 * m3 * tH * uH * tuH2);
1037  } else if (stateSave == 1) {
1038  sig = (20. * M_PI / 27.) * (tH2 + uH2) / (m3 * sH * tuH2);
1039  } else if (stateSave == 2) {
1040  sig = (80. * M_PI / 27.) * ( (7. * tuH + 8. * sH) * (tH2 + uH2)
1041  + 4. * sH * (2. * pow2(s3) - stH2 -usH2) )
1042  / (s3 * m3 * sH * tuH2 * tuH);
1043  }
1044 
1045  // Answer.
1046  sigma = (M_PI/sH2) * pow3(alpS) * oniumME * sig;
1047 
1048 }
1049 
1050 //--------------------------------------------------------------------------
1051 
1052 // Select identity, colour and anticolour.
1053 
1054 void Sigma2qqbar2QQbarX8g::setIdColAcol() {
1055 
1056  // Flavours are trivial.
1057  setId( id1, id2, idHad, 21);
1058 
1059  // Split total contribution into different colour flows just like in
1060  // q qbar -> g g (with kinematics recalculated for massless partons).
1061  double sHr = - (tH + uH);
1062  double sH2r = sHr * sHr;
1063  double sigTS = (4. / 9.) * uH / tH - uH2 / sH2r;
1064  double sigUS = (4. / 9.) * tH / uH - tH2 / sH2r;
1065  double sigSum = sigTS + sigUS;
1066 
1067  // Two colour flow topologies. Swap if first is antiquark.
1068  double sigRand = sigSum * rndmPtr->flat();
1069  if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
1070  else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
1071  if (id1 < 0) swapColAcol();
1072 
1073 }
1074 
1075 //==========================================================================
1076 
1077 } // end namespace Pythia8
1078 
Generated by   doxygen 1.8.5