Line data Source code
1 : // SigmaEW.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 : // electroweak simulation classes (including photon processes).
8 :
9 : #include "Pythia8/SigmaEW.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // Sigma2qg2qgamma class.
16 : // Cross section for q g -> q gamma.
17 :
18 : //--------------------------------------------------------------------------
19 :
20 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
21 :
22 : void Sigma2qg2qgamma::sigmaKin() {
23 :
24 : // Calculate kinematics dependence.
25 0 : sigUS = (1./3.) * (sH2 + uH2) / (-sH * uH);
26 :
27 : // Answer.
28 0 : sigma0 = (M_PI/sH2) * alpS * alpEM * sigUS;
29 :
30 0 : }
31 :
32 : //--------------------------------------------------------------------------
33 :
34 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
35 :
36 : double Sigma2qg2qgamma::sigmaHat() {
37 :
38 : // Incoming flavour gives charge factor.
39 0 : int idNow = (id2 == 21) ? id1 : id2;
40 0 : double eNow = couplingsPtr->ef( abs(idNow) );
41 0 : return sigma0 * pow2(eNow);
42 :
43 0 : }
44 :
45 : //--------------------------------------------------------------------------
46 :
47 : // Select identity, colour and anticolour.
48 :
49 : void Sigma2qg2qgamma::setIdColAcol() {
50 :
51 : // Construct outgoing flavours.
52 0 : id3 = (id1 == 21) ? 22 : id1;
53 0 : id4 = (id2 == 21) ? 22 : id2;
54 0 : setId( id1, id2, id3, id4);
55 :
56 : // Colour flow topology. Swap if first is gluon, or when antiquark.
57 0 : setColAcol( 1, 0, 2, 1, 2, 0, 0, 0);
58 0 : if (id1 == 21) swapCol1234();
59 0 : if (id1 < 0 || id2 < 0) swapColAcol();
60 :
61 0 : }
62 :
63 : //==========================================================================
64 :
65 : // Sigma2qqbar2ggamma class.
66 : // Cross section for q qbar -> g gamma.
67 :
68 : //--------------------------------------------------------------------------
69 :
70 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
71 :
72 : void Sigma2qqbar2ggamma::sigmaKin() {
73 :
74 : // Calculate kinematics dependence.
75 0 : double sigTU = (8./9.) * (tH2 + uH2) / (tH * uH);
76 :
77 : // Answer.
78 0 : sigma0 = (M_PI/sH2) * alpS * alpEM * sigTU;
79 :
80 0 : }
81 :
82 : //--------------------------------------------------------------------------
83 :
84 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
85 :
86 : double Sigma2qqbar2ggamma::sigmaHat() {
87 :
88 : // Incoming flavour gives charge factor.
89 0 : double eNow = couplingsPtr->ef( abs(id1) );
90 0 : return sigma0 * pow2(eNow);
91 :
92 0 : }
93 :
94 : //--------------------------------------------------------------------------
95 :
96 : // Select identity, colour and anticolour.
97 :
98 : void Sigma2qqbar2ggamma::setIdColAcol() {
99 :
100 : // Outgoing flavours trivial.
101 0 : setId( id1, id2, 21, 22);
102 :
103 : // One colour flow topology. Swap if first is antiquark.
104 0 : setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
105 0 : if (id1 < 0) swapColAcol();
106 :
107 0 : }
108 :
109 : //==========================================================================
110 :
111 : // Sigma2gg2ggamma class.
112 : // Cross section for g g -> g gamma.
113 : // Proceeds through a quark box, by default using 5 massless quarks.
114 :
115 : //--------------------------------------------------------------------------
116 :
117 : // Initialize process, especially parton-flux object.
118 :
119 : void Sigma2gg2ggamma::initProc() {
120 :
121 : // Maximum quark flavour in loop.
122 0 : int nQuarkLoop = settingsPtr->mode("PromptPhoton:nQuarkLoop");
123 :
124 : // Calculate charge factor from the allowed quarks in the box.
125 0 : chargeSum = - 1./3. + 2./3. - 1./3.;
126 0 : if (nQuarkLoop >= 4) chargeSum += 2./3.;
127 0 : if (nQuarkLoop >= 5) chargeSum -= 1./3.;
128 0 : if (nQuarkLoop >= 6) chargeSum += 2./3.;
129 :
130 0 : }
131 :
132 : //--------------------------------------------------------------------------
133 :
134 : // Evaluate d(sigmaHat)/d(tHat).
135 :
136 : void Sigma2gg2ggamma::sigmaKin() {
137 :
138 : // Logarithms of Mandelstam variable ratios.
139 0 : double logST = log( -sH / tH );
140 0 : double logSU = log( -sH / uH );
141 0 : double logTU = log( tH / uH );
142 :
143 : // Real and imaginary parts of separate amplitudes.
144 0 : double b0stuRe = 1. + (tH - uH) / sH * logTU
145 0 : + 0.5 * (tH2 + uH2) / sH2 * (pow2(logTU) + pow2(M_PI));
146 0 : double b0stuIm = 0.;
147 0 : double b0tsuRe = 1. + (sH - uH) / tH * logSU
148 0 : + 0.5 * (sH2 + uH2) / tH2 * pow2(logSU);
149 0 : double b0tsuIm = -M_PI * ( (sH - uH) / tH + (sH2 + uH2) / tH2 * logSU);
150 0 : double b0utsRe = 1. + (sH - tH) / uH * logST
151 0 : + 0.5 * (sH2 + tH2) / uH2 * pow2(logST);
152 0 : double b0utsIm = -M_PI * ( (sH - tH) / uH + (sH2 + tH2) / uH2 * logST);
153 0 : double b1stuRe = -1.;
154 0 : double b1stuIm = 0.;
155 0 : double b2stuRe = -1.;
156 0 : double b2stuIm = 0.;
157 :
158 : // Calculate kinematics dependence.
159 0 : double sigBox = pow2(b0stuRe) + pow2(b0stuIm) + pow2(b0tsuRe)
160 0 : + pow2(b0tsuIm) + pow2(b0utsRe) + pow2(b0utsIm) + 4. * pow2(b1stuRe)
161 0 : + 4. * pow2(b1stuIm) + pow2(b2stuRe) + pow2(b2stuIm);
162 :
163 : // Answer.
164 0 : sigma = (5. / (192. * M_PI * sH2)) * pow2(chargeSum)
165 0 : * pow3(alpS) * alpEM * sigBox;
166 :
167 0 : }
168 :
169 : //--------------------------------------------------------------------------
170 :
171 : // Select identity, colour and anticolour.
172 :
173 : void Sigma2gg2ggamma::setIdColAcol() {
174 :
175 : // Flavours and colours are trivial.
176 0 : setId( id1, id2, 21, 22);
177 0 : setColAcol( 1, 2, 2, 3, 1, 3, 0, 0);
178 0 : if (rndmPtr->flat() > 0.5) swapColAcol();
179 :
180 0 : }
181 :
182 : //==========================================================================
183 :
184 : // Sigma2ffbar2gammagamma class.
185 : // Cross section for q qbar -> gamma gamma.
186 :
187 : //--------------------------------------------------------------------------
188 :
189 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
190 :
191 : void Sigma2ffbar2gammagamma::sigmaKin() {
192 :
193 : // Calculate kinematics dependence.
194 0 : sigTU = 2. * (tH2 + uH2) / (tH * uH);
195 :
196 : // Answer contains factor 1/2 from identical photons.
197 0 : sigma0 = (M_PI/sH2) * pow2(alpEM) * 0.5 * sigTU;
198 :
199 0 : }
200 :
201 : //--------------------------------------------------------------------------
202 :
203 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
204 :
205 : double Sigma2ffbar2gammagamma::sigmaHat() {
206 :
207 : // Incoming flavour gives charge and colour factors.
208 0 : double eNow = couplingsPtr->ef( abs(id1) );
209 0 : double colFac = (abs(id1) < 9) ? 1. / 3. : 1.;
210 0 : return sigma0 * pow4(eNow) * colFac;
211 :
212 0 : }
213 :
214 : //--------------------------------------------------------------------------
215 :
216 : // Select identity, colour and anticolour.
217 :
218 : void Sigma2ffbar2gammagamma::setIdColAcol() {
219 :
220 : // Outgoing flavours trivial.
221 0 : setId( id1, id2, 22, 22);
222 :
223 : // No colours at all or one flow topology. Swap if first is antiquark.
224 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
225 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
226 0 : if (id1 < 0) swapColAcol();
227 :
228 0 : }
229 :
230 : //==========================================================================
231 :
232 : // Sigma2gg2gammagamma class.
233 : // Cross section for g g -> gamma gamma.
234 : // Proceeds through a quark box, by default using 5 massless quarks.
235 :
236 : //--------------------------------------------------------------------------
237 :
238 : // Initialize process.
239 :
240 : void Sigma2gg2gammagamma::initProc() {
241 :
242 : // Maximum quark flavour in loop.
243 0 : int nQuarkLoop = settingsPtr->mode("PromptPhoton:nQuarkLoop");
244 :
245 : // Calculate charge factor from the allowed quarks in the box.
246 0 : charge2Sum = 1./9. + 4./9. + 1./9.;
247 0 : if (nQuarkLoop >= 4) charge2Sum += 4./9.;
248 0 : if (nQuarkLoop >= 5) charge2Sum += 1./9.;
249 0 : if (nQuarkLoop >= 6) charge2Sum += 4./9.;
250 :
251 0 : }
252 :
253 : //--------------------------------------------------------------------------
254 :
255 : // Evaluate d(sigmaHat)/d(tHat).
256 :
257 : void Sigma2gg2gammagamma::sigmaKin() {
258 :
259 : // Logarithms of Mandelstam variable ratios.
260 0 : double logST = log( -sH / tH );
261 0 : double logSU = log( -sH / uH );
262 0 : double logTU = log( tH / uH );
263 :
264 : // Real and imaginary parts of separate amplitudes.
265 0 : double b0stuRe = 1. + (tH - uH) / sH * logTU
266 0 : + 0.5 * (tH2 + uH2) / sH2 * (pow2(logTU) + pow2(M_PI));
267 0 : double b0stuIm = 0.;
268 0 : double b0tsuRe = 1. + (sH - uH) / tH * logSU
269 0 : + 0.5 * (sH2 + uH2) / tH2 * pow2(logSU);
270 0 : double b0tsuIm = -M_PI * ( (sH - uH) / tH + (sH2 + uH2) / tH2 * logSU);
271 0 : double b0utsRe = 1. + (sH - tH) / uH * logST
272 0 : + 0.5 * (sH2 + tH2) / uH2 * pow2(logST);
273 0 : double b0utsIm = -M_PI * ( (sH - tH) / uH + (sH2 + tH2) / uH2 * logST);
274 0 : double b1stuRe = -1.;
275 0 : double b1stuIm = 0.;
276 0 : double b2stuRe = -1.;
277 0 : double b2stuIm = 0.;
278 :
279 : // Calculate kinematics dependence.
280 0 : double sigBox = pow2(b0stuRe) + pow2(b0stuIm) + pow2(b0tsuRe)
281 0 : + pow2(b0tsuIm) + pow2(b0utsRe) + pow2(b0utsIm) + 4. * pow2(b1stuRe)
282 0 : + 4. * pow2(b1stuIm) + pow2(b2stuRe) + pow2(b2stuIm);
283 :
284 : // Answer contains factor 1/2 from identical photons.
285 0 : sigma = (0.5 / (16. * M_PI * sH2)) * pow2(charge2Sum)
286 0 : * pow2(alpS) * pow2(alpEM) * sigBox;
287 :
288 0 : }
289 :
290 : //--------------------------------------------------------------------------
291 :
292 : // Select identity, colour and anticolour.
293 :
294 : void Sigma2gg2gammagamma::setIdColAcol() {
295 :
296 : // Flavours and colours are trivial.
297 0 : setId( id1, id2, 22, 22);
298 0 : setColAcol( 1, 2, 2, 1, 0, 0, 0, 0);
299 :
300 0 : }
301 :
302 : //==========================================================================
303 :
304 : // Sigma2ff2fftgmZ class.
305 : // Cross section for f f' -> f f' via t-channel gamma*/Z0 exchange
306 : // (f is quark or lepton).
307 :
308 : //--------------------------------------------------------------------------
309 :
310 : // Initialize process.
311 :
312 : void Sigma2ff2fftgmZ::initProc() {
313 :
314 : // Store Z0 mass for propagator. Common coupling factor.
315 0 : gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
316 0 : mZ = particleDataPtr->m0(23);
317 0 : mZS = mZ*mZ;
318 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
319 0 : * couplingsPtr->cos2thetaW());
320 :
321 0 : }
322 :
323 : //--------------------------------------------------------------------------
324 :
325 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
326 :
327 : void Sigma2ff2fftgmZ::sigmaKin() {
328 :
329 : // Cross section part common for all incoming flavours.
330 0 : double sigma0 = (M_PI / sH2) * pow2(alpEM);
331 :
332 : // Kinematical functions for gamma-gamma, gamma-Z and Z-Z parts.
333 0 : sigmagmgm = sigma0 * 2. * (sH2 + uH2) / tH2;
334 0 : sigmagmZ = sigma0 * 4. * thetaWRat * sH2 / (tH * (tH - mZS));
335 0 : sigmaZZ = sigma0 * 2. * pow2(thetaWRat) * sH2 / pow2(tH - mZS);
336 0 : if (gmZmode == 1) {sigmagmZ = 0.; sigmaZZ = 0.;}
337 0 : if (gmZmode == 2) {sigmagmgm = 0.; sigmagmZ = 0.;}
338 :
339 0 : }
340 :
341 : //--------------------------------------------------------------------------
342 :
343 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
344 :
345 : double Sigma2ff2fftgmZ::sigmaHat() {
346 :
347 : // Couplings for current flavour combination.
348 0 : int id1Abs = abs(id1);
349 0 : double e1 = couplingsPtr->ef(id1Abs);
350 0 : double v1 = couplingsPtr->vf(id1Abs);
351 0 : double a1 = couplingsPtr->af(id1Abs);
352 0 : int id2Abs = abs(id2);
353 0 : double e2 = couplingsPtr->ef(id2Abs);
354 0 : double v2 = couplingsPtr->vf(id2Abs);
355 0 : double a2 = couplingsPtr->af(id2Abs);
356 :
357 : // Distinguish same-sign and opposite-sign fermions.
358 0 : double epsi = (id1 * id2 > 0) ? 1. : -1.;
359 :
360 : // Flavour-dependent cross section.
361 0 : double sigma = sigmagmgm * pow2(e1 * e2)
362 0 : + sigmagmZ * e1 * e2 * (v1 * v2 * (1. + uH2 / sH2)
363 0 : + a1 * a2 * epsi * (1. - uH2 / sH2))
364 0 : + sigmaZZ * ((v1*v1 + a1*a1) * (v2*v2 + a2*a2) * (1. + uH2 / sH2)
365 0 : + 4. * v1 * a1 * v2 * a2 * epsi * (1. - uH2 / sH2));
366 :
367 : // Spin-state extra factor 2 per incoming neutrino.
368 0 : if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
369 0 : if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
370 :
371 : // Answer.
372 0 : return sigma;
373 :
374 : }
375 :
376 : //--------------------------------------------------------------------------
377 :
378 : // Select identity, colour and anticolour.
379 :
380 : void Sigma2ff2fftgmZ::setIdColAcol() {
381 :
382 : // Trivial flavours: out = in.
383 0 : setId( id1, id2, id1, id2);
384 :
385 : // Colour flow topologies. Swap when antiquarks.
386 0 : if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
387 0 : setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
388 0 : else if (abs(id1) < 9 && abs(id2) < 9)
389 0 : setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
390 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
391 0 : else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
392 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
393 0 : if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
394 0 : swapColAcol();
395 :
396 0 : }
397 :
398 : //==========================================================================
399 :
400 : // Sigma2ff2fftW class.
401 : // Cross section for f_1 f_2 -> f_3 f_4 via t-channel W+- exchange
402 : // (f is quark or lepton).
403 :
404 : //--------------------------------------------------------------------------
405 :
406 : // Initialize process.
407 :
408 : void Sigma2ff2fftW::initProc() {
409 :
410 : // Store W+- mass for propagator. Common coupling factor.
411 0 : mW = particleDataPtr->m0(24);
412 0 : mWS = mW*mW;
413 0 : thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
414 :
415 0 : }
416 :
417 : //--------------------------------------------------------------------------
418 :
419 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
420 :
421 : void Sigma2ff2fftW::sigmaKin() {
422 :
423 : // Cross section part common for all incoming flavours.
424 0 : sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat)
425 0 : * 4. * sH2 / pow2(tH - mWS);
426 :
427 0 : }
428 :
429 : //--------------------------------------------------------------------------
430 :
431 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
432 :
433 : double Sigma2ff2fftW::sigmaHat() {
434 :
435 : // Some flavour combinations not possible.
436 0 : int id1Abs = abs(id1);
437 0 : int id2Abs = abs(id2);
438 0 : if ( (id1Abs%2 == id2Abs%2 && id1 * id2 > 0)
439 0 : || (id1Abs%2 != id2Abs%2 && id1 * id2 < 0) ) return 0.;
440 :
441 : // Basic cross section.
442 0 : double sigma = sigma0;
443 0 : if (id1 * id2 < 0) sigma *= uH2 / sH2;
444 :
445 : // CKM factors for final states.
446 0 : sigma *= couplingsPtr->V2CKMsum(id1Abs) * couplingsPtr->V2CKMsum(id2Abs);
447 :
448 : // Spin-state extra factor 2 per incoming neutrino.
449 0 : if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
450 0 : if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
451 :
452 : // Answer.
453 : return sigma;
454 :
455 0 : }
456 :
457 : //--------------------------------------------------------------------------
458 :
459 : // Select identity, colour and anticolour.
460 :
461 : void Sigma2ff2fftW::setIdColAcol() {
462 :
463 : // Pick out-flavours by relative CKM weights.
464 0 : id3 = couplingsPtr->V2CKMpick(id1);
465 0 : id4 = couplingsPtr->V2CKMpick(id2);
466 0 : setId( id1, id2, id3, id4);
467 :
468 : // Colour flow topologies. Swap when antiquarks.
469 0 : if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
470 0 : setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
471 0 : else if (abs(id1) < 9 && abs(id2) < 9)
472 0 : setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
473 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 1, 0, 0, 0);
474 0 : else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
475 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
476 0 : if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
477 0 : swapColAcol();
478 :
479 0 : }
480 :
481 :
482 : //==========================================================================
483 :
484 : // Sigma2qq2QqtW class.
485 : // Cross section for q q' -> Q q" via t-channel W+- exchange.
486 : // Related to Sigma2ff2ffViaW class, but with massive matrix elements.
487 :
488 : //--------------------------------------------------------------------------
489 :
490 : // Initialize process.
491 :
492 : void Sigma2qq2QqtW::initProc() {
493 :
494 : // Process name.
495 0 : nameSave = "q q -> Q q (t-channel W+-)";
496 0 : if (idNew == 4) nameSave = "q q -> c q (t-channel W+-)";
497 0 : if (idNew == 5) nameSave = "q q -> b q (t-channel W+-)";
498 0 : if (idNew == 6) nameSave = "q q -> t q (t-channel W+-)";
499 0 : if (idNew == 7) nameSave = "q q -> b' q (t-channel W+-)";
500 0 : if (idNew == 8) nameSave = "q q -> t' q (t-channel W+-)";
501 :
502 : // Store W+- mass for propagator. Common coupling factor.
503 0 : mW = particleDataPtr->m0(24);
504 0 : mWS = mW*mW;
505 0 : thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
506 :
507 : // Secondary open width fractions, relevant for top (or heavier).
508 0 : openFracPos = particleDataPtr->resOpenFrac(idNew);
509 0 : openFracNeg = particleDataPtr->resOpenFrac(-idNew);
510 :
511 0 : }
512 :
513 : //--------------------------------------------------------------------------
514 :
515 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
516 :
517 : void Sigma2qq2QqtW::sigmaKin() {
518 :
519 : // Cross section part common for all incoming flavours.
520 0 : sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat) * 4. / pow2(tH - mWS);
521 :
522 0 : }
523 :
524 : //--------------------------------------------------------------------------
525 :
526 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
527 :
528 : double Sigma2qq2QqtW::sigmaHat() {
529 :
530 : // Some flavour combinations not possible.
531 0 : int id1Abs = abs(id1);
532 0 : int id2Abs = abs(id2);
533 0 : bool diff12 = (id1Abs%2 != id2Abs%2);
534 0 : if ( (!diff12 && id1 * id2 > 0)
535 0 : || ( diff12 && id1 * id2 < 0) ) return 0.;
536 :
537 : // Basic cross section.
538 0 : double sigma = sigma0;
539 0 : sigma *= (id1 * id2 > 0) ? sH * (sH - s3) : uH * (uH - s3);
540 :
541 : // Secondary width if t or tbar produced on either side.
542 0 : double openFrac1 = (id1 > 0) ? openFracPos : openFracNeg;
543 0 : double openFrac2 = (id2 > 0) ? openFracPos : openFracNeg;
544 :
545 : // CKM factors for final states; further impossible case.
546 0 : bool diff1N = (id1Abs%2 != idNew%2);
547 0 : bool diff2N = (id2Abs%2 != idNew%2);
548 0 : if (diff1N && diff2N)
549 0 : sigma *= ( couplingsPtr->V2CKMid(id1Abs, idNew) * openFrac1
550 0 : * couplingsPtr->V2CKMsum(id2Abs) + couplingsPtr->V2CKMsum(id1Abs)
551 0 : * couplingsPtr->V2CKMid(id2Abs, idNew) * openFrac2 );
552 0 : else if (diff1N)
553 0 : sigma *= couplingsPtr->V2CKMid(id1Abs, idNew) * openFrac1
554 0 : * couplingsPtr->V2CKMsum(id2Abs);
555 0 : else if (diff2N)
556 0 : sigma *= couplingsPtr->V2CKMsum(id1Abs)
557 0 : * couplingsPtr->V2CKMid(id2Abs, idNew) * openFrac2;
558 : else sigma = 0.;
559 :
560 : // Spin-state extra factor 2 per incoming neutrino.
561 0 : if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
562 0 : if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
563 :
564 : // Answer.
565 : return sigma;
566 :
567 0 : }
568 :
569 : //--------------------------------------------------------------------------
570 :
571 : // Select identity, colour and anticolour.
572 :
573 : void Sigma2qq2QqtW::setIdColAcol() {
574 :
575 : // For topologies like d dbar -> (t/c/u) (t/c/u)bar pick side.
576 0 : int id1Abs = abs(id1);
577 0 : int id2Abs = abs(id2);
578 : int side = 1;
579 0 : if ( (id1Abs + idNew)%2 == 1 && (id2Abs + idNew)%2 == 1 ) {
580 0 : double prob1 = couplingsPtr->V2CKMid(id1Abs, idNew)
581 0 : * couplingsPtr->V2CKMsum(id2Abs);
582 0 : prob1 *= (id1 > 0) ? openFracPos : openFracNeg;
583 0 : double prob2 = couplingsPtr->V2CKMid(id2Abs, idNew)
584 0 : * couplingsPtr->V2CKMsum(id1Abs);
585 0 : prob2 *= (id2 > 0) ? openFracPos : openFracNeg;
586 0 : if (prob2 > rndmPtr->flat() * (prob1 + prob2)) side = 2;
587 0 : }
588 0 : else if ((id2Abs + idNew)%2 == 1) side = 2;
589 :
590 : // Pick out-flavours by relative CKM weights.
591 0 : if (side == 1) {
592 : // q q' -> t q" : correct order from start.
593 0 : id3 = (id1 > 0) ? idNew : -idNew;
594 0 : id4 = couplingsPtr->V2CKMpick(id2);
595 0 : setId( id1, id2, id3, id4);
596 0 : } else {
597 : // q q' -> q" t : stored as t q" so swap tHat <-> uHat.
598 0 : swapTU = true;
599 0 : id3 = couplingsPtr->V2CKMpick(id1);
600 0 : id4 = (id2 > 0) ? idNew : -idNew;
601 0 : setId( id1, id2, id4, id3);
602 : }
603 :
604 : // Colour flow topologies. Swap when antiquarks on side 1.
605 0 : if (side == 1 && id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
606 0 : else if (id1 * id2 > 0) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
607 0 : else if (side == 1) setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
608 0 : else setColAcol( 1, 0, 0, 2, 0, 2, 1, 0);
609 0 : if (id1 < 0) swapColAcol();
610 :
611 0 : }
612 :
613 : //--------------------------------------------------------------------------
614 :
615 : // Evaluate weight for decay angles of W in top decay.
616 :
617 : double Sigma2qq2QqtW::weightDecay( Event& process, int iResBeg,
618 : int iResEnd) {
619 :
620 : // For top decay hand over to standard routine, else done.
621 0 : if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
622 0 : return weightTopDecay( process, iResBeg, iResEnd);
623 0 : else return 1.;
624 :
625 0 : }
626 :
627 : //==========================================================================
628 :
629 : // Sigma1ffbar2gmZ class.
630 : // Cross section for f fbar -> gamma*/Z0 (f is quark or lepton).
631 :
632 : //--------------------------------------------------------------------------
633 :
634 : // Initialize process.
635 :
636 : void Sigma1ffbar2gmZ::initProc() {
637 :
638 : // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
639 0 : gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
640 :
641 : // Store Z0 mass and width for propagator.
642 0 : mRes = particleDataPtr->m0(23);
643 0 : GammaRes = particleDataPtr->mWidth(23);
644 0 : m2Res = mRes*mRes;
645 0 : GamMRat = GammaRes / mRes;
646 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
647 0 : * couplingsPtr->cos2thetaW());
648 :
649 : // Set pointer to particle properties and decay table.
650 0 : particlePtr = particleDataPtr->particleDataEntryPtr(23);
651 :
652 0 : }
653 :
654 : //--------------------------------------------------------------------------
655 :
656 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
657 :
658 : void Sigma1ffbar2gmZ::sigmaKin() {
659 :
660 : // Common coupling factors.
661 0 : double colQ = 3. * (1. + alpS / M_PI);
662 :
663 : // Reset quantities to sum. Declare variables in loop.
664 0 : gamSum = 0.;
665 0 : intSum = 0.;
666 0 : resSum = 0.;
667 : int idAbs, onMode;
668 0 : double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
669 :
670 : // Loop over all Z0 decay channels.
671 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
672 0 : idAbs = abs( particlePtr->channel(i).product(0) );
673 :
674 : // Only contributions from three fermion generations, except top.
675 0 : if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
676 0 : mf = particleDataPtr->m0(idAbs);
677 :
678 : // Check that above threshold. Phase space.
679 0 : if (mH > 2. * mf + MASSMARGIN) {
680 0 : mr = pow2(mf / mH);
681 0 : betaf = sqrtpos(1. - 4. * mr);
682 0 : psvec = betaf * (1. + 2. * mr);
683 0 : psaxi = pow3(betaf);
684 :
685 : // Combine phase space with couplings.
686 0 : ef2 = couplingsPtr->ef2(idAbs) * psvec;
687 0 : efvf = couplingsPtr->efvf(idAbs) * psvec;
688 0 : vf2af2 = couplingsPtr->vf2(idAbs) * psvec
689 0 : + couplingsPtr->af2(idAbs) * psaxi;
690 0 : colf = (idAbs < 6) ? colQ : 1.;
691 :
692 : // Store sum of combinations. For outstate only open channels.
693 0 : onMode = particlePtr->channel(i).onMode();
694 0 : if (onMode == 1 || onMode == 2) {
695 0 : gamSum += colf * ef2;
696 0 : intSum += colf * efvf;
697 0 : resSum += colf * vf2af2;
698 0 : }
699 :
700 : // End loop over fermions.
701 : }
702 : }
703 : }
704 :
705 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
706 0 : gamProp = 4. * M_PI * pow2(alpEM) / (3. * sH);
707 0 : intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
708 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
709 0 : resProp = gamProp * pow2(thetaWRat * sH)
710 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
711 :
712 : // Optionally only keep gamma* or Z0 term.
713 0 : if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
714 0 : if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
715 :
716 0 : }
717 :
718 : //--------------------------------------------------------------------------
719 :
720 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
721 :
722 : double Sigma1ffbar2gmZ::sigmaHat() {
723 :
724 : // Combine gamma, interference and Z0 parts.
725 0 : int idAbs = abs(id1);
726 0 : double sigma = couplingsPtr->ef2(idAbs) * gamProp * gamSum
727 0 : + couplingsPtr->efvf(idAbs) * intProp * intSum
728 0 : + couplingsPtr->vf2af2(idAbs) * resProp * resSum;
729 :
730 : // Colour factor. Answer.
731 0 : if (idAbs < 9) sigma /= 3.;
732 0 : return sigma;
733 :
734 : }
735 :
736 : //--------------------------------------------------------------------------
737 :
738 : // Select identity, colour and anticolour.
739 :
740 : void Sigma1ffbar2gmZ::setIdColAcol() {
741 :
742 : // Flavours trivial.
743 0 : setId( id1, id2, 23);
744 :
745 : // Colour flow topologies. Swap when antiquarks.
746 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
747 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
748 0 : if (id1 < 0) swapColAcol();
749 :
750 0 : }
751 :
752 : //--------------------------------------------------------------------------
753 :
754 : // Evaluate weight for gamma*/Z0 decay angle.
755 :
756 : double Sigma1ffbar2gmZ::weightDecay( Event& process, int iResBeg,
757 : int iResEnd) {
758 :
759 : // Z should sit in entry 5.
760 0 : if (iResBeg != 5 || iResEnd != 5) return 1.;
761 :
762 : // Couplings for in- and out-flavours.
763 0 : int idInAbs = process[3].idAbs();
764 0 : double ei = couplingsPtr->ef(idInAbs);
765 0 : double vi = couplingsPtr->vf(idInAbs);
766 0 : double ai = couplingsPtr->af(idInAbs);
767 0 : int idOutAbs = process[6].idAbs();
768 0 : double ef = couplingsPtr->ef(idOutAbs);
769 0 : double vf = couplingsPtr->vf(idOutAbs);
770 0 : double af = couplingsPtr->af(idOutAbs);
771 :
772 : // Phase space factors. (One power of beta left out in formulae.)
773 0 : double mf = process[6].m();
774 0 : double mr = mf*mf / sH;
775 0 : double betaf = sqrtpos(1. - 4. * mr);
776 :
777 : // Coefficients of angular expression.
778 0 : double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
779 0 : + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
780 0 : double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
781 0 : + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
782 0 : double coefAsym = betaf * ( ei * ai * intProp * ef * af
783 0 : + 4. * vi * ai * resProp * vf * af );
784 :
785 : // Flip asymmetry for in-fermion + out-antifermion.
786 0 : if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
787 :
788 : // Reconstruct decay angle and weight for it.
789 0 : double cosThe = (process[3].p() - process[4].p())
790 0 : * (process[7].p() - process[6].p()) / (sH * betaf);
791 0 : double wtMax = 2. * (coefTran + abs(coefAsym));
792 0 : double wt = coefTran * (1. + pow2(cosThe))
793 0 : + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
794 :
795 : // Done.
796 0 : return (wt / wtMax);
797 :
798 0 : }
799 :
800 : //==========================================================================
801 :
802 : // Sigma1ffbar2W class.
803 : // Cross section for f fbar' -> W+- (f is quark or lepton).
804 :
805 : //--------------------------------------------------------------------------
806 :
807 : // Initialize process.
808 :
809 : void Sigma1ffbar2W::initProc() {
810 :
811 : // Store W+- mass and width for propagator.
812 0 : mRes = particleDataPtr->m0(24);
813 0 : GammaRes = particleDataPtr->mWidth(24);
814 0 : m2Res = mRes*mRes;
815 0 : GamMRat = GammaRes / mRes;
816 0 : thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
817 :
818 : // Set pointer to particle properties and decay table.
819 0 : particlePtr = particleDataPtr->particleDataEntryPtr(24);
820 :
821 0 : }
822 :
823 : //--------------------------------------------------------------------------
824 :
825 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
826 :
827 : void Sigma1ffbar2W::sigmaKin() {
828 :
829 : // Set up Breit-Wigner. Cross section for W+ and W- separately.
830 0 : double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
831 0 : double preFac = alpEM * thetaWRat * mH;
832 0 : sigma0Pos = preFac * sigBW * particlePtr->resWidthOpen(24, mH);
833 0 : sigma0Neg = preFac * sigBW * particlePtr->resWidthOpen(-24, mH);
834 :
835 0 : }
836 :
837 : //--------------------------------------------------------------------------
838 :
839 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
840 :
841 : double Sigma1ffbar2W::sigmaHat() {
842 :
843 : // Secondary width for W+ or W-. CKM and colour factors.
844 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
845 0 : double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
846 0 : if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
847 :
848 : // Answer.
849 0 : return sigma;
850 :
851 : }
852 :
853 : //--------------------------------------------------------------------------
854 :
855 : // Select identity, colour and anticolour.
856 :
857 : void Sigma1ffbar2W::setIdColAcol() {
858 :
859 : // Sign of outgoing W.
860 0 : int sign = 1 - 2 * (abs(id1)%2);
861 0 : if (id1 < 0) sign = -sign;
862 0 : setId( id1, id2, 24 * sign);
863 :
864 : // Colour flow topologies. Swap when antiquarks.
865 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
866 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
867 0 : if (id1 < 0) swapColAcol();
868 :
869 0 : }
870 :
871 : //--------------------------------------------------------------------------
872 :
873 : // Evaluate weight for W decay angle.
874 :
875 : double Sigma1ffbar2W::weightDecay( Event& process, int iResBeg,
876 : int iResEnd) {
877 :
878 : // W should sit in entry 5.
879 0 : if (iResBeg != 5 || iResEnd != 5) return 1.;
880 :
881 : // Phase space factors.
882 0 : double mr1 = pow2(process[6].m()) / sH;
883 0 : double mr2 = pow2(process[7].m()) / sH;
884 0 : double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
885 :
886 : // Sign of asymmetry.
887 0 : double eps = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
888 :
889 : // Reconstruct decay angle and weight for it.
890 0 : double cosThe = (process[3].p() - process[4].p())
891 0 : * (process[7].p() - process[6].p()) / (sH * betaf);
892 : double wtMax = 4.;
893 0 : double wt = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
894 :
895 : // Done.
896 0 : return (wt / wtMax);
897 :
898 0 : }
899 :
900 : //==========================================================================
901 :
902 : // Sigma2ffbar2ffbarsgm class.
903 : // Cross section f fbar -> gamma* -> f' fbar', for multiparton interactions.
904 :
905 : //--------------------------------------------------------------------------
906 :
907 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
908 :
909 : void Sigma2ffbar2ffbarsgm::sigmaKin() {
910 :
911 : // Pick new flavour. Allow three leptons and five quarks.
912 0 : double colQ = 1. + (alpS / M_PI);
913 0 : double flavWt = 3. + colQ * 11. / 3.;
914 0 : double flavRndm = rndmPtr->flat() * flavWt;
915 0 : if (flavRndm < 3.) {
916 0 : if (flavRndm < 1.) idNew = 11;
917 0 : else if (flavRndm < 2.) idNew = 13;
918 0 : else idNew = 15;
919 : } else {
920 0 : flavRndm = 3. * (flavRndm - 3.) / colQ;
921 0 : if (flavRndm < 4.) idNew = 2;
922 0 : else if (flavRndm < 8.) idNew = 4;
923 0 : else if (flavRndm < 9.) idNew = 1;
924 0 : else if (flavRndm < 10.) idNew = 3;
925 0 : else idNew = 5;
926 : }
927 0 : double mNew = particleDataPtr->m0(idNew);
928 0 : double m2New = mNew*mNew;
929 :
930 : // Calculate kinematics dependence. Give correct mass factors for
931 : // tHat, uHat defined as if massless kinematics, d(sigma)/d(Omega)
932 : // = beta (1 + cos^2(theta) + (1 - beta^2) sin^2(theta)).
933 : // Special case related to phase space form in multiparton interactions.
934 : double sigS = 0.;
935 0 : if (sH > 4. * m2New) {
936 0 : double beta = sqrt(1. - 4. * m2New / sH);
937 0 : sigS = beta * (2.* (tH2 + uH2) + 4. * (1. - beta * beta) * tH * uH)
938 0 : / sH2;
939 0 : }
940 :
941 : // Answer is proportional to number of outgoing flavours.
942 0 : sigma0 = (M_PI/sH2) * pow2(alpEM) * sigS * flavWt;
943 :
944 0 : }
945 :
946 : //--------------------------------------------------------------------------
947 :
948 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
949 :
950 : double Sigma2ffbar2ffbarsgm::sigmaHat() {
951 :
952 : // Charge and colour factors.
953 0 : double eNow = couplingsPtr->ef( abs(id1) );
954 0 : double sigma = sigma0 * pow2(eNow);
955 0 : if (abs(id1) < 9) sigma /= 3.;
956 :
957 : // Answer.
958 0 : return sigma;
959 :
960 0 : }
961 :
962 : //--------------------------------------------------------------------------
963 :
964 : // Select identity, colour and anticolour.
965 :
966 : void Sigma2ffbar2ffbarsgm::setIdColAcol() {
967 :
968 : // Set outgoing flavours.
969 0 : id3 = (id1 > 0) ? idNew : -idNew;
970 0 : setId( id1, id2, id3, -id3);
971 :
972 : // Colour flow topologies. Swap when antiquarks.
973 0 : if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
974 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
975 0 : else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
976 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
977 0 : if (id1 < 0) swapColAcol();
978 :
979 0 : }
980 :
981 : //==========================================================================
982 :
983 : // Sigma2ffbar2ffbarsgmZ class.
984 : // Cross section f fbar -> gamma*/Z0 -> f' fbar',
985 : // i.e. gamma*/Z0 decay as part of the hard process.
986 :
987 : //--------------------------------------------------------------------------
988 :
989 : // Initialize process.
990 :
991 : void Sigma2ffbar2ffbarsgmZ::initProc() {
992 :
993 : // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
994 0 : gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
995 :
996 : // Store Z0 mass and width for propagator.
997 0 : mRes = particleDataPtr->m0(23);
998 0 : GammaRes = particleDataPtr->mWidth(23);
999 0 : m2Res = mRes*mRes;
1000 0 : GamMRat = GammaRes / mRes;
1001 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1002 0 : * couplingsPtr->cos2thetaW());
1003 :
1004 : // Set pointer to particle properties and decay table.
1005 0 : particlePtr = particleDataPtr->particleDataEntryPtr(23);
1006 :
1007 0 : }
1008 :
1009 : //--------------------------------------------------------------------------
1010 :
1011 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1012 :
1013 : void Sigma2ffbar2ffbarsgmZ::sigmaKin() {
1014 :
1015 : // Common coupling factor.
1016 0 : colQ = 3. * (1. + alpS / M_PI);
1017 :
1018 : // Reset vectors and sums. Declare variables in loop.
1019 0 : idVec.resize(0);
1020 0 : gamT.resize(0);
1021 0 : gamL.resize(0);
1022 0 : intT.resize(0);
1023 0 : intL.resize(0);
1024 0 : intA.resize(0);
1025 0 : resT.resize(0);
1026 0 : resL.resize(0);
1027 0 : resA.resize(0);
1028 0 : gamSumT = 0.;
1029 0 : gamSumL = 0.;
1030 0 : intSumT = 0.;
1031 0 : intSumL = 0.;
1032 0 : intSumA = 0.;
1033 0 : resSumT = 0.;
1034 0 : resSumL = 0.;
1035 0 : resSumA = 0.;
1036 0 : int onMode, idAbs;
1037 0 : double mf, mr, betaf, ef, vf, af, colf, gamTf, gamLf, intTf, intLf,
1038 : intAf, resTf, resLf, resAf;
1039 :
1040 : // Loop over all Z0 decay channels.
1041 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1042 0 : onMode = particlePtr->channel(i).onMode();
1043 0 : idAbs = abs( particlePtr->channel(i).product(0) );
1044 :
1045 : // Only contributions from three fermion generations, except top.
1046 0 : if ( (onMode == 1 || onMode == 2) && ((idAbs > 0 && idAbs < 6)
1047 0 : || ( idAbs > 10 && idAbs < 17)) ) {
1048 0 : mf = particleDataPtr->m0(idAbs);
1049 :
1050 : // Check that channel above threshold. Phase space.
1051 0 : if (mH > 2. * mf + MASSMARGIN) {
1052 0 : mr = pow2(mf / mH);
1053 0 : betaf = sqrtpos(1. - 4. * mr);
1054 :
1055 : // Combine couplings (including colour) with phase space.
1056 0 : ef = couplingsPtr->ef(idAbs);
1057 0 : vf = couplingsPtr->vf(idAbs);
1058 0 : af = couplingsPtr->af(idAbs);
1059 0 : colf = (idAbs < 6) ? colQ : 1.;
1060 0 : gamTf = colf * ef * ef * betaf;
1061 0 : gamLf = gamTf * 4. * mr;
1062 0 : intTf = colf * ef * vf * betaf;
1063 0 : intLf = intTf * 4. * mr;
1064 0 : intAf = colf * ef * af * betaf;
1065 0 : resTf = colf * (vf * vf * betaf + af * af * pow3(betaf));
1066 0 : resLf = colf * vf * vf * betaf * 4. * mr;
1067 0 : resAf = colf * vf * af * betaf * 4.;
1068 :
1069 : // Store individual coplings and their sums.
1070 0 : idVec.push_back(idAbs);
1071 0 : gamT.push_back(gamTf);
1072 0 : gamL.push_back(gamLf);
1073 0 : intT.push_back(intTf);
1074 0 : intL.push_back(intLf);
1075 0 : intA.push_back(intAf);
1076 0 : resT.push_back(resTf);
1077 0 : resL.push_back(resLf);
1078 0 : resA.push_back(resAf);
1079 0 : gamSumT += gamTf;
1080 0 : gamSumL += gamLf;
1081 0 : intSumT += intTf;
1082 0 : intSumL += intLf;
1083 0 : intSumA += intAf;
1084 0 : resSumT += resTf;
1085 0 : resSumL += resLf;
1086 0 : resSumA += resAf;
1087 :
1088 : // End loop over Z0 decay channels.
1089 0 : }
1090 : }
1091 : }
1092 :
1093 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
1094 0 : gamProp = M_PI * pow2(alpEM) / sH2;
1095 0 : intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1096 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1097 0 : resProp = gamProp * pow2(thetaWRat * sH)
1098 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1099 :
1100 : // Optionally only keep gamma* or Z0 term.
1101 0 : if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1102 0 : if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1103 :
1104 : // Scattering angle in subsystem rest frame.
1105 0 : cThe = (tH - uH) / sH;
1106 :
1107 0 : }
1108 :
1109 : //--------------------------------------------------------------------------
1110 :
1111 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1112 :
1113 : double Sigma2ffbar2ffbarsgmZ::sigmaHat() {
1114 :
1115 : // Couplings for current in-flavour.
1116 0 : int id1Abs = abs(id1);
1117 0 : double ei = couplingsPtr->ef(id1Abs);
1118 0 : double vi = couplingsPtr->vf(id1Abs);
1119 0 : double ai = couplingsPtr->af(id1Abs);
1120 :
1121 : // Coefficients of angular expression.
1122 0 : double coefT = ei*ei * gamProp * gamSumT + ei*vi * intProp * intSumT
1123 0 : + (vi*vi + ai*ai) * resProp * resSumT;
1124 0 : double coefL = ei*ei * gamProp * gamSumL + ei*vi * intProp * intSumL
1125 0 : + (vi*vi + ai*ai) * resProp * resSumL;
1126 0 : double coefA = ei*ai * intProp * intSumA + vi*ai * resProp * resSumA;
1127 :
1128 : // Colour factor. Answer.
1129 0 : double sigma = coefT * (1. + pow2(cThe)) + coefL * (1. - pow2(cThe))
1130 0 : + 2. * coefA * cThe;
1131 0 : if (id1Abs < 9) sigma /= 3.;
1132 0 : return sigma;
1133 :
1134 : }
1135 :
1136 : //--------------------------------------------------------------------------
1137 :
1138 : // Select identity, colour and anticolour.
1139 :
1140 : void Sigma2ffbar2ffbarsgmZ::setIdColAcol() {
1141 :
1142 : // Couplings for chosen in-flavour.
1143 0 : int id1Abs = abs(id1);
1144 0 : double ei = couplingsPtr->ef(id1Abs);
1145 0 : double vi = couplingsPtr->vf(id1Abs);
1146 0 : double ai = couplingsPtr->af(id1Abs);
1147 :
1148 : // Contribution from each allowed out-flavour.
1149 0 : sigTLA.resize(0);
1150 0 : for (int i = 0; i < int(idVec.size()); ++i) {
1151 0 : double coefT = ei*ei * gamProp * gamT[i] + ei*vi * intProp * intT[i]
1152 0 : + (vi*vi + ai*ai) * resProp * resT[i];
1153 0 : double coefL = ei*ei * gamProp * gamL[i] + ei*vi * intProp * intL[i]
1154 0 : + (vi*vi + ai*ai) * resProp * resL[i];
1155 0 : double coefA = ei*ai * intProp * intA[i] + vi*ai * resProp * resA[i];
1156 0 : double sigma = coefT * (1. + pow2(cThe)) + coefL * (1. - pow2(cThe))
1157 0 : + 2. * coefA * cThe;
1158 0 : sigTLA.push_back(sigma);
1159 0 : }
1160 :
1161 : // Pick outgoing flavours.
1162 0 : int idNew = idVec[ rndmPtr->pick(sigTLA) ];
1163 0 : id3 = (id1 > 0) ? idNew : -idNew;
1164 0 : setId( id1, id2, id3, -id3);
1165 :
1166 : // Colour flow topologies. Swap when antiquarks.
1167 0 : if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1168 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1169 0 : else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1170 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1171 0 : if (id1 < 0) swapColAcol();
1172 :
1173 0 : }
1174 :
1175 : //==========================================================================
1176 :
1177 : // Sigma2ffbar2ffbarsW class.
1178 : // Cross section f_1 fbar_2 -> W+- -> f_3 fbar_4,
1179 : // i.e. W decay as part of the hard process.
1180 :
1181 : //--------------------------------------------------------------------------
1182 :
1183 : // Initialize process.
1184 :
1185 : void Sigma2ffbar2ffbarsW::initProc() {
1186 :
1187 : // Store W+- mass and width for propagator.
1188 0 : mRes = particleDataPtr->m0(24);
1189 0 : GammaRes = particleDataPtr->mWidth(24);
1190 0 : m2Res = mRes*mRes;
1191 0 : GamMRat = GammaRes / mRes;
1192 0 : thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1193 :
1194 : // Set pointer to particle properties and decay table.
1195 0 : particlePtr = particleDataPtr->particleDataEntryPtr(24);
1196 :
1197 0 : }
1198 :
1199 : //--------------------------------------------------------------------------
1200 :
1201 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1202 :
1203 : void Sigma2ffbar2ffbarsW::sigmaKin() {
1204 :
1205 : // Set up Breit-Wigner. Common sigmaHat cross section for W+ and W-.
1206 0 : double sigBW = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1207 0 : double preFac = alpEM * thetaWRat * mH;
1208 0 : sigma0 = preFac * sigBW * particlePtr->resWidthOpen(24, mH);
1209 :
1210 : // Convert to d(sigmaHat)/d(tHat). Fixed so integrates to sigmaHat.
1211 0 : sigma0 *= 3. * uH2 / (sH2 * sH);
1212 :
1213 : // Pick a decay channel.
1214 0 : if (!particlePtr->preparePick(24, mH)) {
1215 0 : sigma0 = 0.;
1216 0 : return;
1217 : }
1218 0 : DecayChannel& channel = particlePtr->pickChannel();
1219 0 : id3New = channel.product(0);
1220 0 : id4New = channel.product(1);
1221 :
1222 0 : }
1223 :
1224 : //--------------------------------------------------------------------------
1225 :
1226 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1227 :
1228 : double Sigma2ffbar2ffbarsW::sigmaHat() {
1229 :
1230 : // Secondary width for W+-. CKM and colour factors.
1231 0 : double sigma = sigma0;
1232 0 : if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1233 :
1234 : // Answer.
1235 0 : return sigma;
1236 :
1237 : }
1238 :
1239 : //--------------------------------------------------------------------------
1240 :
1241 : // Select identity, colour and anticolour.
1242 :
1243 : void Sigma2ffbar2ffbarsW::setIdColAcol() {
1244 :
1245 : // Find sign of W and its decay products. Set outgoing flavours.
1246 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
1247 0 : id3 = (idUp > 0) ? id3New : -id3New;
1248 0 : id4 = (idUp > 0) ? id4New : -id4New;
1249 0 : if (id1 * id3 < 0) swap(id3, id4);
1250 0 : setId( id1, id2, id3, id4);
1251 :
1252 : // Colour flow topologies. Swap when antifermions in 1 and 3.
1253 0 : if (abs(id1) < 9 && abs(id3) < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1254 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1255 0 : else if (abs(id3) < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1256 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1257 0 : if (id1 < 0) swapColAcol();
1258 :
1259 0 : }
1260 :
1261 : //==========================================================================
1262 :
1263 : // Sigma2ffbar2FFbarsgmZ class.
1264 : // Cross section f fbar -> gamma*/Z0 -> F Fbar.
1265 :
1266 : //--------------------------------------------------------------------------
1267 :
1268 : // Initialize process.
1269 :
1270 : void Sigma2ffbar2FFbarsgmZ::initProc() {
1271 :
1272 : // Process name.
1273 0 : nameSave = "f fbar -> F Fbar (s-channel gamma*/Z0)";
1274 0 : if (idNew == 4) nameSave = "f fbar -> c cbar (s-channel gamma*/Z0)";
1275 0 : if (idNew == 5) nameSave = "f fbar -> b bbar (s-channel gamma*/Z0)";
1276 0 : if (idNew == 6) nameSave = "f fbar -> t tbar (s-channel gamma*/Z0)";
1277 0 : if (idNew == 7) nameSave = "f fbar -> b' b'bar (s-channel gamma*/Z0)";
1278 0 : if (idNew == 8) nameSave = "f fbar -> t' t'bar (s-channel gamma*/Z0)";
1279 0 : if (idNew == 15) nameSave = "f fbar -> tau+ tau- (s-channel gamma*/Z0)";
1280 0 : if (idNew == 17) nameSave = "f fbar -> tau'+ tau'- (s-channel gamma*/Z0)";
1281 0 : if (idNew == 18)
1282 0 : nameSave = "f fbar -> nu'_tau nu'bar_tau (s-channel gamma*/Z0)";
1283 :
1284 : // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1285 0 : gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
1286 :
1287 : // Store Z0 mass and width for propagator.
1288 0 : mRes = particleDataPtr->m0(23);
1289 0 : GammaRes = particleDataPtr->mWidth(23);
1290 0 : m2Res = mRes*mRes;
1291 0 : GamMRat = GammaRes / mRes;
1292 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1293 0 : * couplingsPtr->cos2thetaW());
1294 :
1295 : // Store couplings of F.
1296 0 : ef = couplingsPtr->ef(idNew);
1297 0 : vf = couplingsPtr->vf(idNew);
1298 0 : af = couplingsPtr->af(idNew);
1299 :
1300 : // Secondary open width fraction, relevant for top (or heavier).
1301 0 : openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
1302 :
1303 0 : }
1304 :
1305 : //--------------------------------------------------------------------------
1306 :
1307 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1308 :
1309 : void Sigma2ffbar2FFbarsgmZ::sigmaKin() {
1310 :
1311 : // Check that above threshold.
1312 0 : isPhysical = true;
1313 0 : if (mH < m3 + m4 + MASSMARGIN) {
1314 0 : isPhysical = false;
1315 0 : return;
1316 : }
1317 :
1318 : // Define average F, Fbar mass so same beta. Phase space.
1319 0 : double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1320 0 : mr = s34Avg / sH;
1321 0 : betaf = sqrtpos(1. - 4. * mr);
1322 :
1323 : // Final-state colour factor.
1324 0 : double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) : 1.;
1325 :
1326 : // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1327 0 : cosThe = (tH - uH) / (betaf * sH);
1328 :
1329 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
1330 0 : gamProp = colF * M_PI * pow2(alpEM) / sH2;
1331 0 : intProp = gamProp * 2. * thetaWRat * sH * (sH - m2Res)
1332 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1333 0 : resProp = gamProp * pow2(thetaWRat * sH)
1334 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1335 :
1336 : // Optionally only keep gamma* or Z0 term.
1337 0 : if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
1338 0 : if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
1339 :
1340 0 : }
1341 :
1342 : //--------------------------------------------------------------------------
1343 :
1344 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1345 :
1346 : double Sigma2ffbar2FFbarsgmZ::sigmaHat() {
1347 :
1348 : // Fail if below threshold.
1349 0 : if (!isPhysical) return 0.;
1350 :
1351 : // Couplings for in-flavours.
1352 0 : int idAbs = abs(id1);
1353 0 : double ei = couplingsPtr->ef(idAbs);
1354 0 : double vi = couplingsPtr->vf(idAbs);
1355 0 : double ai = couplingsPtr->af(idAbs);
1356 :
1357 : // Coefficients of angular expression.
1358 0 : double coefTran = ei*ei * gamProp * ef*ef + ei * vi * intProp * ef * vf
1359 0 : + (vi*vi + ai*ai) * resProp * (vf*vf + pow2(betaf) * af*af);
1360 0 : double coefLong = 4. * mr * ( ei*ei * gamProp * ef*ef
1361 0 : + ei * vi * intProp * ef * vf + (vi*vi + ai*ai) * resProp * vf*vf );
1362 0 : double coefAsym = betaf * ( ei * ai * intProp * ef * af
1363 0 : + 4. * vi * ai * resProp * vf * af );
1364 :
1365 : // Combine gamma, interference and Z0 parts.
1366 0 : double sigma = coefTran * (1. + pow2(cosThe))
1367 0 : + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
1368 :
1369 : // Top: corrections for closed decay channels.
1370 0 : sigma *= openFracPair;
1371 :
1372 : // Initial-state colour factor. Answer.
1373 0 : if (idAbs < 9) sigma /= 3.;
1374 : return sigma;
1375 :
1376 0 : }
1377 :
1378 : //--------------------------------------------------------------------------
1379 :
1380 : // Select identity, colour and anticolour.
1381 :
1382 : void Sigma2ffbar2FFbarsgmZ::setIdColAcol() {
1383 :
1384 : // Set outgoing flavours.
1385 0 : id3 = (id1 > 0) ? idNew : -idNew;
1386 0 : setId( id1, id2, id3, -id3);
1387 :
1388 : // Colour flow topologies. Swap when antiquarks.
1389 0 : if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1390 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1391 0 : else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1392 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1393 0 : if (id1 < 0) swapColAcol();
1394 :
1395 0 : }
1396 :
1397 : //--------------------------------------------------------------------------
1398 :
1399 : // Evaluate weight for decay angles of W in top decay.
1400 :
1401 : double Sigma2ffbar2FFbarsgmZ::weightDecay( Event& process, int iResBeg,
1402 : int iResEnd) {
1403 :
1404 : // For top decay hand over to standard routine, else done.
1405 0 : if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1406 0 : return weightTopDecay( process, iResBeg, iResEnd);
1407 0 : else return 1.;
1408 :
1409 0 : }
1410 :
1411 : //==========================================================================
1412 :
1413 : // Sigma2ffbar2FfbarsW class.
1414 : // Cross section f fbar' -> W+- -> F fbar".
1415 :
1416 : //--------------------------------------------------------------------------
1417 :
1418 : // Initialize process.
1419 :
1420 : void Sigma2ffbar2FfbarsW::initProc() {
1421 :
1422 : // Process name.
1423 0 : nameSave = "f fbar -> F fbar (s-channel W+-)";
1424 0 : if (idNew == 4) nameSave = "f fbar -> c qbar (s-channel W+-)";
1425 0 : if (idNew == 5) nameSave = "f fbar -> b qbar (s-channel W+-)";
1426 0 : if (idNew == 6) nameSave = "f fbar -> t qbar (s-channel W+-)";
1427 0 : if (idNew == 7) nameSave = "f fbar -> b' qbar (s-channel W+-)";
1428 0 : if (idNew == 8) nameSave = "f fbar -> t' qbar (s-channel W+-)";
1429 0 : if (idNew == 7 && idNew2 == 6)
1430 0 : nameSave = "f fbar -> b' tbar (s-channel W+-)";
1431 0 : if (idNew == 8 && idNew2 == 7)
1432 0 : nameSave = "f fbar -> t' b'bar (s-channel W+-)";
1433 0 : if (idNew == 15 || idNew == 16)
1434 0 : nameSave = "f fbar -> tau nu_taubar (s-channel W+-)";
1435 0 : if (idNew == 17 || idNew == 18)
1436 0 : nameSave = "f fbar -> tau' nu'_taubar (s-channel W+-)";
1437 :
1438 : // Store W+- mass and width for propagator.
1439 0 : mRes = particleDataPtr->m0(24);
1440 0 : GammaRes = particleDataPtr->mWidth(24);
1441 0 : m2Res = mRes*mRes;
1442 0 : GamMRat = GammaRes / mRes;
1443 0 : thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1444 :
1445 : // For t/t' want to use at least b mass.
1446 0 : idPartner = idNew2;
1447 0 : if ( (idNew == 6 || idNew == 8) && idNew2 == 0 ) idPartner = 5;
1448 :
1449 : // Sum of CKM weights for quarks.
1450 0 : V2New = (idNew < 9) ? couplingsPtr->V2CKMsum(idNew) : 1.;
1451 0 : if (idNew2 != 0) V2New = couplingsPtr->V2CKMid(idNew, idNew2);
1452 :
1453 : // Secondary open width fractions, relevant for top or heavier.
1454 0 : openFracPos = particleDataPtr->resOpenFrac( idNew, -idNew2);
1455 0 : openFracNeg = particleDataPtr->resOpenFrac(-idNew, idNew2);
1456 :
1457 0 : }
1458 :
1459 : //--------------------------------------------------------------------------
1460 :
1461 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1462 :
1463 : void Sigma2ffbar2FfbarsW::sigmaKin() {
1464 :
1465 : // Check that above threshold.
1466 0 : isPhysical = true;
1467 0 : if (mH < m3 + m4 + MASSMARGIN) {
1468 0 : isPhysical = false;
1469 0 : return;
1470 : }
1471 :
1472 : // Phase space factors.
1473 0 : double mr1 = s3 / sH;
1474 0 : double mr2 = s4 / sH;
1475 0 : double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
1476 :
1477 : // Reconstruct decay angle so can reuse 2 -> 1 cross section.
1478 0 : double cosThe = (tH - uH) / (betaf * sH);
1479 :
1480 : // Set up Breit-Wigner and in- and out-widths.
1481 0 : double sigBW = 9. * M_PI * pow2(alpEM * thetaWRat)
1482 0 : / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1483 :
1484 : // Initial-state colour factor.
1485 0 : double colF = (idNew < 9) ? 3. * (1. + alpS / M_PI) * V2New : 1.;
1486 :
1487 : // Angular dependence.
1488 0 : double wt = pow2(1. + betaf * cosThe) - pow2(mr1 - mr2);
1489 :
1490 : // Temporary answer.
1491 0 : sigma0 = sigBW * colF * wt;
1492 :
1493 0 : }
1494 :
1495 : //--------------------------------------------------------------------------
1496 :
1497 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1498 :
1499 : double Sigma2ffbar2FfbarsW::sigmaHat() {
1500 :
1501 : // Fail if below threshold.
1502 0 : if (!isPhysical) return 0.;
1503 :
1504 : // CKM and colour factors.
1505 0 : double sigma = sigma0;
1506 0 : if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
1507 :
1508 : // Correction for secondary width in top (or heavier) decay.
1509 0 : int idSame = ((abs(id1) + idNew)%2 == 0) ? id1 : id2;
1510 0 : sigma *= (idSame > 0) ? openFracPos : openFracNeg;
1511 :
1512 : // Answer.
1513 : return sigma;
1514 :
1515 0 : }
1516 :
1517 : //--------------------------------------------------------------------------
1518 :
1519 : // Select identity, colour and anticolour.
1520 :
1521 : void Sigma2ffbar2FfbarsW::setIdColAcol() {
1522 :
1523 : // Set outgoing flavours.
1524 0 : id3 = idNew;
1525 0 : id4 = (idNew2 != 0) ? idNew2 : couplingsPtr->V2CKMpick(idNew);
1526 0 : if (idNew%2 == 0) {
1527 0 : int idInUp = (abs(id1)%2 == 0) ? id1 : id2;
1528 0 : if (idInUp > 0) id4 = -id4;
1529 0 : else id3 = -id3;
1530 0 : } else {
1531 0 : int idInDn = (abs(id1)%2 == 1) ? id1 : id2;
1532 0 : if (idInDn > 0) id4 = -id4;
1533 0 : else id3 = -id3;
1534 : }
1535 0 : setId( id1, id2, id3, id4);
1536 :
1537 : // Swap tHat and uHat for fbar' f -> F f".
1538 0 : if (id1 * id3 < 0) swapTU = true;
1539 :
1540 : // Colour flow topologies. Swap when antiquarks.
1541 0 : if (abs(id1) < 9 && idNew < 9) setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1542 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1543 0 : else if (idNew < 9) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
1544 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1545 0 : if (id1 < 0) swapCol12();
1546 0 : if (id3 < 0) swapCol34();
1547 :
1548 0 : }
1549 :
1550 : //--------------------------------------------------------------------------
1551 :
1552 : // Evaluate weight for decay angles of W in top decay.
1553 :
1554 : double Sigma2ffbar2FfbarsW::weightDecay( Event& process, int iResBeg,
1555 : int iResEnd) {
1556 :
1557 : // For top decay hand over to standard routine, else done.
1558 0 : if (idNew == 6 && process[process[iResBeg].mother1()].idAbs() == 6)
1559 0 : return weightTopDecay( process, iResBeg, iResEnd);
1560 0 : else return 1.;
1561 :
1562 0 : }
1563 :
1564 : //==========================================================================
1565 :
1566 : // Sigma2ffbargmZWgmZW class.
1567 : // Collects common methods for f fbar -> gamma*/Z0/W+- gamma*/Z0/W-+.
1568 :
1569 : //--------------------------------------------------------------------------
1570 :
1571 : // Calculate and store internal products.
1572 :
1573 : void Sigma2ffbargmZWgmZW::setupProd( Event& process, int i1, int i2,
1574 : int i3, int i4, int i5, int i6) {
1575 :
1576 : // Store incoming and outgoing momenta,
1577 0 : pRot[1] = process[i1].p();
1578 0 : pRot[2] = process[i2].p();
1579 0 : pRot[3] = process[i3].p();
1580 0 : pRot[4] = process[i4].p();
1581 0 : pRot[5] = process[i5].p();
1582 0 : pRot[6] = process[i6].p();
1583 :
1584 : // Do random rotation to avoid accidental zeroes in HA expressions.
1585 : bool smallPT = false;
1586 0 : do {
1587 : smallPT = false;
1588 0 : double thetaNow = acos(2. * rndmPtr->flat() - 1.);
1589 0 : double phiNow = 2. * M_PI * rndmPtr->flat();
1590 0 : for (int i = 1; i <= 6; ++i) {
1591 0 : pRot[i].rot( thetaNow, phiNow);
1592 0 : if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
1593 : }
1594 0 : } while (smallPT);
1595 :
1596 : // Calculate internal products.
1597 0 : for (int i = 1; i < 6; ++i) {
1598 0 : for (int j = i + 1; j <= 6; ++j) {
1599 0 : hA[i][j] =
1600 0 : sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
1601 0 : / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
1602 0 : - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
1603 0 : / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
1604 0 : hC[i][j] = conj( hA[i][j] );
1605 0 : if (i <= 2) {
1606 0 : hA[i][j] *= complex( 0., 1.);
1607 0 : hC[i][j] *= complex( 0., 1.);
1608 0 : }
1609 0 : hA[j][i] = - hA[i][j];
1610 0 : hC[j][i] = - hC[i][j];
1611 : }
1612 : }
1613 :
1614 0 : }
1615 :
1616 : //--------------------------------------------------------------------------
1617 :
1618 : // Evaluate the F function of Gunion and Kunszt.
1619 :
1620 : complex Sigma2ffbargmZWgmZW::fGK(int j1, int j2, int j3, int j4, int j5,
1621 : int j6) {
1622 :
1623 0 : return 4. * hA[j1][j3] * hC[j2][j6]
1624 0 : * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
1625 :
1626 : }
1627 :
1628 : //--------------------------------------------------------------------------
1629 :
1630 : // Evaluate the Xi function of Gunion and Kunszt.
1631 :
1632 : double Sigma2ffbargmZWgmZW::xiGK( double tHnow, double uHnow) {
1633 :
1634 0 : return - 4. * s3 * s4 + tHnow * (3. * tHnow + 4. * uHnow)
1635 0 : + tHnow * tHnow * ( tHnow * uHnow / (s3 * s4)
1636 0 : - 2. * (1. / s3 + 1./s4) * (tHnow + uHnow)
1637 0 : + 2. * (s3 / s4 + s4 / s3) );
1638 :
1639 : }
1640 :
1641 : //--------------------------------------------------------------------------
1642 :
1643 : // Evaluate the Xj function of Gunion and Kunszt.
1644 :
1645 : double Sigma2ffbargmZWgmZW::xjGK( double tHnow, double uHnow) {
1646 :
1647 0 : return 8. * pow2(s3 + s4) - 8. * (s3 + s4) * (tHnow + uHnow)
1648 0 : - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
1649 0 : / (s3 * s4) - 2. * (1. / s3 + 1. / s4) * (tHnow + uHnow)
1650 0 : + 2. * (s3 / s4 + s4 / s3) );
1651 :
1652 : }
1653 :
1654 : //==========================================================================
1655 :
1656 : // Sigma2ffbar2gmZgmZ class.
1657 : // Cross section for f fbar -> gamma*/Z0 gamma*/Z0 (f is quark or lepton).
1658 :
1659 : //--------------------------------------------------------------------------
1660 :
1661 : // Initialize process.
1662 :
1663 : void Sigma2ffbar2gmZgmZ::initProc() {
1664 :
1665 : // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
1666 0 : gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
1667 :
1668 : // Store Z0 mass and width for propagator.
1669 0 : mRes = particleDataPtr->m0(23);
1670 0 : GammaRes = particleDataPtr->mWidth(23);
1671 0 : m2Res = mRes*mRes;
1672 0 : GamMRat = GammaRes / mRes;
1673 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
1674 0 : * couplingsPtr->cos2thetaW());
1675 :
1676 : // Set pointer to particle properties and decay table.
1677 0 : particlePtr = particleDataPtr->particleDataEntryPtr(23);
1678 :
1679 0 : }
1680 :
1681 : //--------------------------------------------------------------------------
1682 :
1683 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1684 :
1685 : void Sigma2ffbar2gmZgmZ::sigmaKin() {
1686 :
1687 : // Cross section part common for all incoming flavours.
1688 0 : sigma0 = (M_PI / sH2) * pow2(alpEM) * 0.5
1689 0 : * ( (tH2 + uH2 + 2. * (s3 + s4) * sH) / (tH * uH)
1690 0 : - s3 * s4 * (1./tH2 + 1./uH2) );
1691 :
1692 : // Common coupling factors at the resonance masses
1693 0 : double alpEM3 = couplingsPtr->alphaEM(s3);
1694 0 : double alpS3 = couplingsPtr->alphaS(s3);
1695 0 : double colQ3 = 3. * (1. + alpS3 / M_PI);
1696 0 : double alpEM4 = couplingsPtr->alphaEM(s4);
1697 0 : double alpS4 = couplingsPtr->alphaS(s4);
1698 0 : double colQ4 = 3. * (1. + alpS4 / M_PI);
1699 :
1700 : // Reset quantities to sum. Declare variables in loop.
1701 0 : gamSum3 = 0.;
1702 0 : intSum3 = 0.;
1703 0 : resSum3 = 0.;
1704 0 : gamSum4 = 0.;
1705 0 : intSum4 = 0.;
1706 0 : resSum4 = 0.;
1707 : int onMode;
1708 0 : double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
1709 :
1710 : // Loop over all Z0 decay channels.
1711 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
1712 0 : int idAbs = abs( particlePtr->channel(i).product(0) );
1713 :
1714 : // Only contributions from three fermion generations, except top.
1715 0 : if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
1716 0 : mf = particleDataPtr->m0(idAbs);
1717 0 : onMode = particlePtr->channel(i).onMode();
1718 :
1719 : // First Z0: check that above threshold. Phase space.
1720 0 : if (m3 > 2. * mf + MASSMARGIN) {
1721 0 : mr = pow2(mf / m3);
1722 0 : betaf = sqrtpos(1. - 4. * mr);
1723 0 : psvec = betaf * (1. + 2. * mr);
1724 0 : psaxi = pow3(betaf);
1725 :
1726 : // First Z0: combine phase space with couplings.
1727 0 : ef2 = couplingsPtr->ef2(idAbs) * psvec;
1728 0 : efvf = couplingsPtr->efvf(idAbs) * psvec;
1729 0 : vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1730 0 : + couplingsPtr->af2(idAbs) * psaxi;
1731 0 : colf = (idAbs < 6) ? colQ3 : 1.;
1732 :
1733 : // First Z0: store sum of combinations for open outstate channels.
1734 0 : if (onMode == 1 || onMode == 2) {
1735 0 : gamSum3 += colf * ef2;
1736 0 : intSum3 += colf * efvf;
1737 0 : resSum3 += colf * vf2af2;
1738 0 : }
1739 : }
1740 :
1741 : // Second Z0: check that above threshold. Phase space.
1742 0 : if (m4 > 2. * mf + MASSMARGIN) {
1743 0 : mr = pow2(mf / m4);
1744 0 : betaf = sqrtpos(1. - 4. * mr);
1745 0 : psvec = betaf * (1. + 2. * mr);
1746 0 : psaxi = pow3(betaf);
1747 :
1748 : // Second Z0: combine phase space with couplings.
1749 0 : ef2 = couplingsPtr->ef2(idAbs) * psvec;
1750 0 : efvf = couplingsPtr->efvf(idAbs) * psvec;
1751 0 : vf2af2 = couplingsPtr->vf2(idAbs) * psvec
1752 0 : + couplingsPtr->af2(idAbs) * psaxi;
1753 0 : colf = (idAbs < 6) ? colQ4 : 1.;
1754 :
1755 : // Second Z0: store sum of combinations for open outstate channels.
1756 0 : if (onMode == 1 || onMode == 2) {
1757 0 : gamSum4 += colf * ef2;
1758 0 : intSum4 += colf * efvf;
1759 0 : resSum4 += colf * vf2af2;
1760 0 : }
1761 : }
1762 :
1763 : // End loop over fermions.
1764 : }
1765 : }
1766 :
1767 : // First Z0: calculate prefactors for gamma/interference/Z0 terms.
1768 0 : gamProp3 = 4. * alpEM3 / (3. * M_PI * s3);
1769 0 : intProp3 = gamProp3 * 2. * thetaWRat * s3 * (s3 - m2Res)
1770 0 : / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1771 0 : resProp3 = gamProp3 * pow2(thetaWRat * s3)
1772 0 : / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
1773 :
1774 : // First Z0: optionally only keep gamma* or Z0 term.
1775 0 : if (gmZmode == 1) {intProp3 = 0.; resProp3 = 0.;}
1776 0 : if (gmZmode == 2) {gamProp3 = 0.; intProp3 = 0.;}
1777 :
1778 : // Second Z0: calculate prefactors for gamma/interference/Z0 terms.
1779 0 : gamProp4 = 4. * alpEM4 / (3. * M_PI * s4);
1780 0 : intProp4 = gamProp4 * 2. * thetaWRat * s4 * (s4 - m2Res)
1781 0 : / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1782 0 : resProp4 = gamProp4 * pow2(thetaWRat * s4)
1783 0 : / ( pow2(s4 - m2Res) + pow2(s4 * GamMRat) );
1784 :
1785 : // Second Z0: optionally only keep gamma* or Z0 term.
1786 0 : if (gmZmode == 1) {intProp4 = 0.; resProp4 = 0.;}
1787 0 : if (gmZmode == 2) {gamProp4 = 0.; intProp4 = 0.;}
1788 :
1789 0 : }
1790 :
1791 : //--------------------------------------------------------------------------
1792 :
1793 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1794 :
1795 : double Sigma2ffbar2gmZgmZ::sigmaHat() {
1796 :
1797 : // Charge/2, left- and righthanded couplings for in-fermion.
1798 0 : int idAbs = abs(id1);
1799 0 : double ei = 0.5 * couplingsPtr->ef(idAbs);
1800 0 : double li = couplingsPtr->lf(idAbs);
1801 0 : double ri = couplingsPtr->rf(idAbs);
1802 :
1803 : // Combine left/right gamma, interference and Z0 parts for each Z0.
1804 0 : double left3 = ei * ei * gamProp3 * gamSum3
1805 0 : + ei * li * intProp3 * intSum3
1806 0 : + li * li * resProp3 * resSum3;
1807 : double right3 = ei * ei * gamProp3 * gamSum3
1808 0 : + ei * ri * intProp3 * intSum3
1809 0 : + ri * ri * resProp3 * resSum3;
1810 0 : double left4 = ei * ei * gamProp4 * gamSum4
1811 0 : + ei * li * intProp4 * intSum4
1812 0 : + li * li * resProp4 * resSum4;
1813 : double right4 = ei * ei * gamProp4 * gamSum4
1814 0 : + ei * ri * intProp4 * intSum4
1815 0 : + ri * ri * resProp4 * resSum4;
1816 :
1817 : // Combine left- and right-handed couplings for the two Z0's.
1818 0 : double sigma = sigma0 * (left3 * left4 + right3 * right4);
1819 :
1820 : // Correct for the running-width Z0 propagators weight in PhaseSpace.
1821 0 : sigma /= (runBW3 * runBW4);
1822 :
1823 : // Initial-state colour factor. Answer.
1824 0 : if (idAbs < 9) sigma /= 3.;
1825 0 : return sigma;
1826 :
1827 : }
1828 :
1829 : //--------------------------------------------------------------------------
1830 :
1831 : // Select identity, colour and anticolour.
1832 :
1833 : void Sigma2ffbar2gmZgmZ::setIdColAcol() {
1834 :
1835 : // Flavours trivial.
1836 0 : setId( id1, id2, 23, 23);
1837 :
1838 : // Colour flow topologies. Swap when antiquarks.
1839 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
1840 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
1841 0 : if (id1 < 0) swapColAcol();
1842 :
1843 0 : }
1844 :
1845 : //--------------------------------------------------------------------------
1846 :
1847 : // Evaluate correlated decay flavours of the two gamma*/Z0.
1848 : // Unique complication, caused by gamma*/Z0 mix different left/right.
1849 :
1850 : double Sigma2ffbar2gmZgmZ::weightDecayFlav( Event& process) {
1851 :
1852 : // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
1853 0 : i1 = (process[3].id() < 0) ? 3 : 4;
1854 0 : i2 = 7 - i1;
1855 0 : i3 = (process[7].id() > 0) ? 7 : 8;
1856 0 : i4 = 15 - i3;
1857 0 : i5 = (process[9].id() > 0) ? 9 : 10;
1858 0 : i6 = 19 - i5;
1859 :
1860 : // Charge/2, left- and righthanded couplings for in- and out-fermions.
1861 0 : int idAbs = process[i1].idAbs();
1862 0 : double ei = 0.5 * couplingsPtr->ef(idAbs);
1863 0 : double li = couplingsPtr->lf(idAbs);
1864 0 : double ri = couplingsPtr->rf(idAbs);
1865 0 : idAbs = process[i3].idAbs();
1866 0 : double e3 = 0.5 * couplingsPtr->ef(idAbs);
1867 0 : double l3 = couplingsPtr->lf(idAbs);
1868 0 : double r3 = couplingsPtr->rf(idAbs);
1869 0 : idAbs = process[i5].idAbs();
1870 0 : double e4 = 0.5 * couplingsPtr->ef(idAbs);
1871 0 : double l4 = couplingsPtr->lf(idAbs);
1872 0 : double r4 = couplingsPtr->rf(idAbs);
1873 :
1874 : // Left- and righthanded couplings combined with propagators.
1875 0 : c3LL = ei * ei * gamProp3 * e3 * e3
1876 0 : + ei * li * intProp3 * e3 * l3
1877 0 : + li * li * resProp3 * l3 * l3;
1878 0 : c3LR = ei * ei * gamProp3 * e3 * e3
1879 0 : + ei * li * intProp3 * e3 * r3
1880 0 : + li * li * resProp3 * r3 * r3;
1881 0 : c3RL = ei * ei * gamProp3 * e3 * e3
1882 0 : + ei * ri * intProp3 * e3 * l3
1883 0 : + ri * ri * resProp3 * l3 * l3;
1884 0 : c3RR = ei * ei * gamProp3 * e3 * e3
1885 0 : + ei * ri * intProp3 * e3 * r3
1886 0 : + ri * ri * resProp3 * r3 * r3;
1887 0 : c4LL = ei * ei * gamProp4 * e4 * e4
1888 0 : + ei * li * intProp4 * e4 * l4
1889 0 : + li * li * resProp4 * l4 * l4;
1890 0 : c4LR = ei * ei * gamProp4 * e4 * e4
1891 0 : + ei * li * intProp4 * e4 * r4
1892 0 : + li * li * resProp4 * r4 * r4;
1893 0 : c4RL = ei * ei * gamProp4 * e4 * e4
1894 0 : + ei * ri * intProp4 * e4 * l4
1895 0 : + ri * ri * resProp4 * l4 * l4;
1896 0 : c4RR = ei * ei * gamProp4 * e4 * e4
1897 0 : + ei * ri * intProp4 * e4 * r4
1898 0 : + ri * ri * resProp4 * r4 * r4;
1899 :
1900 : // Flavour weight and maximum.
1901 0 : flavWt = (c3LL + c3LR) * (c4LL + c4LR) + (c3RL + c3RR) * (c4RL + c4RR);
1902 0 : double flavWtMax = (c3LL + c3LR + c3RL + c3RR) * (c4LL + c4LR + c4RL + c4RR);
1903 :
1904 : // Done.
1905 0 : return flavWt / flavWtMax;
1906 :
1907 : }
1908 :
1909 : //--------------------------------------------------------------------------
1910 :
1911 : // Evaluate weight for decay angles of the two gamma*/Z0.
1912 :
1913 : double Sigma2ffbar2gmZgmZ::weightDecay( Event& process, int iResBeg,
1914 : int iResEnd) {
1915 :
1916 : // Two resonance decays, but with common weight.
1917 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
1918 :
1919 : // Set up four-products and internal products.
1920 0 : setupProd( process, i1, i2, i3, i4, i5, i6);
1921 :
1922 : // Flip tHat and uHat if first incoming is fermion.
1923 0 : double tHres = tH;
1924 0 : double uHres = uH;
1925 0 : if (process[3].id() > 0) swap( tHres, uHres);
1926 :
1927 : // Kinematics factors (norm(x) = |x|^2).
1928 0 : double fGK135 = norm( fGK( 1, 2, 3, 4, 5, 6) / tHres
1929 0 : + fGK( 1, 2, 5, 6, 3, 4) / uHres );
1930 0 : double fGK145 = norm( fGK( 1, 2, 4, 3, 5, 6) / tHres
1931 0 : + fGK( 1, 2, 5, 6, 4, 3) / uHres );
1932 0 : double fGK136 = norm( fGK( 1, 2, 3, 4, 6, 5) / tHres
1933 0 : + fGK( 1, 2, 6, 5, 3, 4) / uHres );
1934 0 : double fGK146 = norm( fGK( 1, 2, 4, 3, 6, 5) / tHres
1935 0 : + fGK( 1, 2, 6, 5, 4, 3) / uHres );
1936 0 : double fGK253 = norm( fGK( 2, 1, 5, 6, 3, 4) / tHres
1937 0 : + fGK( 2, 1, 3, 4, 5, 6) / uHres );
1938 0 : double fGK263 = norm( fGK( 2, 1, 6, 5, 3, 4) / tHres
1939 0 : + fGK( 2, 1, 3, 4, 6, 5) / uHres );
1940 0 : double fGK254 = norm( fGK( 2, 1, 5, 6, 4, 3) / tHres
1941 0 : + fGK( 2, 1, 4, 3, 5, 6) / uHres );
1942 0 : double fGK264 = norm( fGK( 2, 1, 6, 5, 4, 3) / tHres
1943 0 : + fGK( 2, 1, 4, 3, 6, 5) / uHres );
1944 :
1945 : // Weight and maximum.
1946 0 : double wt = c3LL * c4LL * fGK135 + c3LR * c4LL * fGK145
1947 0 : + c3LL * c4LR * fGK136 + c3LR * c4LR * fGK146
1948 0 : + c3RL * c4RL * fGK253 + c3RR * c4RL * fGK263
1949 0 : + c3RL * c4RR * fGK254 + c3RR * c4RR * fGK264;
1950 0 : double wtMax = 16. * s3 * s4 * flavWt
1951 0 : * ( (tHres*tHres + uHres*uHres + 2. * sH * (s3 + s4)) / (tHres * uHres)
1952 0 : - s3 * s4 * (1. / (tHres*tHres) + 1. / (uHres*uHres)) );
1953 :
1954 : // Done.
1955 0 : return wt / wtMax;
1956 :
1957 0 : }
1958 :
1959 : //==========================================================================
1960 :
1961 : // Sigma2ffbar2ZW class.
1962 : // Cross section for f fbar' -> Z0 W+- (f is quark or lepton).
1963 :
1964 : //--------------------------------------------------------------------------
1965 :
1966 : // Initialize process.
1967 :
1968 : void Sigma2ffbar2ZW::initProc() {
1969 :
1970 : // Store W+- mass and width for propagator.
1971 0 : mW = particleDataPtr->m0(24);
1972 0 : widW = particleDataPtr->mWidth(24);
1973 0 : mWS = mW*mW;
1974 0 : mwWS = pow2(mW * widW);
1975 :
1976 : // Left-handed couplings for up/nu- and down/e-type quarks.
1977 0 : lun = (hasLeptonBeams) ? couplingsPtr->lf(12) : couplingsPtr->lf(2);
1978 0 : lde = (hasLeptonBeams) ? couplingsPtr->lf(11) : couplingsPtr->lf(1);
1979 :
1980 : // Common weak coupling factor.
1981 0 : sin2thetaW = couplingsPtr->sin2thetaW();
1982 0 : cos2thetaW = couplingsPtr->cos2thetaW();
1983 0 : thetaWRat = 1. / (4. * cos2thetaW);
1984 0 : cotT = sqrt(cos2thetaW / sin2thetaW);
1985 0 : thetaWpt = (9. - 8. * sin2thetaW) / 4.;
1986 0 : thetaWmm = (8. * sin2thetaW - 6.) / 4.;
1987 :
1988 : // Secondary open width fractions.
1989 0 : openFracPos = particleDataPtr->resOpenFrac(23, 24);
1990 0 : openFracNeg = particleDataPtr->resOpenFrac(23, -24);
1991 :
1992 0 : }
1993 :
1994 : //--------------------------------------------------------------------------
1995 :
1996 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1997 :
1998 : void Sigma2ffbar2ZW::sigmaKin() {
1999 :
2000 : // Evaluate cross section, as programmed by Merlin Kole (after tidying),
2001 : // based on Brown, Sahdev, Mikaelian, Phys Rev. D20 (1979) 1069.
2002 : /*
2003 : double resBW = 1. / (pow2(sH - mWS) + mwWS);
2004 : double prefac = 12.0 * M_PI * pow2(alpEM) / (sH2 * 8. * sin2thetaW);
2005 : double temp1 = tH * uH - s3 * s4;
2006 : double temp2 = temp1 / (s3 * s4);
2007 : double temp3 = (s3 + s4) / sH;
2008 : double temp4 = s3 * s4 / sH;
2009 : double partA = temp2 * (0.25 - 0.5 * temp3
2010 : + (pow2(s3 + s4) + 8. * s3 * s4)/(4. * sH2) )
2011 : + (s3 + s4)/(s3 * s4) * (sH/2. - s3 - s4 + pow2(s3 - s4)/(2. * sH));
2012 : double partB1 = lun * (0.25 * temp2 * (1. - temp3 - 4. * temp4 / tH)
2013 : + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / tH) );
2014 : double partB2 = lde * ( 0.25 * temp2 * (1.- temp3 - 4. * temp4 / uH)
2015 : + ((s3 + s4)/(2. * s3 * s4)) * (sH - s3 - s4 + 2. * s3 * s4 / uH) );
2016 : double partE = 0.25 * temp2 + 0.5 *(s3 + s4) / temp4;
2017 : sigma0 = prefac * (pow2(cotT) * sH2 * resBW * partA
2018 : + 2.* sH * cotT * resBW * (sH - mWS) * (partB2 - partB1)
2019 : + pow2(lun - lde) * partE + pow2(lde) * temp1/uH2
2020 : + pow2(lun) * temp1/tH2 + 2. * lun * lde * sH * (s3 + s4) / (uH * tH));
2021 : */
2022 :
2023 : // Evaluate cross section. Expression from EHLQ, with bug fix,
2024 : // but can still give negative cross section so suspect.
2025 0 : double resBW = 1. / (pow2(sH - mWS) + mwWS);
2026 0 : sigma0 = (M_PI / sH2) * 0.5 * pow2(alpEM / sin2thetaW);
2027 0 : sigma0 *= sH * resBW * (thetaWpt * pT2 + thetaWmm * (s3 + s4))
2028 0 : + (sH - mWS) * resBW * sH * (pT2 - s3 - s4) * (lun / tH - lde / uH)
2029 0 : + thetaWRat * sH * pT2 * ( lun*lun / tH2 + lde*lde / uH2 )
2030 0 : + 2. * thetaWRat * sH * (s3 + s4) * lun * lde / (tH * uH);
2031 : // Need to protect against negative cross sections at times.
2032 0 : sigma0 = max(0., sigma0);
2033 :
2034 0 : }
2035 :
2036 : //--------------------------------------------------------------------------
2037 :
2038 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2039 :
2040 : double Sigma2ffbar2ZW::sigmaHat() {
2041 :
2042 : // CKM and colour factors.
2043 0 : double sigma = sigma0;
2044 0 : if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
2045 :
2046 : // Corrections for secondary widths in Z0 and W+- decays.
2047 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2048 0 : sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2049 :
2050 : // Answer.
2051 0 : return sigma;
2052 :
2053 : }
2054 :
2055 : //--------------------------------------------------------------------------
2056 :
2057 : // Select identity, colour and anticolour.
2058 :
2059 : void Sigma2ffbar2ZW::setIdColAcol() {
2060 :
2061 : // Sign of outgoing W.
2062 0 : int sign = 1 - 2 * (abs(id1)%2);
2063 0 : if (id1 < 0) sign = -sign;
2064 0 : setId( id1, id2, 23, 24 * sign);
2065 :
2066 : // tHat is defined between (f, W-) or (fbar, W+),
2067 : // so OK for u/ubar on side 1, but must swap tHat <-> uHat if d/dbar.
2068 0 : if (abs(id1)%2 == 1) swapTU = true;
2069 :
2070 : // Colour flow topologies. Swap when antiquarks.
2071 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2072 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2073 0 : if (id1 < 0) swapColAcol();
2074 :
2075 0 : }
2076 :
2077 : //--------------------------------------------------------------------------
2078 :
2079 : // Evaluate weight for Z0 and W+- decay angles.
2080 :
2081 : double Sigma2ffbar2ZW::weightDecay( Event& process, int iResBeg, int iResEnd) {
2082 :
2083 : // Two resonance decays, but with common weight.
2084 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
2085 :
2086 : // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6)
2087 : // with f' fbar' from W+- and f" fbar" from Z0 (note flip Z0 <-> W+-).
2088 0 : int i1 = (process[3].id() < 0) ? 3 : 4;
2089 0 : int i2 = 7 - i1;
2090 0 : int i3 = (process[9].id() > 0) ? 9 : 10;
2091 0 : int i4 = 19 - i3;
2092 0 : int i5 = (process[7].id() > 0) ? 7 : 8;
2093 0 : int i6 = 15 - i5;
2094 :
2095 : // Set up four-products and internal products.
2096 0 : setupProd( process, i1, i2, i3, i4, i5, i6);
2097 :
2098 : // Swap tHat and uHat if incoming fermion is downtype.
2099 0 : double tHres = tH;
2100 0 : double uHres = uH;
2101 0 : if (process[i2].id()%2 == 1) swap( tHres, uHres);
2102 :
2103 : // Couplings of incoming (anti)fermions and outgoing from Z0.
2104 0 : int idAbs = process[i1].idAbs();
2105 0 : double ai = couplingsPtr->af(idAbs);
2106 0 : double li1 = couplingsPtr->lf(idAbs);
2107 0 : idAbs = process[i2].idAbs();
2108 0 : double li2 = couplingsPtr->lf(idAbs);
2109 0 : idAbs = process[i5].idAbs();
2110 0 : double l4 = couplingsPtr->lf(idAbs);
2111 0 : double r4 = couplingsPtr->rf(idAbs);
2112 :
2113 : // W propagator/interference factor.
2114 0 : double Wint = cos2thetaW * (sH - mWS) / (pow2(sH - mWS) + mwWS);
2115 :
2116 : // Combinations of couplings and kinematics (norm(x) = |x|^2).
2117 0 : double aWZ = li2 / tHres - 2. * Wint * ai;
2118 0 : double bWZ = li1 / uHres + 2. * Wint * ai;
2119 0 : double fGK135 = norm( aWZ * fGK( 1, 2, 3, 4, 5, 6)
2120 0 : + bWZ * fGK( 1, 2, 5, 6, 3, 4) );
2121 0 : double fGK136 = norm( aWZ * fGK( 1, 2, 3, 4, 6, 5)
2122 0 : + bWZ * fGK( 1, 2, 6, 5, 3, 4) );
2123 0 : double xiT = xiGK( tHres, uHres);
2124 0 : double xiU = xiGK( uHres, tHres);
2125 0 : double xjTU = xjGK( tHres, uHres);
2126 :
2127 : // Weight and maximum weight.
2128 0 : double wt = l4*l4 * fGK135 + r4*r4 * fGK136;
2129 0 : double wtMax = 4. * s3 * s4 * (l4*l4 + r4*r4)
2130 0 : * (aWZ * aWZ * xiT + bWZ * bWZ * xiU + aWZ * bWZ * xjTU);
2131 :
2132 : // Done.
2133 0 : return wt / wtMax;
2134 :
2135 0 : }
2136 :
2137 : //==========================================================================
2138 :
2139 : // Sigma2ffbar2WW class.
2140 : // Cross section for f fbar -> W- W+ (f is quark or lepton).
2141 :
2142 : //--------------------------------------------------------------------------
2143 :
2144 : // Initialize process.
2145 :
2146 : void Sigma2ffbar2WW::initProc() {
2147 :
2148 : // Store Z0 mass and width for propagator. Common coupling factor.
2149 0 : mZ = particleDataPtr->m0(23);
2150 0 : widZ = particleDataPtr->mWidth(23);
2151 0 : mZS = mZ*mZ;
2152 0 : mwZS = pow2(mZ * widZ);
2153 0 : thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
2154 :
2155 : // Secondary open width fraction.
2156 0 : openFracPair = particleDataPtr->resOpenFrac(24, -24);
2157 :
2158 0 : }
2159 :
2160 : //--------------------------------------------------------------------------
2161 :
2162 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2163 :
2164 : void Sigma2ffbar2WW::sigmaKin() {
2165 :
2166 : // Cross section part common for all incoming flavours.
2167 0 : sigma0 = (M_PI / sH2) * pow2(alpEM);
2168 :
2169 : // Z0 propagator and gamma*/Z0 interference.
2170 0 : double Zprop = sH2 / (pow2(sH - mZS) + mwZS);
2171 0 : double Zint = Zprop * (1. - mZS / sH);
2172 :
2173 : // Common coupling factors (g = gamma*, Z = Z0, f = t-channel fermion).
2174 0 : cgg = 0.5;
2175 0 : cgZ = thetaWRat * Zint;
2176 0 : cZZ = 0.5 * pow2(thetaWRat) * Zprop;
2177 0 : cfg = thetaWRat;
2178 0 : cfZ = pow2(thetaWRat) * Zint;
2179 0 : cff = pow2(thetaWRat);
2180 :
2181 : // Kinematical functions.
2182 0 : double rat34 = sH * (2. * (s3 + s4) + pT2) / (s3 * s4);
2183 0 : double lambdaS = pow2(sH - s3 - s4) - 4. * s3 * s4;
2184 0 : double intA = (sH - s3 - s4) * rat34 / sH;
2185 0 : double intB = 4. * (s3 + s4 - pT2);
2186 0 : gSS = (lambdaS * rat34 + 12. * sH * pT2) / sH2;
2187 0 : gTT = rat34 + 4. * sH * pT2 / tH2;
2188 0 : gST = intA + intB / tH;
2189 0 : gUU = rat34 + 4. * sH * pT2 / uH2;
2190 0 : gSU = intA + intB / uH;
2191 :
2192 0 : }
2193 :
2194 : //--------------------------------------------------------------------------
2195 :
2196 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2197 :
2198 : double Sigma2ffbar2WW::sigmaHat() {
2199 :
2200 : // Flavour-specific couplings.
2201 0 : int idAbs = abs(id1);
2202 0 : double ei = couplingsPtr->ef(idAbs);
2203 0 : double vi = couplingsPtr->vf(idAbs);
2204 0 : double ai = couplingsPtr->af(idAbs);
2205 :
2206 : // Combine, with different cases for up- and down-type in-flavours.
2207 0 : double sigma = sigma0;
2208 0 : sigma *= (idAbs%2 == 1)
2209 0 : ? (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
2210 0 : + (cfg * ei + cfZ * (vi + ai)) * gST + cff * gTT
2211 : : (cgg * ei*ei + cgZ * ei * vi + cZZ * (vi*vi + ai*ai)) * gSS
2212 0 : - (cfg * ei + cfZ * (vi + ai)) * gSU + cff * gUU;
2213 :
2214 : // Initial-state colour factor. Correction for secondary widths. Answer.
2215 0 : if (idAbs < 9) sigma /= 3.;
2216 0 : sigma *= openFracPair;
2217 0 : return sigma;
2218 :
2219 : }
2220 :
2221 : //--------------------------------------------------------------------------
2222 :
2223 : // Select identity, colour and anticolour.
2224 :
2225 : void Sigma2ffbar2WW::setIdColAcol() {
2226 :
2227 : // Always order W- W+, i.e. W- first.
2228 0 : setId( id1, id2, -24, 24);
2229 :
2230 : // tHat is defined between (f, W-) or (fbar, W+),
2231 0 : if (id1 < 0) swapTU = true;
2232 :
2233 : // Colour flow topologies. Swap when antiquarks.
2234 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2235 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2236 0 : if (id1 < 0) swapColAcol();
2237 :
2238 0 : }
2239 :
2240 : //--------------------------------------------------------------------------
2241 :
2242 : // Evaluate weight for W+ and W- decay angles.
2243 :
2244 : double Sigma2ffbar2WW::weightDecay( Event& process, int iResBeg, int iResEnd) {
2245 :
2246 : // Two resonance decays, but with common weight.
2247 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
2248 :
2249 : // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
2250 : // with f' fbar' from W- and f" fbar" from W+.
2251 0 : int i1 = (process[3].id() < 0) ? 3 : 4;
2252 0 : int i2 = 7 - i1;
2253 0 : int i3 = (process[7].id() > 0) ? 7 : 8;
2254 0 : int i4 = 15 - i3;
2255 0 : int i5 = (process[9].id() > 0) ? 9 : 10;
2256 0 : int i6 = 19 - i5;
2257 :
2258 : // Set up four-products and internal products.
2259 0 : setupProd( process, i1, i2, i3, i4, i5, i6);
2260 :
2261 : // tHat and uHat of fbar f -> W- W+ opposite to previous convention.
2262 0 : double tHres = uH;
2263 0 : double uHres = tH;
2264 :
2265 : // Couplings of incoming (anti)fermion.
2266 0 : int idAbs = process[i1].idAbs();
2267 0 : double ai = couplingsPtr->af(idAbs);
2268 0 : double li = couplingsPtr->lf(idAbs);
2269 0 : double ri = couplingsPtr->rf(idAbs);
2270 :
2271 : // gamma*/Z0 propagator/interference factor.
2272 0 : double Zint = mZS * (sH - mZS) / (pow2(sH - mZS) + mwZS);
2273 :
2274 : // Combinations of couplings and kinematics (norm(x) = |x|^2).
2275 0 : double dWW = (li * Zint + ai) / sH;
2276 0 : double aWW = dWW + 0.5 * (ai + 1.) / tHres;
2277 0 : double bWW = dWW + 0.5 * (ai - 1.) / uHres;
2278 0 : double cWW = ri * Zint / sH;
2279 0 : double fGK135 = norm( aWW * fGK( 1, 2, 3, 4, 5, 6)
2280 0 : - bWW * fGK( 1, 2, 5, 6, 3, 4) );
2281 0 : double fGK253 = norm( cWW * ( fGK( 2, 1, 5, 6, 3, 4)
2282 0 : - fGK( 2, 1, 3, 4, 5, 6) ) );
2283 0 : double xiT = xiGK( tHres, uHres);
2284 0 : double xiU = xiGK( uHres, tHres);
2285 0 : double xjTU = xjGK( tHres, uHres);
2286 :
2287 : // Weight and maximum weight.
2288 0 : double wt = fGK135 + fGK253;
2289 0 : double wtMax = 4. * s3 * s4
2290 0 : * ( aWW * aWW * xiT + bWW * bWW * xiU - aWW * bWW * xjTU
2291 0 : + cWW * cWW * (xiT + xiU - xjTU) );
2292 :
2293 : // Done.
2294 0 : return wt / wtMax;
2295 0 : }
2296 :
2297 : //==========================================================================
2298 :
2299 : // Sigma2ffbargmZggm class.
2300 : // Collects common methods for f fbar -> gamma*/Z0 g/gamma and permutations.
2301 :
2302 : //--------------------------------------------------------------------------
2303 :
2304 : // Initialize process.
2305 :
2306 : void Sigma2ffbargmZggm::initProc() {
2307 :
2308 : // Allow to pick only gamma* or Z0 part of full gamma*/Z0 expression.
2309 0 : gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
2310 :
2311 : // Store Z0 mass and width for propagator.
2312 0 : mRes = particleDataPtr->m0(23);
2313 0 : GammaRes = particleDataPtr->mWidth(23);
2314 0 : m2Res = mRes*mRes;
2315 0 : GamMRat = GammaRes / mRes;
2316 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
2317 0 : * couplingsPtr->cos2thetaW());
2318 :
2319 : // Set pointer to particle properties and decay table.
2320 0 : particlePtr = particleDataPtr->particleDataEntryPtr(23);
2321 :
2322 0 : }
2323 :
2324 : //--------------------------------------------------------------------------
2325 :
2326 : // Evaluate sum of flavour couplings times phase space.
2327 :
2328 : void Sigma2ffbargmZggm::flavSum() {
2329 :
2330 : // Coupling factors for Z0 subsystem.
2331 0 : double alpSZ = couplingsPtr->alphaS(s3);
2332 0 : double colQZ = 3. * (1. + alpSZ / M_PI);
2333 :
2334 : // Reset quantities to sum. Declare variables in loop.
2335 0 : gamSum = 0.;
2336 0 : intSum = 0.;
2337 0 : resSum = 0.;
2338 : int onMode;
2339 0 : double mf, mr, psvec, psaxi, betaf, ef2, efvf, vf2af2, colf;
2340 :
2341 : // Loop over all Z0 decay channels.
2342 0 : for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
2343 0 : int idAbs = abs( particlePtr->channel(i).product(0) );
2344 :
2345 : // Only contributions from three fermion generations, except top.
2346 0 : if ( (idAbs > 0 && idAbs < 6) || ( idAbs > 10 && idAbs < 17)) {
2347 0 : mf = particleDataPtr->m0(idAbs);
2348 :
2349 : // Check that above threshold. Phase space.
2350 0 : if (m3 > 2. * mf + MASSMARGIN) {
2351 0 : mr = pow2(mf / m3);
2352 0 : betaf = sqrtpos(1. - 4. * mr);
2353 0 : psvec = betaf * (1. + 2. * mr);
2354 0 : psaxi = pow3(betaf);
2355 :
2356 : // Combine phase space with couplings.
2357 0 : ef2 = couplingsPtr->ef2(idAbs) * psvec;
2358 0 : efvf = couplingsPtr->efvf(idAbs) * psvec;
2359 0 : vf2af2 = couplingsPtr->vf2(idAbs) * psvec
2360 0 : + couplingsPtr->af2(idAbs) * psaxi;
2361 0 : colf = (idAbs < 6) ? colQZ : 1.;
2362 :
2363 : // Store sum of combinations. For outstate only open channels.
2364 0 : onMode = particlePtr->channel(i).onMode();
2365 0 : if (onMode == 1 || onMode == 2) {
2366 0 : gamSum += colf * ef2;
2367 0 : intSum += colf * efvf;
2368 0 : resSum += colf * vf2af2;
2369 0 : }
2370 :
2371 : // End loop over fermions.
2372 : }
2373 : }
2374 : }
2375 :
2376 : // Done. Return values in gamSum, intSum and resSum.
2377 :
2378 0 : }
2379 :
2380 : //--------------------------------------------------------------------------
2381 :
2382 : // Calculate common parts of gamma/interference/Z0 propagator terms.
2383 :
2384 : void Sigma2ffbargmZggm::propTerm() {
2385 :
2386 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
2387 0 : gamProp = 4. * alpEM / (3. * M_PI * s3);
2388 0 : intProp = gamProp * 2. * thetaWRat * s3 * (s3 - m2Res)
2389 0 : / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2390 0 : resProp = gamProp * pow2(thetaWRat * s3)
2391 0 : / ( pow2(s3 - m2Res) + pow2(s3 * GamMRat) );
2392 :
2393 : // Optionally only keep gamma* or Z0 term.
2394 0 : if (gmZmode == 1) {intProp = 0.; resProp = 0.;}
2395 0 : if (gmZmode == 2) {gamProp = 0.; intProp = 0.;}
2396 :
2397 0 : }
2398 :
2399 : //--------------------------------------------------------------------------
2400 :
2401 : // Evaluate weight for gamma*/Z0 decay angle.
2402 :
2403 : double Sigma2ffbargmZggm::weightDecay( Event& process, int iResBeg,
2404 : int iResEnd) {
2405 :
2406 : // Z should sit in entry 5 and one more parton in entry 6.
2407 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
2408 :
2409 : // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2410 : // where f' fbar' come from gamma*/Z0 decay.
2411 : int i1, i2;
2412 0 : int i3 = (process[7].id() > 0) ? 7 : 8;
2413 0 : int i4 = 15 - i3;
2414 :
2415 : // Order so that fbar(1) f(2) -> gamma*/Z0 g/gamma.
2416 0 : if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2417 0 : i1 = (process[3].id() < 0) ? 3 : 4;
2418 0 : i2 = 7 - i1;
2419 :
2420 : // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) gamma*/Z0.
2421 0 : } else if (process[3].idAbs() < 20) {
2422 0 : i1 = (process[3].id() < 0) ? 3 : 6;
2423 0 : i2 = 9 - i1;
2424 0 : } else {
2425 0 : i1 = (process[4].id() < 0) ? 4 : 6;
2426 0 : i2 = 10 - i1;
2427 : }
2428 :
2429 : // Charge/2, left- and righthanded couplings for in- and out-fermion.
2430 0 : int id1Abs = process[i1].idAbs();
2431 0 : double ei = 0.5 * couplingsPtr->ef(id1Abs);
2432 0 : double li = couplingsPtr->lf(id1Abs);
2433 0 : double ri = couplingsPtr->rf(id1Abs);
2434 0 : int id3Abs = process[i3].idAbs();
2435 0 : double ef = 0.5 * couplingsPtr->ef(id3Abs);
2436 0 : double lf = couplingsPtr->lf(id3Abs);
2437 0 : double rf = couplingsPtr->rf(id3Abs);
2438 :
2439 : // Combinations of left/right for in/out, gamma*/interference/Z0.
2440 0 : double clilf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*lf
2441 0 : + li*li * resProp * lf*lf;
2442 0 : double clirf = ei*ei * gamProp * ef*ef + ei*li * intProp * ef*rf
2443 0 : + li*li * resProp * rf*rf;
2444 0 : double crilf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*lf
2445 0 : + ri*ri * resProp * lf*lf;
2446 0 : double crirf = ei*ei * gamProp * ef*ef + ei*ri * intProp * ef*rf
2447 0 : + ri*ri * resProp * rf*rf;
2448 :
2449 : // Evaluate four-vector products.
2450 0 : double p13 = process[i1].p() * process[i3].p();
2451 0 : double p14 = process[i1].p() * process[i4].p();
2452 0 : double p23 = process[i2].p() * process[i3].p();
2453 0 : double p24 = process[i2].p() * process[i4].p();
2454 :
2455 : // Calculate weight and its maximum.
2456 0 : double wt = (clilf + crirf) * (p13*p13 + p24*p24)
2457 0 : + (clirf + crilf) * (p14*p14 + p23*p23) ;
2458 0 : double wtMax = (clilf + clirf + crilf + crirf)
2459 0 : * (pow2(p13 + p14) + pow2(p23 + p24));
2460 :
2461 : // Done.
2462 0 : return (wt / wtMax);
2463 :
2464 0 : }
2465 :
2466 : //==========================================================================
2467 :
2468 : // Sigma2qqbar2gmZg class.
2469 : // Cross section for q qbar -> gamma*/Z0 g.
2470 :
2471 : //--------------------------------------------------------------------------
2472 :
2473 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2474 :
2475 : void Sigma2qqbar2gmZg::sigmaKin() {
2476 :
2477 : // Cross section part common for all incoming flavours.
2478 0 : sigma0 = (M_PI / sH2) * (alpEM * alpS)
2479 0 : * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2480 :
2481 : // Calculate flavour sums for final state.
2482 0 : flavSum();
2483 :
2484 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
2485 0 : propTerm();
2486 :
2487 0 : }
2488 :
2489 : //--------------------------------------------------------------------------
2490 :
2491 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2492 :
2493 : double Sigma2qqbar2gmZg::sigmaHat() {
2494 :
2495 : // Combine gamma, interference and Z0 parts.
2496 0 : int idAbs = abs(id1);
2497 0 : double sigma = sigma0
2498 0 : * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2499 0 : + couplingsPtr->efvf(idAbs) * intProp * intSum
2500 0 : + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2501 :
2502 : // Correct for the running-width Z0 propagater weight in PhaseSpace.
2503 0 : sigma /= runBW3;
2504 :
2505 : // Answer.
2506 0 : return sigma;
2507 :
2508 : }
2509 :
2510 : //--------------------------------------------------------------------------
2511 :
2512 : // Select identity, colour and anticolour.
2513 :
2514 : void Sigma2qqbar2gmZg::setIdColAcol() {
2515 :
2516 : // Flavours trivial.
2517 0 : setId( id1, id2, 23, 21);
2518 :
2519 : // Colour flow topologies. Swap when antiquarks.
2520 0 : setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2521 0 : if (id1 < 0) swapColAcol();
2522 :
2523 0 : }
2524 :
2525 : //==========================================================================
2526 :
2527 : // Sigma2qg2gmZq class.
2528 : // Cross section for q g -> gamma*/Z0 q.
2529 :
2530 : //--------------------------------------------------------------------------
2531 :
2532 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2533 :
2534 : void Sigma2qg2gmZq::sigmaKin() {
2535 :
2536 : // Cross section part common for all incoming flavours.
2537 0 : sigma0 = (M_PI / sH2) * (alpEM * alpS)
2538 0 : * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2539 :
2540 : // Calculate flavour sums for final state.
2541 0 : flavSum();
2542 :
2543 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
2544 0 : propTerm();
2545 :
2546 0 : }
2547 :
2548 : //--------------------------------------------------------------------------
2549 :
2550 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2551 :
2552 : double Sigma2qg2gmZq::sigmaHat() {
2553 :
2554 : // Combine gamma, interference and Z0 parts.
2555 0 : int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2556 0 : double sigma = sigma0
2557 0 : * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2558 0 : + couplingsPtr->efvf(idAbs) * intProp * intSum
2559 0 : + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2560 :
2561 : // Correct for the running-width Z0 propagater weight in PhaseSpace.
2562 0 : sigma /= runBW3;
2563 :
2564 : // Answer.
2565 0 : return sigma;
2566 :
2567 : }
2568 :
2569 : //--------------------------------------------------------------------------
2570 :
2571 : // Select identity, colour and anticolour.
2572 :
2573 : void Sigma2qg2gmZq::setIdColAcol() {
2574 :
2575 : // Flavour set up for q g -> gamma*/Z0 q.
2576 0 : int idq = (id2 == 21) ? id1 : id2;
2577 0 : setId( id1, id2, 23, idq);
2578 :
2579 : // tH defined between f and f': must swap tHat <-> uHat if q g in.
2580 0 : swapTU = (id2 == 21);
2581 :
2582 : // Colour flow topologies. Swap when antiquarks.
2583 0 : if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2584 0 : else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2585 0 : if (idq < 0) swapColAcol();
2586 :
2587 0 : }
2588 :
2589 : //==========================================================================
2590 :
2591 : // Sigma2ffbar2gmZgm class.
2592 : // Cross section for f fbar -> gamma*/Z0 gamma.
2593 :
2594 : //--------------------------------------------------------------------------
2595 :
2596 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2597 :
2598 : void Sigma2ffbar2gmZgm::sigmaKin() {
2599 :
2600 : // Cross section part common for all incoming flavours.
2601 0 : sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2602 0 : * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2603 :
2604 : // Calculate flavour sums for final state.
2605 0 : flavSum();
2606 :
2607 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
2608 0 : propTerm();
2609 :
2610 :
2611 0 : }
2612 :
2613 : //--------------------------------------------------------------------------
2614 :
2615 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2616 :
2617 : double Sigma2ffbar2gmZgm::sigmaHat() {
2618 :
2619 : // Combine gamma, interference and Z0 parts.
2620 0 : int idAbs = abs(id1);
2621 0 : double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2622 0 : * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2623 0 : + couplingsPtr->efvf(idAbs) * intProp * intSum
2624 0 : + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2625 :
2626 : // Correct for the running-width Z0 propagater weight in PhaseSpace.
2627 0 : sigma /= runBW3;
2628 :
2629 : // Colour factor. Answer.
2630 0 : if (idAbs < 9) sigma /= 3.;
2631 0 : return sigma;
2632 :
2633 : }
2634 :
2635 : //--------------------------------------------------------------------------
2636 :
2637 : // Select identity, colour and anticolour.
2638 :
2639 : void Sigma2ffbar2gmZgm::setIdColAcol() {
2640 :
2641 : // Flavours trivial.
2642 0 : setId( id1, id2, 23, 22);
2643 :
2644 : // Colour flow topologies. Swap when antiquarks.
2645 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2646 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2647 0 : if (id1 < 0) swapColAcol();
2648 :
2649 0 : }
2650 :
2651 : //==========================================================================
2652 :
2653 : // Sigma2fgm2gmZf class.
2654 : // Cross section for f gamma -> gamma*/Z0 f'.
2655 :
2656 : //--------------------------------------------------------------------------
2657 :
2658 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2659 :
2660 : void Sigma2fgm2gmZf::sigmaKin() {
2661 :
2662 : // Cross section part common for all incoming flavours.
2663 0 : sigma0 = (M_PI / sH2) * (alpEM*alpEM)
2664 0 : * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (- sH * uH);
2665 :
2666 : // Calculate flavour sums for final state.
2667 0 : flavSum();
2668 :
2669 : // Calculate prefactors for gamma/interference/Z0 cross section terms.
2670 0 : propTerm();
2671 :
2672 0 : }
2673 :
2674 : //--------------------------------------------------------------------------
2675 :
2676 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2677 :
2678 : double Sigma2fgm2gmZf::sigmaHat() {
2679 :
2680 : // Combine gamma, interference and Z0 parts.
2681 0 : int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
2682 0 : double sigma = sigma0 * couplingsPtr->ef2(idAbs)
2683 0 : * ( couplingsPtr->ef2(idAbs) * gamProp * gamSum
2684 0 : + couplingsPtr->efvf(idAbs) * intProp * intSum
2685 0 : + couplingsPtr->vf2af2(idAbs) * resProp * resSum);
2686 :
2687 : // Correct for the running-width Z0 propagater weight in PhaseSpace.
2688 0 : sigma /= runBW3;
2689 :
2690 : // Answer.
2691 0 : return sigma;
2692 :
2693 : }
2694 :
2695 : //--------------------------------------------------------------------------
2696 :
2697 : // Select identity, colour and anticolour.
2698 :
2699 : void Sigma2fgm2gmZf::setIdColAcol() {
2700 :
2701 : // Flavour set up for q gamma -> gamma*/Z0 q.
2702 0 : int idq = (id2 == 22) ? id1 : id2;
2703 0 : setId( id1, id2, 23, idq);
2704 :
2705 : // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
2706 0 : swapTU = (id2 == 22);
2707 :
2708 : // Colour flow topologies. Swap when antiquarks.
2709 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
2710 0 : else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
2711 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2712 0 : if (idq < 0) swapColAcol();
2713 :
2714 0 : }
2715 :
2716 : //==========================================================================
2717 :
2718 : // Sigma2ffbarWggm class.
2719 : // Collects common methods for f fbar -> W+- g/gamma and permutations.
2720 :
2721 : //--------------------------------------------------------------------------
2722 :
2723 : // Evaluate weight for W+- decay angle.
2724 :
2725 : double Sigma2ffbarWggm::weightDecay( Event& process, int iResBeg,
2726 : int iResEnd) {
2727 :
2728 : // W should sit in entry 5 and one more parton in entry 6.
2729 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
2730 :
2731 : // In an outgoing sense fermions are labelled f(1) fbar(2) f'(3) fbar'(4)
2732 : // where f' fbar' come from W+- decay.
2733 : int i1, i2;
2734 0 : int i3 = (process[7].id() > 0) ? 7 : 8;
2735 0 : int i4 = 15 - i3;
2736 :
2737 : // Order so that fbar(1) f(2) -> W+- g/gamma.
2738 0 : if (process[3].idAbs() < 20 && process[4].idAbs() < 20) {
2739 0 : i1 = (process[3].id() < 0) ? 3 : 4;
2740 0 : i2 = 7 - i1;
2741 :
2742 : // Order so that f(2)/fbar(1) g/gamma -> f(1)/fbar(2) f'(3) W+-.
2743 0 : } else if (process[3].idAbs() < 20) {
2744 0 : i1 = (process[3].id() < 0) ? 3 : 6;
2745 0 : i2 = 9 - i1;
2746 0 : } else {
2747 0 : i1 = (process[4].id() < 0) ? 4 : 6;
2748 0 : i2 = 10 - i1;
2749 : }
2750 :
2751 : // Evaluate four-vector products.
2752 0 : double p13 = process[i1].p() * process[i3].p();
2753 0 : double p14 = process[i1].p() * process[i4].p();
2754 0 : double p23 = process[i2].p() * process[i3].p();
2755 0 : double p24 = process[i2].p() * process[i4].p();
2756 :
2757 : // Calculate weight and its maximum.
2758 0 : double wt = pow2(p13) + pow2(p24);
2759 0 : double wtMax = pow2(p13 + p14) + pow2(p23 + p24);
2760 :
2761 : // Done.
2762 0 : return (wt / wtMax);
2763 :
2764 0 : }
2765 :
2766 : //==========================================================================
2767 :
2768 : // Sigma2qqbar2Wg class.
2769 : // Cross section for q qbar' -> W+- g.
2770 :
2771 : //--------------------------------------------------------------------------
2772 :
2773 : // Initialize process.
2774 :
2775 : void Sigma2qqbar2Wg::initProc() {
2776 :
2777 : // Secondary open width fractions, relevant for top (or heavier).
2778 0 : openFracPos = particleDataPtr->resOpenFrac(24);
2779 0 : openFracNeg = particleDataPtr->resOpenFrac(-24);
2780 :
2781 0 : }
2782 :
2783 : //--------------------------------------------------------------------------
2784 :
2785 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2786 :
2787 : void Sigma2qqbar2Wg::sigmaKin() {
2788 :
2789 : // Cross section part common for all incoming flavours.
2790 0 : sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2791 0 : * (2./9.) * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2792 :
2793 0 : }
2794 :
2795 : //--------------------------------------------------------------------------
2796 :
2797 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2798 :
2799 : double Sigma2qqbar2Wg::sigmaHat() {
2800 :
2801 : // CKM factor. Secondary width for W+ or W-.
2802 0 : double sigma = sigma0 * couplingsPtr->V2CKMid(abs(id1), abs(id2));
2803 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2804 0 : sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2805 :
2806 : // Answer.
2807 0 : return sigma;
2808 :
2809 : }
2810 :
2811 : //--------------------------------------------------------------------------
2812 :
2813 : // Select identity, colour and anticolour.
2814 :
2815 : void Sigma2qqbar2Wg::setIdColAcol() {
2816 :
2817 : // Sign of outgoing W.
2818 0 : int sign = 1 - 2 * (abs(id1)%2);
2819 0 : if (id1 < 0) sign = -sign;
2820 0 : setId( id1, id2, 24 * sign, 21);
2821 :
2822 : // Colour flow topologies. Swap when antiquarks.
2823 0 : setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
2824 0 : if (id1 < 0) swapColAcol();
2825 :
2826 0 : }
2827 :
2828 : //==========================================================================
2829 :
2830 : // Sigma2qg2Wq class.
2831 : // Cross section for q g -> W+- q'.
2832 :
2833 : //--------------------------------------------------------------------------
2834 :
2835 : // Initialize process.
2836 :
2837 : void Sigma2qg2Wq::initProc() {
2838 :
2839 : // Secondary open width fractions, relevant for top (or heavier).
2840 0 : openFracPos = particleDataPtr->resOpenFrac(24);
2841 0 : openFracNeg = particleDataPtr->resOpenFrac(-24);
2842 :
2843 0 : }
2844 :
2845 : //--------------------------------------------------------------------------
2846 :
2847 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2848 :
2849 : void Sigma2qg2Wq::sigmaKin() {
2850 :
2851 : // Cross section part common for all incoming flavours.
2852 0 : sigma0 = (M_PI / sH2) * (alpEM * alpS / couplingsPtr->sin2thetaW())
2853 0 : * (1./12.) * (sH2 + uH2 + 2. * tH * s3) / (-sH * uH);
2854 :
2855 0 : }
2856 :
2857 : //--------------------------------------------------------------------------
2858 :
2859 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2860 :
2861 : double Sigma2qg2Wq::sigmaHat() {
2862 :
2863 : // CKM factor. Secondary width for W+ or W-.
2864 0 : int idAbs = (id2 == 21) ? abs(id1) : abs(id2);
2865 0 : double sigma = sigma0 * couplingsPtr->V2CKMsum(idAbs);
2866 0 : int idUp = (id2 == 21) ? id1 : id2;
2867 0 : if (idAbs%2 == 1) idUp = -idUp;
2868 0 : sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2869 :
2870 : // Answer.
2871 0 : return sigma;
2872 :
2873 : }
2874 :
2875 : //--------------------------------------------------------------------------
2876 :
2877 : // Select identity, colour and anticolour.
2878 :
2879 : void Sigma2qg2Wq::setIdColAcol() {
2880 :
2881 : // Sign of outgoing W. Flavour of outgoing quark.
2882 0 : int idq = (id2 == 21) ? id1 : id2;
2883 0 : int sign = 1 - 2 * (abs(idq)%2);
2884 0 : if (idq < 0) sign = -sign;
2885 0 : id4 = couplingsPtr->V2CKMpick(idq);
2886 :
2887 : // Flavour set up for q g -> W q.
2888 0 : setId( id1, id2, 24 * sign, id4);
2889 :
2890 : // tH defined between f and f': must swap tHat <-> uHat if q g in.
2891 0 : swapTU = (id2 == 21);
2892 :
2893 : // Colour flow topologies. Swap when antiquarks.
2894 0 : if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2895 0 : else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2896 0 : if (idq < 0) swapColAcol();
2897 :
2898 0 : }
2899 :
2900 : //==========================================================================
2901 :
2902 : // Sigma2ffbar2Wgm class.
2903 : // Cross section for f fbar' -> W+- gamma.
2904 :
2905 : //--------------------------------------------------------------------------
2906 :
2907 : // Initialize process.
2908 :
2909 : void Sigma2ffbar2Wgm::initProc() {
2910 :
2911 : // Secondary open width fractions, relevant for top (or heavier).
2912 0 : openFracPos = particleDataPtr->resOpenFrac(24);
2913 0 : openFracNeg = particleDataPtr->resOpenFrac(-24);
2914 :
2915 0 : }
2916 :
2917 : //--------------------------------------------------------------------------
2918 :
2919 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2920 :
2921 : void Sigma2ffbar2Wgm::sigmaKin() {
2922 :
2923 : // Cross section part common for all incoming flavours.
2924 0 : sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2925 0 : * 0.5 * (tH2 + uH2 + 2. * sH * s3) / (tH * uH);
2926 0 : }
2927 :
2928 : //--------------------------------------------------------------------------
2929 :
2930 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2931 :
2932 : double Sigma2ffbar2Wgm::sigmaHat() {
2933 :
2934 : // Extrafactor different for e nu and q qbar' instate.
2935 0 : int id1Abs = abs(id1);
2936 0 : int id2Abs = abs(id2);
2937 0 : double chgUp = (id1Abs > 10) ? 0. : 2./3.;
2938 0 : double sigma = sigma0 * pow2( chgUp - tH / (tH + uH) );
2939 :
2940 : // CKM and colour factors. Secondary width for W+ or W-.
2941 0 : if (id1Abs < 9) sigma *= couplingsPtr->V2CKMid(id1Abs, id2Abs) / 3.;
2942 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2943 0 : sigma *= (idUp > 0) ? openFracPos : openFracNeg;
2944 :
2945 : // Answer.
2946 0 : return sigma;
2947 :
2948 : }
2949 :
2950 : //--------------------------------------------------------------------------
2951 :
2952 : // Select identity, colour and anticolour.
2953 :
2954 : void Sigma2ffbar2Wgm::setIdColAcol() {
2955 :
2956 : // Sign of outgoing W.
2957 0 : int sign = 1 - 2 * (abs(id1)%2);
2958 0 : if (id1 < 0) sign = -sign;
2959 0 : setId( id1, id2, 24 * sign, 22);
2960 :
2961 : // tH defined between (f,W-) or (fbar',W+).
2962 0 : swapTU = (sign * id1 > 0);
2963 :
2964 : // Colour flow topologies. Swap when antiquarks.
2965 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2966 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
2967 0 : if (id1 < 0) swapColAcol();
2968 :
2969 0 : }
2970 :
2971 : //==========================================================================
2972 :
2973 : // Sigma2fgm2Wf class.
2974 : // Cross section for f gamma -> W+- f'.
2975 :
2976 : //--------------------------------------------------------------------------
2977 :
2978 : // Initialize process.
2979 :
2980 : void Sigma2fgm2Wf::initProc() {
2981 :
2982 : // Secondary open width fractions, relevant for top (or heavier).
2983 0 : openFracPos = particleDataPtr->resOpenFrac(24);
2984 0 : openFracNeg = particleDataPtr->resOpenFrac(-24);
2985 :
2986 0 : }
2987 :
2988 : //--------------------------------------------------------------------------
2989 :
2990 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2991 :
2992 : void Sigma2fgm2Wf::sigmaKin() {
2993 :
2994 : // Cross section part common for all incoming flavours.
2995 0 : sigma0 = (M_PI / sH2) * (alpEM*alpEM / couplingsPtr->sin2thetaW())
2996 0 : * 0.5 * (sH2 + uH2 + 2. * tH * s3) / (pT2 * s3 - sH * uH);
2997 :
2998 0 : }
2999 :
3000 : //--------------------------------------------------------------------------
3001 :
3002 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
3003 :
3004 : double Sigma2fgm2Wf::sigmaHat() {
3005 :
3006 : // Extrafactor dependent on charge of incoming fermion.
3007 0 : int idAbs = (id2 == 22) ? abs(id1) : abs(id2);
3008 0 : double charge = (idAbs > 10) ? 1. : ( (idAbs%2 == 1) ? 1./3. : 2./3. );
3009 0 : double sigma = sigma0 * pow2( charge - sH / (sH + uH) );
3010 :
3011 : // CKM factor. Secondary width for W+ or W-.
3012 0 : sigma *= couplingsPtr->V2CKMsum(idAbs);
3013 0 : int idUp = (id2 == 22) ? id1 : id2;
3014 0 : if (idAbs%2 == 1) idUp = -idUp;
3015 0 : sigma *= (idUp > 0) ? openFracPos : openFracNeg;
3016 :
3017 : // Answer.
3018 0 : return sigma;
3019 :
3020 : }
3021 :
3022 : //--------------------------------------------------------------------------
3023 :
3024 : // Select identity, colour and anticolour.
3025 :
3026 : void Sigma2fgm2Wf::setIdColAcol() {
3027 :
3028 : // Sign of outgoing W. Flavour of outgoing fermion.
3029 0 : int idq = (id2 == 22) ? id1 : id2;
3030 0 : int sign = 1 - 2 * (abs(idq)%2);
3031 0 : if (idq < 0) sign = -sign;
3032 0 : id4 = couplingsPtr->V2CKMpick(idq);
3033 :
3034 : // Flavour set up for q gamma -> W q.
3035 0 : setId( id1, id2, 24 * sign, id4);
3036 :
3037 : // tH defined between f and f': must swap tHat <-> uHat if q gamma in.
3038 0 : swapTU = (id2 == 22);
3039 :
3040 : // Colour flow topologies. Swap when antiquarks.
3041 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0);
3042 0 : else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 1, 0);
3043 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
3044 0 : if (idq < 0) swapColAcol();
3045 :
3046 0 : }
3047 :
3048 : //==========================================================================
3049 :
3050 : // Sigma2gmgm2ffbar class.
3051 : // Cross section for gamma gamma -> l lbar.
3052 :
3053 : //--------------------------------------------------------------------------
3054 :
3055 : // Initialize process wrt flavour.
3056 :
3057 : void Sigma2gmgm2ffbar::initProc() {
3058 :
3059 : // Process name.
3060 0 : nameSave = "gamma gamma -> f fbar";
3061 0 : if (idNew == 1) nameSave = "gamma gamma -> q qbar (uds)";
3062 0 : if (idNew == 4) nameSave = "gamma gamma -> c cbar";
3063 0 : if (idNew == 5) nameSave = "gamma gamma -> b bbar";
3064 0 : if (idNew == 6) nameSave = "gamma gamma -> t tbar";
3065 0 : if (idNew == 11) nameSave = "gamma gamma -> e+ e-";
3066 0 : if (idNew == 13) nameSave = "gamma gamma -> mu+ mu-";
3067 0 : if (idNew == 15) nameSave = "gamma gamma -> tau+ tau-";
3068 :
3069 : // Generate massive phase space, except for u+d+s.
3070 0 : idMass = 0;
3071 0 : if (idNew > 3) idMass = idNew;
3072 :
3073 : // Charge and colour factor.
3074 0 : ef4 = 1.;
3075 0 : if (idNew == 1) ef4 = 3. * (pow4(2./3.) + 2. * pow4(1./3.));
3076 0 : if (idNew == 4 || idNew == 6) ef4 = 3. * pow4(2./3.);
3077 0 : if (idNew == 5) ef4 = 3. * pow4(1./3.);
3078 :
3079 : // Secondary open width fraction.
3080 0 : openFracPair = particleDataPtr->resOpenFrac(idNew, -idNew);
3081 :
3082 0 : }
3083 :
3084 : //--------------------------------------------------------------------------
3085 :
3086 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
3087 :
3088 : void Sigma2gmgm2ffbar::sigmaKin() {
3089 :
3090 : // Pick current flavour for u+d+s mix by e_q^4 weights.
3091 0 : if (idNew == 1) {
3092 0 : double rId = 18. * rndmPtr->flat();
3093 0 : idNow = 1;
3094 0 : if (rId > 1.) idNow = 2;
3095 0 : if (rId > 17.) idNow = 3;
3096 0 : s34Avg = pow2(particleDataPtr->m0(idNow));
3097 0 : } else {
3098 0 : idNow = idNew;
3099 0 : s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
3100 : }
3101 :
3102 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
3103 0 : double tHQ = -0.5 * (sH - tH + uH);
3104 0 : double uHQ = -0.5 * (sH + tH - uH);
3105 0 : double tHQ2 = tHQ * tHQ;
3106 0 : double uHQ2 = uHQ * uHQ;
3107 :
3108 : // Calculate kinematics dependence.
3109 0 : if (sH < 4. * s34Avg) sigTU = 0.;
3110 0 : else sigTU = 2. * (tHQ * uHQ - s34Avg * sH)
3111 0 : * (tHQ2 + uHQ2 + 2. * s34Avg * sH) / (tHQ2 * uHQ2);
3112 :
3113 : // Answer.
3114 0 : sigma = (M_PI / sH2) * pow2(alpEM) * ef4 * sigTU * openFracPair;
3115 :
3116 0 : }
3117 :
3118 : //--------------------------------------------------------------------------
3119 :
3120 : // Select identity, colour and anticolour.
3121 :
3122 : void Sigma2gmgm2ffbar::setIdColAcol() {
3123 :
3124 : // Flavours are trivial.
3125 0 : setId( id1, id2, idNow, -idNow);
3126 :
3127 : // Colour flow in singlet state.
3128 0 : if (idNow < 10) setColAcol( 0, 0, 0, 0, 1, 0, 0, 1);
3129 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
3130 :
3131 0 : }
3132 :
3133 : //==========================================================================
3134 :
3135 : } // end namespace Pythia8
|
