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