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