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