Line data Source code
1 : // SigmaQCD.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 : // QCD simulation classes.
8 :
9 : #include "Pythia8/SigmaQCD.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // Sigma0AB2AB class.
16 : // Cross section for elastic scattering A B -> A B.
17 :
18 : //--------------------------------------------------------------------------
19 :
20 : // Select identity, colour and anticolour.
21 :
22 : void Sigma0AB2AB::setIdColAcol() {
23 :
24 : // Flavours and colours are trivial.
25 0 : setId( idA, idB, idA, idB);
26 0 : setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
27 0 : }
28 :
29 : //==========================================================================
30 :
31 : // Sigma0AB2XB class.
32 : // Cross section for single diffractive scattering A B -> X B.
33 :
34 : //--------------------------------------------------------------------------
35 :
36 : // Select identity, colour and anticolour.
37 :
38 : void Sigma0AB2XB::setIdColAcol() {
39 :
40 : // Flavours and colours are trivial.
41 0 : int idX = 10* (abs(idA) / 10) + 9900000;
42 0 : if (idA < 0) idX = -idX;
43 0 : setId( idA, idB, idX, idB);
44 0 : setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
45 :
46 0 : }
47 :
48 : //==========================================================================
49 :
50 : // Sigma0AB2AX class.
51 : // Cross section for single diffractive scattering A B -> A X.
52 :
53 : //--------------------------------------------------------------------------
54 :
55 : // Select identity, colour and anticolour.
56 :
57 : void Sigma0AB2AX::setIdColAcol() {
58 :
59 : // Flavours and colours are trivial.
60 0 : int idX = 10* (abs(idB) / 10) + 9900000;
61 0 : if (idB < 0) idX = -idX;
62 0 : setId( idA, idB, idA, idX);
63 0 : setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
64 :
65 0 : }
66 :
67 : //==========================================================================
68 :
69 : // Sigma0AB2XX class.
70 : // Cross section for double diffractive scattering A B -> X X.
71 :
72 : //--------------------------------------------------------------------------
73 :
74 : // Select identity, colour and anticolour.
75 :
76 : void Sigma0AB2XX::setIdColAcol() {
77 :
78 : // Flavours and colours are trivial.
79 0 : int idX1 = 10* (abs(idA) / 10) + 9900000;
80 0 : if (idA < 0) idX1 = -idX1;
81 0 : int idX2 = 10* (abs(idB) / 10) + 9900000;
82 0 : if (idB < 0) idX2 = -idX2;
83 0 : setId( idA, idB, idX1, idX2);
84 0 : setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
85 :
86 0 : }
87 :
88 : //==========================================================================
89 :
90 : // Sigma0AB2AXB class.
91 : // Cross section for central scattering A B -> A X B.
92 :
93 : //--------------------------------------------------------------------------
94 :
95 : // Select identity, colour and anticolour.
96 :
97 : void Sigma0AB2AXB::setIdColAcol() {
98 :
99 : // Central diffractive state represented by rho_diffr0. Colours trivial.
100 : int idX = 9900110;
101 0 : setId( idA, idB, idA, idB,idX);
102 0 : setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
103 :
104 0 : }
105 :
106 : //==========================================================================
107 :
108 : // Sigma2gg2gg class.
109 : // Cross section for g g -> g g.
110 :
111 : //--------------------------------------------------------------------------
112 :
113 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
114 :
115 : void Sigma2gg2gg::sigmaKin() {
116 :
117 : // Calculate kinematics dependence.
118 0 : sigTS = (9./4.) * (tH2 / sH2 + 2. * tH / sH + 3. + 2. * sH / tH
119 0 : + sH2 / tH2);
120 0 : sigUS = (9./4.) * (uH2 / sH2 + 2. * uH / sH + 3. + 2. * sH / uH
121 0 : + sH2 / uH2);
122 0 : sigTU = (9./4.) * (tH2 / uH2 + 2. * tH / uH + 3. + 2. * uH / tH
123 0 : + uH2 / tH2);
124 0 : sigSum = sigTS + sigUS + sigTU;
125 :
126 : // Answer contains factor 1/2 from identical gluons.
127 0 : sigma = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
128 :
129 0 : }
130 :
131 : //--------------------------------------------------------------------------
132 :
133 : // Select identity, colour and anticolour.
134 :
135 : void Sigma2gg2gg::setIdColAcol() {
136 :
137 : // Flavours are trivial.
138 0 : setId( id1, id2, 21, 21);
139 :
140 : // Three colour flow topologies, each with two orientations.
141 0 : double sigRand = sigSum * rndmPtr->flat();
142 0 : if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
143 0 : else if (sigRand < sigTS + sigUS)
144 0 : setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
145 0 : else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
146 0 : if (rndmPtr->flat() > 0.5) swapColAcol();
147 :
148 0 : }
149 :
150 : //==========================================================================
151 :
152 : // Sigma2gg2qqbar class.
153 : // Cross section for g g -> q qbar (q = u, d, s, i.e. almost massless).
154 :
155 : //--------------------------------------------------------------------------
156 :
157 : // Initialize process.
158 :
159 : void Sigma2gg2qqbar::initProc() {
160 :
161 : // Read number of quarks to be considered in massless approximation.
162 0 : nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
163 :
164 0 : }
165 :
166 : //--------------------------------------------------------------------------
167 :
168 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
169 :
170 : void Sigma2gg2qqbar::sigmaKin() {
171 :
172 : // Pick new flavour.
173 0 : idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
174 0 : mNew = particleDataPtr->m0(idNew);
175 0 : m2New = mNew*mNew;
176 :
177 : // Calculate kinematics dependence.
178 0 : sigTS = 0.;
179 0 : sigUS = 0.;
180 0 : if (sH > 4. * m2New) {
181 0 : sigTS = (1./6.) * uH / tH - (3./8.) * uH2 / sH2;
182 0 : sigUS = (1./6.) * tH / uH - (3./8.) * tH2 / sH2;
183 0 : }
184 0 : sigSum = sigTS + sigUS;
185 :
186 : // Answer is proportional to number of outgoing flavours.
187 0 : sigma = (M_PI / sH2) * pow2(alpS) * nQuarkNew * sigSum;
188 :
189 0 : }
190 :
191 : //--------------------------------------------------------------------------
192 :
193 : // Select identity, colour and anticolour.
194 :
195 : void Sigma2gg2qqbar::setIdColAcol() {
196 :
197 : // Flavours are trivial.
198 0 : setId( id1, id2, idNew, -idNew);
199 :
200 : // Two colour flow topologies.
201 0 : double sigRand = sigSum * rndmPtr->flat();
202 0 : if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
203 0 : else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
204 :
205 0 : }
206 :
207 : //==========================================================================
208 :
209 : // Sigma2qg2qg class.
210 : // Cross section for q g -> q g.
211 :
212 : //--------------------------------------------------------------------------
213 :
214 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
215 :
216 : void Sigma2qg2qg::sigmaKin() {
217 :
218 : // Calculate kinematics dependence.
219 0 : sigTS = uH2 / tH2 - (4./9.) * uH / sH;
220 0 : sigTU = sH2 / tH2 - (4./9.) * sH / uH;
221 0 : sigSum = sigTS + sigTU;
222 :
223 : // Answer.
224 0 : sigma = (M_PI / sH2) * pow2(alpS) * sigSum;
225 :
226 0 : }
227 :
228 : //--------------------------------------------------------------------------
229 :
230 : // Select identity, colour and anticolour.
231 :
232 : void Sigma2qg2qg::setIdColAcol() {
233 :
234 : // Outgoing = incoming flavours.
235 0 : setId( id1, id2, id1, id2);
236 :
237 : // Two colour flow topologies. Swap if first is gluon, or when antiquark.
238 0 : double sigRand = sigSum * rndmPtr->flat();
239 0 : if (sigRand < sigTS) setColAcol( 1, 0, 2, 1, 3, 0, 2, 3);
240 0 : else setColAcol( 1, 0, 2, 3, 2, 0, 1, 3);
241 0 : if (id1 == 21) swapCol1234();
242 0 : if (id1 < 0 || id2 < 0) swapColAcol();
243 :
244 0 : }
245 :
246 : //==========================================================================
247 :
248 : // Sigma2qq2qq class.
249 : // Cross section for q qbar' -> q qbar' or q q' -> q q'
250 : // (qbar qbar' -> qbar qbar'), q' may be same as q.
251 :
252 : //--------------------------------------------------------------------------
253 :
254 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
255 :
256 : void Sigma2qq2qq::sigmaKin() {
257 :
258 : // Calculate kinematics dependence for different terms.
259 0 : sigT = (4./9.) * (sH2 + uH2) / tH2;
260 0 : sigU = (4./9.) * (sH2 + tH2) / uH2;
261 0 : sigTU = - (8./27.) * sH2 / (tH * uH);
262 0 : sigST = - (8./27.) * uH2 / (sH * tH);
263 :
264 0 : }
265 :
266 : //--------------------------------------------------------------------------
267 :
268 :
269 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
270 :
271 : double Sigma2qq2qq::sigmaHat() {
272 :
273 : // Combine cross section terms; factor 1/2 when identical quarks.
274 0 : if (id2 == id1) sigSum = 0.5 * (sigT + sigU + sigTU);
275 0 : else if (id2 == -id1) sigSum = sigT + sigST;
276 0 : else sigSum = sigT;
277 :
278 : // Answer.
279 0 : return (M_PI/sH2) * pow2(alpS) * sigSum;
280 :
281 : }
282 :
283 : //--------------------------------------------------------------------------
284 :
285 : // Select identity, colour and anticolour.
286 :
287 : void Sigma2qq2qq::setIdColAcol() {
288 :
289 : // Outgoing = incoming flavours.
290 0 : setId( id1, id2, id1, id2);
291 :
292 : // Colour flow topologies. Swap when antiquarks.
293 0 : if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
294 0 : else setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
295 0 : if (id2 == id1 && (sigT + sigU) * rndmPtr->flat() > sigT)
296 0 : setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
297 0 : if (id1 < 0) swapColAcol();
298 :
299 0 : }
300 :
301 : //==========================================================================
302 :
303 : // Sigma2qqbar2gg class.
304 : // Cross section for q qbar -> g g.
305 :
306 : //--------------------------------------------------------------------------
307 :
308 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
309 :
310 : void Sigma2qqbar2gg::sigmaKin() {
311 :
312 : // Calculate kinematics dependence.
313 0 : sigTS = (32./27.) * uH / tH - (8./3.) * uH2 / sH2;
314 0 : sigUS = (32./27.) * tH / uH - (8./3.) * tH2 / sH2;
315 0 : sigSum = sigTS + sigUS;
316 :
317 : // Answer contains factor 1/2 from identical gluons.
318 0 : sigma = (M_PI / sH2) * pow2(alpS) * 0.5 * sigSum;
319 :
320 0 : }
321 :
322 : //--------------------------------------------------------------------------
323 :
324 : // Select identity, colour and anticolour.
325 :
326 : void Sigma2qqbar2gg::setIdColAcol() {
327 :
328 : // Outgoing flavours trivial.
329 0 : setId( id1, id2, 21, 21);
330 :
331 : // Two colour flow topologies. Swap if first is antiquark.
332 0 : double sigRand = sigSum * rndmPtr->flat();
333 0 : if (sigRand < sigTS) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
334 0 : else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
335 0 : if (id1 < 0) swapColAcol();
336 :
337 0 : }
338 :
339 : //==========================================================================
340 :
341 : // Sigma2qqbar2qqbarNew class.
342 : // Cross section q qbar -> q' qbar'.
343 :
344 : //--------------------------------------------------------------------------
345 :
346 : // Initialize process.
347 :
348 : void Sigma2qqbar2qqbarNew::initProc() {
349 :
350 : // Read number of quarks to be considered in massless approximation.
351 0 : nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
352 :
353 0 : }
354 :
355 : //--------------------------------------------------------------------------
356 :
357 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
358 :
359 : void Sigma2qqbar2qqbarNew::sigmaKin() {
360 :
361 : // Pick new flavour.
362 0 : idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
363 0 : mNew = particleDataPtr->m0(idNew);
364 0 : m2New = mNew*mNew;
365 :
366 : // Calculate kinematics dependence.
367 0 : sigS = 0.;
368 0 : if (sH > 4. * m2New) sigS = (4./9.) * (tH2 + uH2) / sH2;
369 :
370 : // Answer is proportional to number of outgoing flavours.
371 0 : sigma = (M_PI / sH2) * pow2(alpS) * nQuarkNew * sigS;
372 :
373 0 : }
374 :
375 : //--------------------------------------------------------------------------
376 :
377 : // Select identity, colour and anticolour.
378 :
379 : void Sigma2qqbar2qqbarNew::setIdColAcol() {
380 :
381 : // Set outgoing flavours ones.
382 0 : id3 = (id1 > 0) ? idNew : -idNew;
383 0 : setId( id1, id2, id3, -id3);
384 :
385 : // Colour flow topologies. Swap when antiquarks.
386 0 : setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
387 0 : if (id1 < 0) swapColAcol();
388 :
389 0 : }
390 :
391 : //==========================================================================
392 :
393 : // Sigma2gg2QQbar class.
394 : // Cross section g g -> Q Qbar (Q = c, b or t).
395 : // Only provided for fixed m3 = m4 so do some gymnastics:
396 : // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
397 : // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
398 : // but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
399 : // tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
400 :
401 : //--------------------------------------------------------------------------
402 :
403 : // Initialize process.
404 :
405 : void Sigma2gg2QQbar::initProc() {
406 :
407 : // Process name.
408 0 : nameSave = "g g -> Q Qbar";
409 0 : if (idNew == 4) nameSave = "g g -> c cbar";
410 0 : if (idNew == 5) nameSave = "g g -> b bbar";
411 0 : if (idNew == 6) nameSave = "g g -> t tbar";
412 0 : if (idNew == 7) nameSave = "g g -> b' b'bar";
413 0 : if (idNew == 8) nameSave = "g g -> t' t'bar";
414 :
415 : // Secondary open width fraction.
416 0 : openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
417 :
418 0 : }
419 :
420 : //--------------------------------------------------------------------------
421 :
422 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
423 :
424 : void Sigma2gg2QQbar::sigmaKin() {
425 :
426 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
427 0 : double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
428 0 : double tHQ = -0.5 * (sH - tH + uH);
429 0 : double uHQ = -0.5 * (sH + tH - uH);
430 0 : double tHQ2 = tHQ * tHQ;
431 0 : double uHQ2 = uHQ * uHQ;
432 :
433 : // Calculate kinematics dependence.
434 0 : double tumHQ = tHQ * uHQ - s34Avg * sH;
435 0 : sigTS = ( uHQ / tHQ - 2.25 * uHQ2 / sH2 + 4.5 * s34Avg * tumHQ
436 0 : / ( sH * tHQ2) + 0.5 * s34Avg * (tHQ + s34Avg) / tHQ2
437 0 : - s34Avg*s34Avg / (sH * tHQ) ) / 6.;
438 0 : sigUS = ( tHQ / uHQ - 2.25 * tHQ2 / sH2 + 4.5 * s34Avg * tumHQ
439 0 : / ( sH * uHQ2) + 0.5 * s34Avg * (uHQ + s34Avg) / uHQ2
440 0 : - s34Avg*s34Avg / (sH * uHQ) ) / 6.;
441 0 : sigSum = sigTS + sigUS;
442 :
443 : // Answer.
444 0 : sigma = (M_PI / sH2) * pow2(alpS) * sigSum * openFracPair;
445 :
446 0 : }
447 :
448 : //--------------------------------------------------------------------------
449 :
450 : // Select identity, colour and anticolour.
451 :
452 : void Sigma2gg2QQbar::setIdColAcol() {
453 :
454 : // Flavours are trivial.
455 0 : setId( id1, id2, idNew, -idNew);
456 :
457 : // Two colour flow topologies.
458 0 : double sigRand = sigSum * rndmPtr->flat();
459 0 : if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
460 0 : else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
461 :
462 0 : }
463 :
464 : //--------------------------------------------------------------------------
465 :
466 : // Evaluate weight for decay angles of W in top decay.
467 :
468 : double Sigma2gg2QQbar::weightDecay( Event& process, int iResBeg,
469 : int iResEnd) {
470 :
471 : // For top decay hand over to standard routine, else done.
472 0 : if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
473 0 : return weightTopDecay( process, iResBeg, iResEnd);
474 0 : else return 1.;
475 :
476 0 : }
477 :
478 : //==========================================================================
479 :
480 : // Sigma2qqbar2QQbar class.
481 : // Cross section q qbar -> Q Qbar (Q = c, b or t).
482 : // Only provided for fixed m3 = m4 so do some gymnastics:
483 : // i) s34Avg picked so that beta34 same when s3, s4 -> s34Avg.
484 : // ii) tHQ = tH - mQ^2 = -0.5 sH (1 - beta34 cos(thetaH)) for m3 = m4 = mQ,
485 : // but tH - uH = sH beta34 cos(thetaH) also for m3 != m4, so use
486 : // tH, uH selected for m3 != m4 to derive tHQ, uHQ valid for m3 = m4.
487 :
488 : //--------------------------------------------------------------------------
489 :
490 : // Initialize process, especially parton-flux object.
491 :
492 : void Sigma2qqbar2QQbar::initProc() {
493 :
494 : // Process name.
495 0 : nameSave = "q qbar -> Q Qbar";
496 0 : if (idNew == 4) nameSave = "q qbar -> c cbar";
497 0 : if (idNew == 5) nameSave = "q qbar -> b bbar";
498 0 : if (idNew == 6) nameSave = "q qbar -> t tbar";
499 0 : if (idNew == 7) nameSave = "q qbar -> b' b'bar";
500 0 : if (idNew == 8) nameSave = "q qbar -> t' t'bar";
501 :
502 : // Secondary open width fraction.
503 0 : openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
504 :
505 0 : }
506 :
507 : //--------------------------------------------------------------------------
508 :
509 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
510 :
511 : void Sigma2qqbar2QQbar::sigmaKin() {
512 :
513 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
514 0 : double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
515 0 : double tHQ = -0.5 * (sH - tH + uH);
516 0 : double uHQ = -0.5 * (sH + tH - uH);
517 0 : double tHQ2 = tHQ * tHQ;
518 0 : double uHQ2 = uHQ * uHQ;
519 :
520 : // Calculate kinematics dependence.
521 0 : double sigS = (4./9.) * ((tHQ2 + uHQ2) / sH2 + 2. * s34Avg / sH);
522 :
523 : // Answer.
524 0 : sigma = (M_PI / sH2) * pow2(alpS) * sigS * openFracPair;
525 :
526 0 : }
527 :
528 : //--------------------------------------------------------------------------
529 :
530 : // Select identity, colour and anticolour.
531 :
532 : void Sigma2qqbar2QQbar::setIdColAcol() {
533 :
534 : // Set outgoing flavours.
535 0 : id3 = (id1 > 0) ? idNew : -idNew;
536 0 : setId( id1, id2, id3, -id3);
537 :
538 : // Colour flow topologies. Swap when antiquarks.
539 0 : setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
540 0 : if (id1 < 0) swapColAcol();
541 :
542 0 : }
543 :
544 : //--------------------------------------------------------------------------
545 :
546 : // Evaluate weight for decay angles of W in top decay.
547 :
548 : double Sigma2qqbar2QQbar::weightDecay( Event& process, int iResBeg,
549 : int iResEnd) {
550 :
551 : // For top decay hand over to standard routine, else done.
552 0 : if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
553 0 : return weightTopDecay( process, iResBeg, iResEnd);
554 0 : else return 1.;
555 :
556 0 : }
557 :
558 :
559 : //==========================================================================
560 :
561 : // Sigma3gg2ggg class.
562 : // Cross section for g g -> g g g.
563 :
564 : //--------------------------------------------------------------------------
565 :
566 : // Evaluate |M|^2 - no incoming flavour dependence.
567 :
568 : void Sigma3gg2ggg::sigmaKin() {
569 :
570 : // Calculate all four-vector products.
571 0 : Vec4 p1cm( 0., 0., 0.5 * mH, 0.5 * mH);
572 0 : Vec4 p2cm( 0., 0., -0.5 * mH, 0.5 * mH);
573 0 : pp[1][2] = p1cm * p2cm;
574 0 : pp[1][3] = p1cm * p3cm;
575 0 : pp[1][4] = p1cm * p4cm;
576 0 : pp[1][5] = p1cm * p5cm;
577 0 : pp[2][3] = p2cm * p3cm;
578 0 : pp[2][4] = p2cm * p4cm;
579 0 : pp[2][5] = p2cm * p5cm;
580 0 : pp[3][4] = p3cm * p4cm;
581 0 : pp[3][5] = p3cm * p5cm;
582 0 : pp[4][5] = p4cm * p5cm;
583 0 : for (int i = 1; i < 5; ++i)
584 0 : for (int j = i + 1; j < 6; ++j) pp[j][i] = pp[i][j];
585 :
586 : // Cross section, in three main sections.
587 0 : double num1 = cycle(1,2,3,4,5) + cycle(1,2,3,5,4) + cycle(1,2,4,3,5)
588 0 : + cycle(1,2,4,5,3) + cycle(1,2,5,3,4) + cycle(1,2,5,4,3)
589 0 : + cycle(1,3,2,4,5) + cycle(1,3,2,5,4) + cycle(1,3,4,2,5)
590 0 : + cycle(1,3,5,2,4) + cycle(1,4,2,3,5) + cycle(1,4,3,2,5);
591 0 : double num2 = pow4(pp[1][2]) + pow4(pp[1][3]) + pow4(pp[1][4])
592 0 : + pow4(pp[1][5]) + pow4(pp[2][3]) + pow4(pp[2][4])
593 0 : + pow4(pp[2][5]) + pow4(pp[3][4]) + pow4(pp[3][5])
594 0 : + pow4(pp[4][5]);
595 0 : double den = pp[1][2] * pp[1][3] * pp[1][4] * pp[1][5] * pp[2][3]
596 0 : * pp[2][4] * pp[2][5] * pp[3][4] * pp[3][5] * pp[4][5];
597 :
598 : // Answer has a factor 6 due to identical gluons
599 : // This is cancelled by phase space factor (1 / 6)
600 0 : sigma = pow3(4. * M_PI * alpS) * (27./16.) * num1 * num2 / den;
601 :
602 0 : }
603 :
604 : //--------------------------------------------------------------------------
605 :
606 : // Select identity, colour and anticolour.
607 :
608 : void Sigma3gg2ggg::setIdColAcol() {
609 :
610 : // Flavours are trivial.
611 0 : setId( id1, id2, 21, 21, 21);
612 :
613 : // Three colour flow topologies, each with two orientations.
614 : /*
615 : double sigRand = sigSum * rndmPtr->flat();
616 : if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
617 : else if (sigRand < sigTS + sigUS)
618 : setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
619 : else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
620 : if (rndmPtr->flat() > 0.5) swapColAcol();
621 : */
622 :
623 : // Temporary solution.
624 0 : setColAcol( 1, 2, 2, 3, 1, 4, 4, 5, 5, 3);
625 0 : }
626 :
627 :
628 : //==========================================================================
629 :
630 : // Sigma3qqbar2ggg class.
631 : // Cross section for q qbar -> g g g.
632 :
633 : //--------------------------------------------------------------------------
634 :
635 : // Evaluate |M|^2 - no incoming flavour dependence.
636 : void Sigma3qqbar2ggg::sigmaKin() {
637 :
638 : // Setup four-vectors
639 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
640 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
641 0 : pCM[2] = p3cm;
642 0 : pCM[3] = p4cm;
643 0 : pCM[4] = p5cm;
644 :
645 : // Calculate |M|^2
646 : // Answer has a factor 6 due to identical gluons,
647 : // which is cancelled by phase space factor (1 / 6)
648 0 : sigma = m2Calc();
649 :
650 0 : }
651 :
652 : //--------------------------------------------------------------------------
653 :
654 : // |M|^2
655 :
656 : inline double Sigma3qqbar2ggg::m2Calc() {
657 :
658 : // Calculate four-products
659 0 : double sHnow = (pCM[0] + pCM[1]).m2Calc();
660 0 : double sHhalf = sH / 2.;
661 :
662 : // qbar (p+) + q(p-) -> g(k1) g(k2) g(k3)
663 : // a_i = (p+ . ki), i = 1, 2, 3
664 : // b_i = (p- . ki), i = 1, 2, 3
665 0 : a[0] = pCM[0] * pCM[2];
666 0 : a[1] = pCM[0] * pCM[3];
667 0 : a[2] = pCM[0] * pCM[4];
668 0 : b[0] = pCM[1] * pCM[2];
669 0 : b[1] = pCM[1] * pCM[3];
670 0 : b[2] = pCM[1] * pCM[4];
671 :
672 0 : pp[0][1] = pCM[2] * pCM[3];
673 0 : pp[1][2] = pCM[3] * pCM[4];
674 0 : pp[2][0] = pCM[4] * pCM[2];
675 :
676 : // ab[i][j] = a_i * b_j + a_j * b_i
677 0 : ab[0][1] = a[0] * b[1] + a[1] * b[0];
678 0 : ab[1][2] = a[1] * b[2] + a[2] * b[1];
679 0 : ab[2][0] = a[2] * b[0] + a[0] * b[2];
680 :
681 : // Cross section
682 0 : double num1 = a[0] * b[0] * (a[0] * a[0] + b[0] * b[0]) +
683 0 : a[1] * b[1] * (a[1] * a[1] + b[1] * b[1]) +
684 0 : a[2] * b[2] * (a[2] * a[2] + b[2] * b[2]);
685 0 : double den1 = a[0] * a[1] * a[2] * b[0] * b[1] * b[2];
686 0 : double num2 = - ( ab[0][1] / pp[0][1] )
687 0 : - ( ab[1][2] / pp[1][2] )
688 0 : - ( ab[2][0] / pp[2][0] );
689 0 : double num3 = a[2] * b[2] * ab[0][1] / (pp[1][2] * pp[2][0] ) +
690 0 : a[0] * b[0] * ab[1][2] / (pp[2][0] * pp[0][1] ) +
691 0 : a[1] * b[1] * ab[2][0] / (pp[0][1] * pp[1][2] );
692 :
693 : // Final answer
694 0 : return pow3(4. * M_PI * alpS) * (8. / 324.) * (num1 / den1) *
695 0 : ( sHhalf + 9. * (sHhalf + num2) + (2. * 81. / sHnow) * num3 );
696 :
697 : }
698 :
699 : //--------------------------------------------------------------------------
700 :
701 : // Select identity, colour and anticolour.
702 :
703 : void Sigma3qqbar2ggg::setIdColAcol(){
704 :
705 : // Flavours are trivial.
706 0 : setId( id1, id2, 21, 21, 21);
707 :
708 : // Temporary solution.
709 0 : setColAcol( 1, 0, 0, 2, 1, 3, 3, 4, 4, 2);
710 0 : if (id1 < 0) swapColAcol();
711 0 : }
712 :
713 : //--------------------------------------------------------------------------
714 :
715 : // Map a final state configuration
716 :
717 : inline void Sigma3qqbar2ggg::mapFinal() {
718 0 : switch (config) {
719 0 : case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
720 0 : case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
721 0 : case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
722 0 : case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
723 0 : case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
724 0 : case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
725 : }
726 0 : }
727 :
728 : //==========================================================================
729 :
730 : // Sigma3qg2qgg class.
731 : // Cross section for q g -> q g g.
732 : // Crossed relation from q qbar -> g g g:
733 : // qbar(p+) q(p-) -> g(k1) g(k2) g(k3)
734 :
735 : //--------------------------------------------------------------------------
736 :
737 : // Evaluate |M|^2 - no incoming flavour dependence
738 : // Note: two different contributions from gq and qg incoming
739 :
740 : void Sigma3qg2qgg::sigmaKin() {
741 :
742 : // Pick a final state configuration
743 0 : pickFinal();
744 :
745 : // gq and qg incoming
746 0 : for (int i = 0; i < 2; i++) {
747 :
748 : // Map incoming four-vectors to p+, p-, k1, k2, k3
749 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
750 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
751 0 : mapFinal();
752 :
753 : // Crossing
754 0 : swap(pCM[i], pCM[2]);
755 :
756 : // |M|^2
757 : // XXX - Extra factor of (3) from selecting a final state
758 : // configuration (already a factor of 2 in the original
759 : // answer due to two identical final state gluons)???
760 : // Extra factor of (3 / 8) as average over incoming gluon
761 0 : sigma[i] = 3. * (3. / 8.) * m2Calc();
762 :
763 : } // for (int i = 0; i < 2; i++)
764 :
765 0 : }
766 :
767 : //--------------------------------------------------------------------------
768 :
769 : // Evaluate |M|^2 - incoming flavour dependence
770 : // Pick from two configurations calculated previously
771 :
772 : double Sigma3qg2qgg::sigmaHat() {
773 : // gq or qg incoming
774 0 : return (id1 == 21) ? sigma[0] : sigma[1];
775 : }
776 :
777 : //--------------------------------------------------------------------------
778 :
779 : // Select identity, colour and anticolour.
780 :
781 : void Sigma3qg2qgg::setIdColAcol(){
782 : // Outgoing flavours; only need to know where the quark is
783 0 : int qIdx = config / 2;
784 0 : int idTmp[3] = { 21, 21, 21 };
785 0 : idTmp[qIdx] = (id1 == 21) ? id2 : id1;
786 0 : setId( id1, id2, idTmp[0], idTmp[1], idTmp[2]);
787 :
788 : // Temporary solution
789 0 : if (qIdx == 0) setColAcol(1, 0, 2, 1, 4, 0, 3, 4, 2, 3);
790 0 : else if (qIdx == 1) setColAcol(1, 0, 2, 1, 3, 4, 4, 0, 2, 3);
791 0 : else setColAcol(1, 0, 2, 1, 3, 4, 2, 3, 4, 0);
792 : // gq or qg incoming
793 0 : if (id1 == 21) {
794 0 : swap( colSave[1], colSave[2]);
795 0 : swap(acolSave[1], acolSave[2]);
796 0 : }
797 : // qbar rather than q incoming
798 0 : if (id1 < 0 || id2 < 0) swapColAcol();
799 :
800 0 : }
801 :
802 : //==========================================================================
803 :
804 : // Sigma3gg2qqbarg class.
805 : // Cross section for g g -> q qbar g
806 : // Crossed relation from q qbar -> g g g
807 :
808 : //--------------------------------------------------------------------------
809 :
810 : // Initialize process.
811 :
812 : void Sigma3gg2qqbarg::initProc() {
813 :
814 : // Read number of quarks to be considered in massless approximation.
815 0 : nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
816 :
817 0 : }
818 :
819 : //--------------------------------------------------------------------------
820 :
821 : // Evaluate |M|^2 - no incoming flavour dependence.
822 :
823 : void Sigma3gg2qqbarg::sigmaKin() {
824 :
825 : // Incoming four-vectors
826 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
827 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
828 :
829 : // Pick and map a final state configuration
830 0 : pickFinal();
831 0 : mapFinal();
832 :
833 : // Crossing
834 0 : swap(pCM[0], pCM[2]);
835 0 : swap(pCM[1], pCM[3]);
836 :
837 : // |M|^2
838 : // Extra factor of (6.) from picking a final state configuration
839 : // Extra factor of nQuarkNew
840 : // Extra factor of (3. / 8.) ^ 2 as averaging over two incoming gluons
841 0 : sigma = 6. * nQuarkNew * (3. / 8.) * (3. / 8.) * m2Calc();
842 :
843 0 : }
844 :
845 : //--------------------------------------------------------------------------
846 :
847 : // Select identity, colour and anticolour.
848 :
849 : void Sigma3gg2qqbarg::setIdColAcol(){
850 :
851 : // Pick new flavour
852 0 : int idNew = 1 + int( nQuarkNew * rndmPtr->flat() );
853 :
854 : // Outgoing flavours; easiest just to map by hand
855 0 : switch (config) {
856 0 : case 0: id3 = idNew; id4 = -idNew; id5 = 21; break;
857 0 : case 1: id3 = idNew; id4 = 21; id5 = -idNew; break;
858 0 : case 2: id3 = -idNew; id4 = idNew; id5 = 21; break;
859 0 : case 3: id3 = 21; id4 = idNew; id5 = -idNew; break;
860 0 : case 4: id3 = -idNew; id4 = 21; id5 = idNew; break;
861 0 : case 5: id3 = 21; id4 = -idNew; id5 = idNew; break;
862 : }
863 0 : setId(id1, id2, id3, id4, id5);
864 :
865 : // Temporary solution
866 0 : switch (config) {
867 0 : case 0: setColAcol( 1, 2, 2, 3, 4, 0, 0, 3, 1, 4 ); break;
868 0 : case 1: setColAcol( 1, 2, 2, 3, 4, 0, 1, 4, 0, 3 ); break;
869 0 : case 2: setColAcol( 1, 2, 2, 3, 0, 3, 4, 0, 1, 4 ); break;
870 0 : case 3: setColAcol( 1, 2, 2, 3, 1, 4, 4, 0, 0, 3 ); break;
871 0 : case 4: setColAcol( 1, 2, 2, 3, 0, 3, 1, 4, 4, 0 ); break;
872 0 : case 5: setColAcol( 1, 2, 2, 3, 1, 4, 0, 3, 4, 0 ); break;
873 : }
874 :
875 0 : }
876 :
877 : //==========================================================================
878 :
879 : // Sigma3qq2qqgDiff class.
880 : // Cross section for q q' -> q q' g, q != q'
881 :
882 : //--------------------------------------------------------------------------
883 :
884 : // Evaluate |M|^2 - no incoming flavour dependence
885 :
886 : void Sigma3qq2qqgDiff::sigmaKin() {
887 :
888 : // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
889 :
890 : // Incoming four-vectors
891 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
892 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
893 : // Pick and map a final state configuration
894 0 : pickFinal();
895 0 : mapFinal();
896 :
897 : // |M|^2
898 : // Extra factor of (6.) from picking a final state configuration
899 0 : sigma = 6. * m2Calc();
900 0 : }
901 :
902 : //--------------------------------------------------------------------------
903 :
904 : // |M|^2
905 :
906 : inline double Sigma3qq2qqgDiff::m2Calc() {
907 :
908 : // Four-products
909 0 : s = (pCM[0] + pCM[1]).m2Calc();
910 0 : t = (pCM[0] - pCM[2]).m2Calc();
911 0 : u = (pCM[0] - pCM[3]).m2Calc();
912 0 : up = (pCM[1] - pCM[2]).m2Calc();
913 0 : sp = (pCM[2] + pCM[3]).m2Calc();
914 0 : tp = (pCM[1] - pCM[3]).m2Calc();
915 :
916 : // |M|^2
917 0 : double num1 = (s * s + sp * sp + u * u + up * up) / (t * tp);
918 0 : double den1 = (pCM[0] * pCM[4]) * (pCM[1] * pCM[4]) *
919 0 : (pCM[2] * pCM[4]) * (pCM[3] * pCM[4]);
920 0 : double num2 = (u + up) * (s * sp + t * tp - u * up) +
921 0 : u * (s * t + sp * tp) + up * (s * tp + sp * t);
922 0 : double num3 = (s + sp) * (s * sp - t * tp - u * up) +
923 0 : 2. * t * tp * (u + up) + 2. * u * up * (t + tp);
924 :
925 : // (N^2 - 1)^2 / 4N^3 = 16. / 27.
926 : // (N^2 - 1) / 4N^3 = 2. / 27.
927 0 : return (1. / 8.) * pow3(4. * M_PI * alpS) * num1 / den1 *
928 0 : ( (16. / 27.) * num2 - (2. / 27.) * num3 );
929 :
930 : }
931 :
932 : //--------------------------------------------------------------------------
933 :
934 : // Evaluate |M|^2 - incoming flavour dependence
935 :
936 : double Sigma3qq2qqgDiff::sigmaHat() {
937 : // Different incoming flavours only
938 0 : if (abs(id1) == abs(id2)) return 0.;
939 0 : return sigma;
940 0 : }
941 :
942 : //--------------------------------------------------------------------------
943 :
944 : // Select identity, colour and anticolour.
945 :
946 : void Sigma3qq2qqgDiff::setIdColAcol(){
947 :
948 : // Outgoing flavours; easiest just to map by hand
949 0 : switch (config) {
950 0 : case 0: id3 = id1; id4 = id2; id5 = 21; break;
951 0 : case 1: id3 = id1; id4 = 21; id5 = id2; break;
952 0 : case 2: id3 = id2; id4 = id1; id5 = 21; break;
953 0 : case 3: id3 = 21; id4 = id1; id5 = id2; break;
954 0 : case 4: id3 = id2; id4 = 21; id5 = id1; break;
955 0 : case 5: id3 = 21; id4 = id2; id5 = id1; break;
956 : }
957 0 : setId(id1, id2, id3, id4, id5);
958 :
959 : // Temporary solution; id1 and id2 can be q/qbar independently
960 0 : int cols[5][2];
961 0 : if (id1 > 0) {
962 0 : cols[0][0] = 1; cols[0][1] = 0;
963 0 : cols[2][0] = 1; cols[2][1] = 0;
964 0 : } else {
965 0 : cols[0][0] = 0; cols[0][1] = 1;
966 0 : cols[2][0] = 0; cols[2][1] = 1;
967 : }
968 0 : if (id2 > 0) {
969 0 : cols[1][0] = 2; cols[1][1] = 0;
970 0 : cols[3][0] = 3; cols[3][1] = 0;
971 0 : cols[4][0] = 2; cols[4][1] = 3;
972 0 : } else {
973 0 : cols[1][0] = 0; cols[1][1] = 2;
974 0 : cols[3][0] = 0; cols[3][1] = 3;
975 0 : cols[4][0] = 3; cols[4][1] = 2;
976 : }
977 : // Map correct final state configuration
978 : int i3 = 0, i4 = 0, i5 = 0;
979 0 : switch (config) {
980 0 : case 0: i3 = 2; i4 = 3; i5 = 4; break;
981 0 : case 1: i3 = 2; i4 = 4; i5 = 3; break;
982 0 : case 2: i3 = 3; i4 = 2; i5 = 4; break;
983 0 : case 3: i3 = 4; i4 = 2; i5 = 3; break;
984 0 : case 4: i3 = 3; i4 = 4; i5 = 2; break;
985 0 : case 5: i3 = 4; i4 = 3; i5 = 2; break;
986 : }
987 : // Put colours in place
988 0 : setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
989 0 : cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
990 0 : cols[i5][0], cols[i5][1]);
991 :
992 0 : }
993 :
994 : //--------------------------------------------------------------------------
995 :
996 : // Map a final state configuration
997 :
998 : inline void Sigma3qq2qqgDiff::mapFinal() {
999 0 : switch (config) {
1000 0 : case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
1001 0 : case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
1002 0 : case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
1003 0 : case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
1004 0 : case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
1005 0 : case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
1006 : }
1007 0 : }
1008 :
1009 :
1010 : //==========================================================================
1011 :
1012 : // Sigma3qqbar2qqbargDiff
1013 : // Cross section for q qbar -> q' qbar' g
1014 : // Crossed relation from q q' -> q q' g, q != q'
1015 :
1016 : //--------------------------------------------------------------------------
1017 :
1018 : // Initialize process.
1019 :
1020 : void Sigma3qqbar2qqbargDiff::initProc() {
1021 :
1022 : // Read number of quarks to be considered in massless approximation.
1023 0 : nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
1024 :
1025 0 : }
1026 :
1027 : //--------------------------------------------------------------------------
1028 :
1029 : // Evaluate |M|^2 - no incoming flavour dependence.
1030 :
1031 : void Sigma3qqbar2qqbargDiff::sigmaKin() {
1032 : // Overall 6 possibilities for final state ordering
1033 : // To keep symmetry between final states, always map to:
1034 : // 1) q1(p+) qbar1(p-) -> qbar2(q+) q2(q-) g(k)
1035 : // 2) qbar1(p+) q1(p-) -> q2(q+) qbar2(q-) g(k)
1036 : // Crossing p- and q+ gives:
1037 : // 1) q1(p+) q2(-q+) -> q1(-p-) q2(q-) g(k)
1038 : // 2) qbar1(p+) qbar2(-q+) -> qbar1(-p-) qbar2(q-) g(k)
1039 :
1040 : // Incoming four-vectors
1041 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1042 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1043 : // Pick and map a final state configuration
1044 0 : pickFinal();
1045 0 : mapFinal();
1046 :
1047 : // Crossing
1048 0 : swap(pCM[1], pCM[2]);
1049 0 : pCM[1] = -pCM[1];
1050 0 : pCM[2] = -pCM[2];
1051 :
1052 : // |M|^2
1053 : // Extra factor of (6.) from picking a final state configuration
1054 : // Extra factor of (nQuarkNew - 1) from new q/qbar pairs
1055 : // XXX - Extra factor of (2.) from second possible crossing???
1056 0 : sigma = 6. * (nQuarkNew - 1) * 2. * m2Calc();
1057 :
1058 0 : }
1059 :
1060 : //--------------------------------------------------------------------------
1061 :
1062 : // Select identity, colour and anticolour.
1063 :
1064 : void Sigma3qqbar2qqbargDiff::setIdColAcol(){
1065 :
1066 : // Pick new q qbar flavour with incoming flavour disallowed
1067 0 : int idNew = 1 + int( (nQuarkNew - 1) * rndmPtr->flat() );
1068 0 : if (idNew >= abs(id1)) ++idNew;
1069 : // For qbar q incoming, q+ is always mapped to q2
1070 : // For q qbar incoming, q+ is always mapped to qbar2
1071 0 : if (id1 > 0) idNew = -idNew;
1072 :
1073 : // Outgoing flavours; easiest just to map by hand
1074 0 : switch (config) {
1075 0 : case 0: id3 = idNew; id4 = -idNew; id5 = 21; break;
1076 0 : case 1: id3 = idNew; id4 = 21; id5 = -idNew; break;
1077 0 : case 2: id3 = -idNew; id4 = idNew; id5 = 21; break;
1078 0 : case 3: id3 = 21; id4 = idNew; id5 = -idNew; break;
1079 0 : case 4: id3 = -idNew; id4 = 21; id5 = idNew; break;
1080 0 : case 5: id3 = 21; id4 = -idNew; id5 = idNew; break;
1081 : }
1082 0 : setId(id1, id2, id3, id4, id5);
1083 :
1084 : // Temporary solution; start with q qbar -> qbar q g
1085 0 : int cols[5][2];
1086 0 : cols[0][0] = 1; cols[0][1] = 0;
1087 0 : cols[1][0] = 0; cols[1][1] = 2;
1088 0 : cols[2][0] = 0; cols[2][1] = 3;
1089 0 : cols[3][0] = 1; cols[3][1] = 0;
1090 0 : cols[4][0] = 3; cols[4][1] = 2;
1091 : // Map into correct place
1092 : int i3 = 0, i4 = 0, i5 = 0;
1093 0 : switch (config) {
1094 0 : case 0: i3 = 2; i4 = 3; i5 = 4; break;
1095 0 : case 1: i3 = 2; i4 = 4; i5 = 3; break;
1096 0 : case 2: i3 = 3; i4 = 2; i5 = 4; break;
1097 0 : case 3: i3 = 4; i4 = 2; i5 = 3; break;
1098 0 : case 4: i3 = 3; i4 = 4; i5 = 2; break;
1099 0 : case 5: i3 = 4; i4 = 3; i5 = 2; break;
1100 : }
1101 0 : setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1102 0 : cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1103 0 : cols[i5][0], cols[i5][1]);
1104 : // Swap for qbar q incoming
1105 0 : if (id1 < 0) swapColAcol();
1106 :
1107 0 : }
1108 :
1109 : //==========================================================================
1110 :
1111 : // Sigma3qg2qqqbarDiff class.
1112 : // Cross section for q g -> q q' qbar'
1113 : // Crossed relation from q q' -> q q' g, q != q'
1114 : // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1115 :
1116 : //--------------------------------------------------------------------------
1117 :
1118 : // Initialize process.
1119 :
1120 : void Sigma3qg2qqqbarDiff::initProc() {
1121 :
1122 : // Read number of quarks to be considered in massless approximation.
1123 0 : nQuarkNew = settingsPtr->mode("HardQCD:nQuarkNew");
1124 :
1125 0 : }
1126 :
1127 : //--------------------------------------------------------------------------
1128 :
1129 : // Evaluate |M|^2 - no incoming flavour dependence
1130 :
1131 : void Sigma3qg2qqqbarDiff::sigmaKin() {
1132 :
1133 : // Pick a final state configuration
1134 0 : pickFinal();
1135 :
1136 : // gq or qg incoming
1137 0 : for (int i = 0; i < 2; i++) {
1138 :
1139 : // Map incoming four-vectors
1140 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1141 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1142 0 : mapFinal();
1143 :
1144 : // Crossing (note extra -ve sign in total sigma)
1145 0 : swap(pCM[i], pCM[4]);
1146 0 : pCM[i] = -pCM[i];
1147 0 : pCM[4] = -pCM[4];
1148 :
1149 : // |M|^2
1150 : // Extra factor of (6) from picking a final state configuration
1151 : // Extra factor of (3 / 8) as averaging over incoming gluon
1152 : // Extra factor of (nQuarkNew - 1) due to new q/qbar pair
1153 0 : sigma[i] = -6. * (3. / 8.) * (nQuarkNew - 1) * m2Calc();
1154 :
1155 : }
1156 :
1157 0 : }
1158 :
1159 : //--------------------------------------------------------------------------
1160 :
1161 : // Evaluate |M|^2 - incoming flavour dependence
1162 :
1163 : double Sigma3qg2qqqbarDiff::sigmaHat() {
1164 : // gq or qg incoming
1165 0 : return (id1 == 21) ? sigma[0] : sigma[1];
1166 : }
1167 :
1168 : //--------------------------------------------------------------------------
1169 :
1170 : // Select identity, colour and anticolour.
1171 :
1172 : void Sigma3qg2qqqbarDiff::setIdColAcol(){
1173 : // Pick new q qbar flavour with incoming flavour disallowed
1174 0 : int sigmaIdx = (id1 == 21) ? 0 : 1;
1175 0 : int idIn = (id1 == 21) ? id2 : id1;
1176 0 : int idNew = 1 + int( (nQuarkNew - 1) * rndmPtr->flat() );
1177 0 : if (idNew >= abs(idIn)) ++idNew;
1178 :
1179 : // qbar instead of q incoming means swap outgoing q/qbar pair
1180 0 : int id3Tmp = idIn, id4Tmp = idNew, id5Tmp = -idNew;
1181 0 : if (idIn < 0) swap(id4Tmp, id5Tmp);
1182 : // If g q incoming rather than q g, idIn and idNew
1183 : // should be exchanged (see sigmaKin)
1184 0 : if (sigmaIdx == 0) swap(id3Tmp, id4Tmp);
1185 : // Outgoing flavours; now just map as if q g incoming
1186 0 : switch (config) {
1187 0 : case 0: id3 = id3Tmp; id4 = id4Tmp; id5 = id5Tmp; break;
1188 0 : case 1: id3 = id3Tmp; id4 = id5Tmp; id5 = id4Tmp; break;
1189 0 : case 2: id3 = id4Tmp; id4 = id3Tmp; id5 = id5Tmp; break;
1190 0 : case 3: id3 = id5Tmp; id4 = id3Tmp; id5 = id4Tmp; break;
1191 0 : case 4: id3 = id4Tmp; id4 = id5Tmp; id5 = id3Tmp; break;
1192 0 : case 5: id3 = id5Tmp; id4 = id4Tmp; id5 = id3Tmp; break;
1193 : }
1194 0 : setId(id1, id2, id3, id4, id5);
1195 :
1196 : // Temporary solution; start with either
1197 : // g q1 -> q1 q2 qbar2
1198 : // g qbar1 -> qbar1 qbar2 q2
1199 0 : int cols[5][2];
1200 0 : cols[0][0] = 1; cols[0][1] = 2;
1201 0 : if (idIn > 0) {
1202 0 : cols[1][0] = 3; cols[1][1] = 0;
1203 0 : cols[2][0] = 1; cols[2][1] = 0;
1204 0 : cols[3][0] = 3; cols[3][1] = 0;
1205 0 : cols[4][0] = 0; cols[4][1] = 2;
1206 0 : } else {
1207 0 : cols[1][0] = 0; cols[1][1] = 3;
1208 0 : cols[2][0] = 0; cols[2][1] = 2;
1209 0 : cols[3][0] = 0; cols[3][1] = 3;
1210 0 : cols[4][0] = 1; cols[4][1] = 0;
1211 : }
1212 : // Swap incoming if q/qbar g instead
1213 0 : if (id2 == 21) {
1214 0 : swap(cols[0][0], cols[1][0]);
1215 0 : swap(cols[0][1], cols[1][1]);
1216 0 : }
1217 : // Map final state
1218 : int i3 = 0, i4 = 0, i5 = 0;
1219 0 : if (sigmaIdx == 0) {
1220 0 : switch (config) {
1221 0 : case 0: i3 = 3; i4 = 2; i5 = 4; break;
1222 0 : case 1: i3 = 3; i4 = 4; i5 = 2; break;
1223 0 : case 2: i3 = 2; i4 = 3; i5 = 4; break;
1224 0 : case 3: i3 = 4; i4 = 3; i5 = 2; break;
1225 0 : case 4: i3 = 2; i4 = 4; i5 = 3; break;
1226 0 : case 5: i3 = 4; i4 = 2; i5 = 3; break;
1227 : }
1228 : } else {
1229 0 : switch (config) {
1230 0 : case 0: i3 = 2; i4 = 3; i5 = 4; break;
1231 0 : case 1: i3 = 2; i4 = 4; i5 = 3; break;
1232 0 : case 2: i3 = 3; i4 = 2; i5 = 4; break;
1233 0 : case 3: i3 = 4; i4 = 2; i5 = 3; break;
1234 0 : case 4: i3 = 3; i4 = 4; i5 = 2; break;
1235 0 : case 5: i3 = 4; i4 = 3; i5 = 2; break;
1236 : }
1237 : }
1238 0 : setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1239 0 : cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1240 0 : cols[i5][0], cols[i5][1]);
1241 0 : }
1242 :
1243 : //==========================================================================
1244 :
1245 : // Sigma3qq2qqgSame class.
1246 : // Cross section for q q' -> q q' g, q == q'.
1247 :
1248 : //--------------------------------------------------------------------------
1249 :
1250 : // Evaluate |M|^2 - no incoming flavour dependence
1251 :
1252 : void Sigma3qq2qqgSame::sigmaKin() {
1253 : // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1254 :
1255 : // Incoming four-vectors
1256 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1257 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1258 : // Pick/map a final state configuration
1259 0 : pickFinal();
1260 0 : mapFinal();
1261 :
1262 : // |M|^2
1263 : // Extra factor (3) from picking final state configuration
1264 : // (original answer already has a factor 2 from identical
1265 : // quarks in the final state)
1266 0 : sigma = 3. * m2Calc();
1267 :
1268 0 : }
1269 :
1270 : //--------------------------------------------------------------------------
1271 :
1272 : // |M|^2
1273 :
1274 : inline double Sigma3qq2qqgSame::m2Calc() {
1275 :
1276 : // Four-products
1277 0 : s = (pCM[0] + pCM[1]).m2Calc();
1278 0 : t = (pCM[0] - pCM[2]).m2Calc();
1279 0 : u = (pCM[0] - pCM[3]).m2Calc();
1280 0 : sp = (pCM[2] + pCM[3]).m2Calc();
1281 0 : tp = (pCM[1] - pCM[3]).m2Calc();
1282 0 : up = (pCM[1] - pCM[2]).m2Calc();
1283 :
1284 : // |M|^2
1285 0 : ssp = s * sp;
1286 0 : ttp = t * tp;
1287 0 : uup = u * up;
1288 0 : s_sp = s + sp;
1289 0 : t_tp = t + tp;
1290 0 : u_up = u + up;
1291 :
1292 0 : double den1 = (pCM[0] * pCM[4]) * (pCM[1] * pCM[4]) *
1293 0 : (pCM[2] * pCM[4]) * (pCM[3] * pCM[4]);
1294 :
1295 0 : double fac1 = s * (t * u + tp * up) + sp * (t * up + tp * u);
1296 0 : double fac2 = ssp - ttp - uup;
1297 0 : double fac3 = 2. * (ttp * u_up + uup * t_tp);
1298 :
1299 0 : double num1 = u_up * (ssp + ttp - uup) + fac1;
1300 0 : double num2 = s_sp * fac2 + fac3;
1301 0 : double num3 = (s * s + sp * sp + u * u + up * up) / (t * tp);
1302 :
1303 0 : double num4 = t_tp * (ssp - ttp + uup) + fac1;
1304 0 : double num5 = (s * s + sp * sp + t * t + tp * tp) / (u * up);
1305 :
1306 0 : double num6 = s_sp * fac2 - fac3 - 2. * fac1;
1307 0 : double num7 = (s * s + sp * sp) * fac2;
1308 0 : double den7 = (ttp * uup);
1309 :
1310 : // C1 = (N^2 - 1)^2 / 4N^3 = 16. / 27.
1311 : // C2 = (N^2 - 1) / 4N^3 = 2. / 27.
1312 : // C3 = (N^4 - 1) / 8N^4 = 10. / 81.
1313 : // C4 = (N^2 - 1)^2 / 8N^4 = 8. / 81.
1314 0 : return (1. / 8.) * pow3(4. * M_PI * alpS) *
1315 0 : ( ( (16. / 27.) * num1 - (2. / 27.) * num2 ) * num3 +
1316 0 : ( (16. / 27.) * num4 - (2. / 27.) * num2 ) * num5 +
1317 0 : ( (10. / 81.) * num2 + (8. / 81.) * num6 ) *
1318 0 : ( num7 / den7 ) ) / den1;
1319 :
1320 : }
1321 :
1322 : //--------------------------------------------------------------------------
1323 :
1324 : // Evaluate |M|^2 - incoming flavour dependence
1325 :
1326 : double Sigma3qq2qqgSame::sigmaHat() {
1327 : // q q / qbar qbar incoming states only
1328 0 : if (id1 != id2) return 0.;
1329 0 : return sigma;
1330 0 : }
1331 :
1332 : //--------------------------------------------------------------------------
1333 :
1334 : // Select identity, colour and anticolour.
1335 :
1336 : void Sigma3qq2qqgSame::setIdColAcol(){
1337 :
1338 : // Need to know where the gluon was mapped (pCM[4])
1339 : int gIdx = 0;
1340 0 : switch (config) {
1341 0 : case 3: case 5: gIdx = 0; break;
1342 0 : case 1: case 4: gIdx = 1; break;
1343 0 : case 0: case 2: gIdx = 2; break;
1344 : }
1345 :
1346 : // Outgoing flavours
1347 0 : int idTmp[3] = { id1, id1, id1 };
1348 0 : idTmp[gIdx] = 21;
1349 0 : setId(id1, id2, idTmp[0], idTmp[1], idTmp[2]);
1350 :
1351 : // Temporary solution; start with q q -> q q g
1352 0 : setColAcol(1, 0, 2, 0, 1, 0, 3, 0, 2, 3);
1353 : // Map gluon
1354 0 : swap( colSave[5], colSave[gIdx + 3]);
1355 0 : swap(acolSave[5], acolSave[gIdx + 3]);
1356 : // Swap if qbar qbar incoming
1357 0 : if (id1 < 0) swapColAcol();
1358 :
1359 0 : }
1360 :
1361 : //--------------------------------------------------------------------------
1362 :
1363 : // Map a final state configuration
1364 : inline void Sigma3qq2qqgSame::mapFinal() {
1365 0 : switch (config) {
1366 0 : case 0: pCM[2] = p3cm; pCM[3] = p4cm; pCM[4] = p5cm; break;
1367 0 : case 1: pCM[2] = p3cm; pCM[3] = p5cm; pCM[4] = p4cm; break;
1368 0 : case 2: pCM[2] = p4cm; pCM[3] = p3cm; pCM[4] = p5cm; break;
1369 0 : case 3: pCM[2] = p4cm; pCM[3] = p5cm; pCM[4] = p3cm; break;
1370 0 : case 4: pCM[2] = p5cm; pCM[3] = p3cm; pCM[4] = p4cm; break;
1371 0 : case 5: pCM[2] = p5cm; pCM[3] = p4cm; pCM[4] = p3cm; break;
1372 : }
1373 0 : }
1374 :
1375 : //==========================================================================
1376 :
1377 : // Sigma3qqbar2qqbargSame class.
1378 : // Cross section for q qbar -> q qbar g
1379 : // Crossed relation from q(bar) q(bar) -> q(bar) q(bar) g:
1380 : // q1(p+) q2(p-) -> q1(q+) q2(q-) g(k)
1381 : //--------------------------------------------------------------------------
1382 :
1383 : // Evaluate |M|^2 - no incoming flavour dependence
1384 :
1385 : void Sigma3qqbar2qqbargSame::sigmaKin() {
1386 :
1387 : // Incoming four-vectors
1388 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1389 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1390 :
1391 : // Pick and map a final state configuration
1392 0 : pickFinal();
1393 0 : mapFinal();
1394 :
1395 : // Crossing
1396 0 : swap(pCM[1], pCM[3]);
1397 0 : pCM[1] = -pCM[1];
1398 0 : pCM[3] = -pCM[3];
1399 :
1400 : // |M|^2
1401 : // Extra factor of (6) from picking a final state configuration
1402 0 : sigma = 6. * m2Calc();
1403 :
1404 0 : }
1405 :
1406 : //--------------------------------------------------------------------------
1407 :
1408 : // Select identity, colour and anticolour.
1409 :
1410 : void Sigma3qqbar2qqbargSame::setIdColAcol(){
1411 : // Outgoing flavours; easiest to map by hand
1412 0 : switch (config) {
1413 0 : case 0: id3 = id1; id4 = id2; id5 = 21; break;
1414 0 : case 1: id3 = id1; id4 = 21; id5 = id2; break;
1415 0 : case 2: id3 = id2; id4 = id1; id5 = 21; break;
1416 0 : case 3: id3 = 21; id4 = id1; id5 = id2; break;
1417 0 : case 4: id3 = id2; id4 = 21; id5 = id1; break;
1418 0 : case 5: id3 = 21; id4 = id2; id5 = id1; break;
1419 : }
1420 0 : setId(id1, id2, id3, id4, id5);
1421 :
1422 : // Temporary solution; start with q qbar -> q qbar g
1423 0 : int cols[5][2];
1424 0 : cols[0][0] = 1; cols[0][1] = 0;
1425 0 : cols[1][0] = 0; cols[1][1] = 2;
1426 0 : cols[2][0] = 1; cols[2][1] = 0;
1427 0 : cols[3][0] = 0; cols[3][1] = 3;
1428 0 : cols[4][0] = 3; cols[4][1] = 2;
1429 : // Map final state
1430 : int i3 = 0, i4 = 0, i5 = 0;
1431 0 : switch (config) {
1432 0 : case 0: i3 = 2; i4 = 3; i5 = 4; break;
1433 0 : case 1: i3 = 2; i4 = 4; i5 = 3; break;
1434 0 : case 2: i3 = 3; i4 = 2; i5 = 4; break;
1435 0 : case 3: i3 = 4; i4 = 2; i5 = 3; break;
1436 0 : case 4: i3 = 3; i4 = 4; i5 = 2; break;
1437 0 : case 5: i3 = 4; i4 = 3; i5 = 2; break;
1438 : }
1439 0 : setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1440 0 : cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1441 0 : cols[i5][0], cols[i5][1]);
1442 : // Swap for qbar q incoming
1443 0 : if (id1 < 0) swapColAcol();
1444 0 : }
1445 :
1446 : //==========================================================================
1447 :
1448 : // Sigma3qg2qqqbarSame class.
1449 : // Cross section for q g -> q q qbar.
1450 : // Crossed relation from q(bar) q(bar) -> q(bar) q(bar) g:
1451 : // q1(p+) q1(p-) -> q1(q+) q1(q-) g(k)
1452 :
1453 : //--------------------------------------------------------------------------
1454 :
1455 : // Evaluate |M|^2 - no incoming flavour dependence
1456 :
1457 : void Sigma3qg2qqqbarSame::sigmaKin() {
1458 :
1459 : // Pick a final state configuration
1460 0 : pickFinal();
1461 :
1462 : // gq and qg incoming
1463 0 : for (int i = 0; i < 2; i++) {
1464 :
1465 : // Map incoming four-vectors
1466 0 : pCM[0] = Vec4( 0., 0., 0.5 * mH, 0.5 * mH);
1467 0 : pCM[1] = Vec4( 0., 0., -0.5 * mH, 0.5 * mH);
1468 0 : mapFinal();
1469 :
1470 : // Crossing (note extra -ve sign in total sigma)
1471 0 : swap(pCM[i], pCM[4]);
1472 0 : pCM[i] = -pCM[i];
1473 0 : pCM[4] = -pCM[4];
1474 :
1475 : // |M|^2
1476 : // XXX - Extra factor of (3) from picking a final state configuration???
1477 : // Extra factor of (3 / 8) as averaging over incoming gluon
1478 0 : sigma[i] = -3. * (3. / 8.) * m2Calc();
1479 :
1480 : }
1481 :
1482 0 : }
1483 :
1484 : //--------------------------------------------------------------------------
1485 :
1486 : // Evaluate |M|^2 - incoming flavour dependence
1487 :
1488 : double Sigma3qg2qqqbarSame::sigmaHat() {
1489 : // gq or qg incoming
1490 0 : return (id1 == 21) ? sigma[0] : sigma[1];
1491 : }
1492 :
1493 : //--------------------------------------------------------------------------
1494 :
1495 : // Select identity, colour and anticolour.
1496 :
1497 : void Sigma3qg2qqqbarSame::setIdColAcol(){
1498 :
1499 : // Pick outgoing flavour configuration
1500 0 : int idIn = (id1 == 21) ? id2 : id1;
1501 :
1502 : // Outgoing flavours; easiest just to map by hand
1503 0 : switch (config) {
1504 0 : case 0: id3 = idIn; id4 = idIn; id5 = -idIn; break;
1505 0 : case 1: id3 = idIn; id4 = -idIn; id5 = idIn; break;
1506 0 : case 2: id3 = idIn; id4 = idIn; id5 = -idIn; break;
1507 0 : case 3: id3 = -idIn; id4 = idIn; id5 = idIn; break;
1508 0 : case 4: id3 = idIn; id4 = -idIn; id5 = idIn; break;
1509 0 : case 5: id3 = -idIn; id4 = idIn; id5 = idIn; break;
1510 : }
1511 0 : setId(id1, id2, id3, id4, id5);
1512 :
1513 : // Temporary solution; start with either
1514 : // g q1 -> q1 q2 qbar2
1515 : // g qbar1 -> qbar1 qbar2 q2
1516 0 : int cols[5][2];
1517 0 : cols[0][0] = 1; cols[0][1] = 2;
1518 0 : if (idIn > 0) {
1519 0 : cols[1][0] = 3; cols[1][1] = 0;
1520 0 : cols[2][0] = 1; cols[2][1] = 0;
1521 0 : cols[3][0] = 3; cols[3][1] = 0;
1522 0 : cols[4][0] = 0; cols[4][1] = 2;
1523 0 : } else {
1524 0 : cols[1][0] = 0; cols[1][1] = 3;
1525 0 : cols[2][0] = 0; cols[2][1] = 2;
1526 0 : cols[3][0] = 0; cols[3][1] = 3;
1527 0 : cols[4][0] = 1; cols[4][1] = 0;
1528 : }
1529 : // Swap incoming if q/qbar g instead
1530 0 : if (id2 == 21) {
1531 0 : swap(cols[0][0], cols[1][0]);
1532 0 : swap(cols[0][1], cols[1][1]);
1533 0 : }
1534 : // Map final state
1535 : int i3 = 0, i4 = 0, i5 = 0;
1536 0 : switch (config) {
1537 0 : case 0: i3 = 2; i4 = 3; i5 = 4; break;
1538 0 : case 1: i3 = 2; i4 = 4; i5 = 3; break;
1539 0 : case 2: i3 = 3; i4 = 2; i5 = 4; break;
1540 0 : case 3: i3 = 4; i4 = 2; i5 = 3; break;
1541 0 : case 4: i3 = 3; i4 = 4; i5 = 2; break;
1542 0 : case 5: i3 = 4; i4 = 3; i5 = 2; break;
1543 : }
1544 0 : setColAcol(cols[0][0], cols[0][1], cols[1][0], cols[1][1],
1545 0 : cols[i3][0], cols[i3][1], cols[i4][0], cols[i4][1],
1546 0 : cols[i5][0], cols[i5][1]);
1547 0 : }
1548 :
1549 : //==========================================================================
1550 :
1551 : } // end namespace Pythia8
|
