Line data Source code
1 : // SigmaHiggs.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // Part of code written by Marc Montull, CERN summer student 2007.
4 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 : // Function definitions (not found in the header) for the
7 : // Higgs simulation classes.
8 :
9 : #include "Pythia8/SigmaHiggs.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // Sigma1ffbar2H class.
16 : // Cross section for f fbar -> H0 , H1, H2 or A3.
17 : // (f is quark or lepton, H0 SM Higgs and H1, H2, A3 BSM Higgses ).
18 :
19 : //--------------------------------------------------------------------------
20 :
21 : // Initialize process.
22 :
23 : void Sigma1ffbar2H::initProc() {
24 :
25 : // Properties specific to Higgs state.
26 0 : if (higgsType == 0) {
27 0 : nameSave = "f fbar -> H (SM)";
28 0 : codeSave = 901;
29 0 : idRes = 25;
30 0 : }
31 0 : else if (higgsType == 1) {
32 0 : nameSave = "f fbar -> h0(H1)";
33 0 : codeSave = 1001;
34 0 : idRes = 25;
35 0 : }
36 0 : else if (higgsType == 2) {
37 0 : nameSave = "f fbar -> H0(H2)";
38 0 : codeSave = 1021;
39 0 : idRes = 35;
40 0 : }
41 0 : else if (higgsType == 3) {
42 0 : nameSave = "f fbar -> A0(A3)";
43 0 : codeSave = 1041;
44 0 : idRes = 36;
45 0 : }
46 :
47 : // Find pointer to H0, H1, H2 or A3 depending on the value of idRes.
48 0 : HResPtr = particleDataPtr->particleDataEntryPtr(idRes);
49 :
50 : // Store H0, H1, H2 or A3 mass and width for propagator.
51 0 : mRes = HResPtr->m0();
52 0 : GammaRes = HResPtr->mWidth();
53 0 : m2Res = mRes*mRes;
54 0 : GamMRat = GammaRes / mRes;
55 :
56 0 : }
57 :
58 :
59 : //--------------------------------------------------------------------------
60 :
61 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
62 :
63 : void Sigma1ffbar2H::sigmaKin() {
64 :
65 : // Set up Breit-Wigner.
66 0 : double width = HResPtr->resWidth(idRes, mH);
67 0 : sigBW = 4. * M_PI/ ( pow2(sH - m2Res) + pow2(mH * width) );
68 :
69 : // Width out only includes open channels.
70 0 : widthOut = width * HResPtr->resOpenFrac(idRes);
71 :
72 0 : }
73 :
74 : //--------------------------------------------------------------------------
75 :
76 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
77 :
78 : double Sigma1ffbar2H::sigmaHat() {
79 :
80 : // Calculate mass-dependent incoming width, including colour factor.
81 0 : int idAbs = abs(id1);
82 0 : double widthIn = HResPtr->resWidthChan( mH, idAbs, -idAbs);
83 0 : if (idAbs < 9) widthIn /= 9.;
84 :
85 : // Done.
86 0 : return widthIn * sigBW * widthOut;
87 :
88 : }
89 :
90 : //--------------------------------------------------------------------------
91 :
92 : // Select identity, colour and anticolour.
93 :
94 : void Sigma1ffbar2H::setIdColAcol() {
95 :
96 : // Flavours trivial.
97 0 : setId( id1, id2, idRes);
98 :
99 : // Colour flow topologies. Swap when antiquarks.
100 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
101 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
102 0 : if (id1 < 0) swapColAcol();
103 :
104 0 : }
105 :
106 : //--------------------------------------------------------------------------
107 :
108 : // Evaluate weight for decay angles.
109 :
110 : double Sigma1ffbar2H::weightDecay( Event& process, int iResBeg,
111 : int iResEnd) {
112 :
113 : // Identity of mother of decaying reseonance(s).
114 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
115 :
116 : // For Higgs decay hand over to standard routine.
117 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
118 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
119 :
120 : // For top decay hand over to standard routine.
121 0 : if (idMother == 6)
122 0 : return weightTopDecay( process, iResBeg, iResEnd);
123 :
124 : // Else done.
125 0 : return 1.;
126 :
127 0 : }
128 :
129 : //==========================================================================
130 :
131 : // Sigma1gg2H class.
132 : // Cross section for g g -> H0, H1, H2 or A3 (H0 SM Higgs, H1, H2, A3 BSM).
133 :
134 : //--------------------------------------------------------------------------
135 :
136 : // Initialize process.
137 :
138 : void Sigma1gg2H::initProc() {
139 :
140 : // Properties specific to Higgs state.
141 0 : if (higgsType == 0) {
142 0 : nameSave = "g g -> H (SM)";
143 0 : codeSave = 902;
144 0 : idRes = 25;
145 0 : }
146 0 : else if (higgsType == 1) {
147 0 : nameSave = "g g -> h0(H1)";
148 0 : codeSave = 1002;
149 0 : idRes = 25;
150 0 : }
151 0 : else if (higgsType == 2) {
152 0 : nameSave = "g g -> H0(H2)";
153 0 : codeSave = 1022;
154 0 : idRes = 35;
155 0 : }
156 0 : else if (higgsType == 3) {
157 0 : nameSave = "g g -> A0(A3)";
158 0 : codeSave = 1042;
159 0 : idRes = 36;
160 0 : }
161 :
162 : // Find pointer to H0, H1, H2 or A3 depending on idRes.
163 0 : HResPtr = particleDataPtr->particleDataEntryPtr(idRes);
164 :
165 : // Store H0, H1, H2 or A3 mass and width for propagator.
166 0 : mRes = HResPtr->m0();
167 0 : GammaRes = HResPtr->mWidth();
168 0 : m2Res = mRes*mRes;
169 0 : GamMRat = GammaRes / mRes;
170 :
171 0 : }
172 :
173 : //--------------------------------------------------------------------------
174 :
175 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
176 :
177 : void Sigma1gg2H::sigmaKin() {
178 :
179 : // Incoming width for gluons, gives colour factor of 1/8 * 1/8.
180 0 : double widthIn = HResPtr->resWidthChan( mH, 21, 21) / 64.;
181 :
182 : // Set up Breit-Wigner.
183 0 : double width = HResPtr->resWidth(idRes, mH);
184 0 : double sigBW = 8. * M_PI/ ( pow2(sH - m2Res) + pow2(mH * width) );
185 :
186 : // Width out only includes open channels.
187 0 : double widthOut = width * HResPtr->resOpenFrac(idRes);
188 :
189 : // Done.
190 0 : sigma = widthIn * sigBW * widthOut;
191 :
192 0 : }
193 :
194 : //--------------------------------------------------------------------------
195 :
196 : // Select identity, colour and anticolour.
197 :
198 : void Sigma1gg2H::setIdColAcol() {
199 :
200 : // Flavours trivial.
201 0 : setId( 21, 21, idRes);
202 :
203 : // Colour flow topology.
204 0 : setColAcol( 1, 2, 2, 1, 0, 0);
205 :
206 0 : }
207 :
208 : //--------------------------------------------------------------------------
209 :
210 : // Evaluate weight for decay angles.
211 :
212 : double Sigma1gg2H::weightDecay( Event& process, int iResBeg,
213 : int iResEnd) {
214 :
215 : // Identity of mother of decaying reseonance(s).
216 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
217 :
218 : // For Higgs decay hand over to standard routine.
219 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
220 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
221 :
222 : // For top decay hand over to standard routine.
223 0 : if (idMother == 6)
224 0 : return weightTopDecay( process, iResBeg, iResEnd);
225 :
226 : // Else done.
227 0 : return 1.;
228 :
229 0 : }
230 :
231 : //==========================================================================
232 :
233 : // Sigma1gmgm2H class.
234 : // Cross section for gamma gamma -> H0, H1, H2 or H3.
235 : // (H0 SM Higgs, H1, H2 and A3 BSM Higgses).
236 :
237 : //--------------------------------------------------------------------------
238 :
239 : // Initialize process.
240 :
241 : void Sigma1gmgm2H::initProc() {
242 :
243 : // Properties specific to Higgs state.
244 0 : if (higgsType == 0) {
245 0 : nameSave = "gamma gamma -> H (SM)";
246 0 : codeSave = 903;
247 0 : idRes = 25;
248 0 : }
249 0 : else if (higgsType == 1) {
250 0 : nameSave = "gamma gamma -> h0(H1)";
251 0 : codeSave = 1003;
252 0 : idRes = 25;
253 0 : }
254 0 : else if (higgsType == 2) {
255 0 : nameSave = "gamma gamma -> H0(H2)";
256 0 : codeSave = 1023;
257 0 : idRes = 35;
258 0 : }
259 0 : else if (higgsType == 3) {
260 0 : nameSave = "gamma gamma -> A0(A3)";
261 0 : codeSave = 1043;
262 0 : idRes = 36;
263 0 : }
264 :
265 : // Find pointer to H0, H1, H2 or A3.
266 0 : HResPtr = particleDataPtr->particleDataEntryPtr(idRes);
267 :
268 : // Store H0, H1, H2 or A3 mass and width for propagator.
269 0 : mRes = HResPtr->m0();
270 0 : GammaRes = HResPtr->mWidth();
271 0 : m2Res = mRes*mRes;
272 0 : GamMRat = GammaRes / mRes;
273 :
274 0 : }
275 :
276 : //--------------------------------------------------------------------------
277 :
278 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
279 :
280 : void Sigma1gmgm2H::sigmaKin() {
281 :
282 : // Incoming width for photons.
283 0 : double widthIn = HResPtr->resWidthChan( mH, 22, 22);
284 :
285 : // Set up Breit-Wigner.
286 0 : double width = HResPtr->resWidth(idRes, mH);
287 0 : double sigBW = 8. * M_PI/ ( pow2(sH - m2Res) + pow2(mH * width) );
288 :
289 : // Width out only includes open channels.
290 0 : double widthOut = width * HResPtr->resOpenFrac(idRes);
291 :
292 : // Done.
293 0 : sigma = widthIn * sigBW * widthOut;
294 :
295 0 : }
296 :
297 : //--------------------------------------------------------------------------
298 :
299 : // Select identity, colour and anticolour.
300 :
301 : void Sigma1gmgm2H::setIdColAcol() {
302 :
303 : // Flavours trivial.
304 0 : setId( 22, 22, idRes);
305 :
306 : // Colour flow trivial.
307 0 : setColAcol( 0, 0, 0, 0, 0, 0);
308 :
309 0 : }
310 :
311 : //--------------------------------------------------------------------------
312 :
313 : // Evaluate weight for decay angles.
314 :
315 : double Sigma1gmgm2H::weightDecay( Event& process, int iResBeg,
316 : int iResEnd) {
317 :
318 : // Identity of mother of decaying reseonance(s).
319 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
320 :
321 : // For Higgs decay hand over to standard routine.
322 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
323 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
324 :
325 : // For top decay hand over to standard routine.
326 0 : if (idMother == 6)
327 0 : return weightTopDecay( process, iResBeg, iResEnd);
328 :
329 : // Else done.
330 0 : return 1.;
331 :
332 0 : }
333 :
334 : //==========================================================================
335 :
336 : // Sigma2ffbar2HZ class.
337 : // Cross section for f fbar -> H0 Z0, H1 Z0, H2 Z0 or A3 Z0.
338 : // (H0 SM Higgs, H1, H2 and A3 BSM Higgses).
339 :
340 : //--------------------------------------------------------------------------
341 :
342 : // Initialize process.
343 :
344 : void Sigma2ffbar2HZ::initProc() {
345 :
346 : // Properties specific to Higgs state.
347 0 : if (higgsType == 0) {
348 0 : nameSave = "f fbar -> H0 Z0 (SM)";
349 0 : codeSave = 904;
350 0 : idRes = 25;
351 0 : coup2Z = 1.;
352 0 : }
353 0 : else if (higgsType == 1) {
354 0 : nameSave = "f fbar -> h0(H1) Z0";
355 0 : codeSave = 1004;
356 0 : idRes = 25;
357 0 : coup2Z = settingsPtr->parm("HiggsH1:coup2Z");
358 0 : }
359 0 : else if (higgsType == 2) {
360 0 : nameSave = "f fbar -> H0(H2) Z0";
361 0 : codeSave = 1024;
362 0 : idRes = 35;
363 0 : coup2Z = settingsPtr->parm("HiggsH2:coup2Z");
364 0 : }
365 0 : else if (higgsType == 3) {
366 0 : nameSave = "f fbar -> A0(A3) ZO";
367 0 : codeSave = 1044;
368 0 : idRes = 36;
369 0 : coup2Z = settingsPtr->parm("HiggsA3:coup2Z");
370 0 : }
371 :
372 : // Store Z0 mass and width for propagator. Common coupling factor.
373 0 : mZ = particleDataPtr->m0(23);
374 0 : widZ = particleDataPtr->mWidth(23);
375 0 : mZS = mZ*mZ;
376 0 : mwZS = pow2(mZ * widZ);
377 0 : thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
378 0 : * couplingsPtr->cos2thetaW());
379 :
380 : // Secondary open width fraction.
381 0 : openFracPair = particleDataPtr->resOpenFrac(idRes, 23);
382 :
383 0 : }
384 :
385 : //--------------------------------------------------------------------------
386 :
387 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
388 :
389 : void Sigma2ffbar2HZ::sigmaKin() {
390 :
391 : // Evaluate differential cross section.
392 0 : sigma0 = (M_PI / sH2) * 8. * pow2(alpEM * thetaWRat * coup2Z)
393 0 : * (tH * uH - s3 * s4 + 2. * sH * s4) / (pow2(sH - mZS) + mwZS);
394 :
395 0 : }
396 :
397 : //--------------------------------------------------------------------------
398 :
399 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
400 :
401 : double Sigma2ffbar2HZ::sigmaHat() {
402 :
403 : // Coupling a_f^2 + v_f^2 to s-channel Z0 and colour factor.
404 0 : int idAbs = abs(id1);
405 0 : double sigma = sigma0 * couplingsPtr->vf2af2(idAbs);
406 0 : if (idAbs < 9) sigma /= 3.;
407 :
408 : // Secondary width for H0 and Z0 or H1 and Z0 or H2 and Z0 or A3 and Z0.
409 0 : sigma *= openFracPair;
410 :
411 : // Answer.
412 0 : return sigma;
413 :
414 : }
415 :
416 : //--------------------------------------------------------------------------
417 :
418 : // Select identity, colour and anticolour.
419 :
420 : void Sigma2ffbar2HZ::setIdColAcol() {
421 :
422 : // Flavours trivial.
423 0 : setId( id1, id2, idRes, 23);
424 :
425 : // Colour flow topologies. Swap when antiquarks.
426 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
427 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
428 0 : if (id1 < 0) swapColAcol();
429 :
430 0 : }
431 :
432 : //--------------------------------------------------------------------------
433 :
434 : // Evaluate weight for decay angles.
435 :
436 : double Sigma2ffbar2HZ::weightDecay( Event& process, int iResBeg,
437 : int iResEnd) {
438 :
439 : // Identity of mother of decaying reseonance(s).
440 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
441 :
442 : // For Higgs decay hand over to standard routine.
443 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
444 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
445 :
446 : // For top decay hand over to standard routine.
447 0 : if (idMother == 6)
448 0 : return weightTopDecay( process, iResBeg, iResEnd);
449 :
450 : // If not decay of Z0 created along with Higgs then done.
451 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
452 :
453 : // Order so that fbar(1) f(2) -> H() f'(3) fbar'(4).
454 0 : int i1 = (process[3].id() < 0) ? 3 : 4;
455 0 : int i2 = 7 - i1;
456 0 : int i3 = process[6].daughter1();
457 0 : int i4 = process[6].daughter2();
458 0 : if (process[i3].id() < 0) swap( i3, i4);
459 :
460 : // Find left- and righthanded couplings of fermion pairs.
461 0 : int idAbs = process[i1].idAbs();
462 0 : double liS = pow2( couplingsPtr->lf(idAbs) );
463 0 : double riS = pow2( couplingsPtr->rf(idAbs) );
464 0 : idAbs = process[i3].idAbs();
465 0 : double lfS = pow2( couplingsPtr->lf(idAbs) );
466 0 : double rfS = pow2( couplingsPtr->rf(idAbs) );
467 :
468 : // Evaluate relevant four-products.
469 0 : double pp13 = process[i1].p() * process[i3].p();
470 0 : double pp14 = process[i1].p() * process[i4].p();
471 0 : double pp23 = process[i2].p() * process[i3].p();
472 0 : double pp24 = process[i2].p() * process[i4].p();
473 :
474 : // Weight and maximum.
475 0 : double wt = (liS * lfS + riS * rfS) * pp13 * pp24
476 0 : + (liS * rfS + riS * lfS) * pp14 * pp23;
477 0 : double wtMax = (liS + riS) * (lfS + rfS) * (pp13 + pp14) * (pp23 + pp24);
478 :
479 : // Done.
480 0 : return wt / wtMax;
481 :
482 0 : }
483 :
484 : //==========================================================================
485 :
486 : // Sigma2ffbar2HW class.
487 : // Cross section for f fbar -> H0 W+-, H1 W+-, H2 W+- or A3 W+-.
488 : // (H0 SM Higgs, H1, H2 and A3 BSM Higgses).
489 :
490 : //--------------------------------------------------------------------------
491 :
492 : // Initialize process.
493 :
494 : void Sigma2ffbar2HW::initProc() {
495 :
496 : // Properties specific to Higgs state.
497 0 : if (higgsType == 0) {
498 0 : nameSave = "f fbar -> H0 W+- (SM)";
499 0 : codeSave = 905;
500 0 : idRes = 25;
501 0 : coup2W = 1.;
502 0 : }
503 0 : else if (higgsType == 1) {
504 0 : nameSave = "f fbar -> h0(H1) W+-";
505 0 : codeSave = 1005;
506 0 : idRes = 25;
507 0 : coup2W = settingsPtr->parm("HiggsH1:coup2W");
508 0 : }
509 0 : else if (higgsType == 2) {
510 0 : nameSave = "f fbar -> H0(H2) W+-";
511 0 : codeSave = 1025;
512 0 : idRes = 35;
513 0 : coup2W = settingsPtr->parm("HiggsH2:coup2W");
514 0 : }
515 0 : else if (higgsType == 3) {
516 0 : nameSave = "f fbar -> A0(A3) W+-";
517 0 : codeSave = 1045;
518 0 : idRes = 36;
519 0 : coup2W = settingsPtr->parm("HiggsA3:coup2W");
520 0 : }
521 :
522 : // Store W+- mass and width for propagator. Common coupling factor.
523 0 : mW = particleDataPtr->m0(24);
524 0 : widW = particleDataPtr->mWidth(24);
525 0 : mWS = mW*mW;
526 0 : mwWS = pow2(mW * widW);
527 0 : thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
528 :
529 : // Secondary open width fractions.
530 0 : openFracPairPos = particleDataPtr->resOpenFrac(idRes, 24);
531 0 : openFracPairNeg = particleDataPtr->resOpenFrac(idRes, -24);
532 :
533 0 : }
534 :
535 : //--------------------------------------------------------------------------
536 :
537 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
538 :
539 : void Sigma2ffbar2HW::sigmaKin() {
540 :
541 : // Evaluate differential cross section.
542 0 : sigma0 = (M_PI / sH2) * 2. * pow2(alpEM * thetaWRat * coup2W)
543 0 : * (tH * uH - s3 * s4 + 2. * sH * s4) / (pow2(sH - mWS) + mwWS);
544 :
545 0 : }
546 :
547 : //--------------------------------------------------------------------------
548 :
549 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
550 :
551 : double Sigma2ffbar2HW::sigmaHat() {
552 :
553 : // CKM and colour factors.
554 0 : double sigma = sigma0;
555 0 : if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
556 :
557 : // Secondary width for H0 and W+-.
558 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
559 0 : sigma *= (idUp > 0) ? openFracPairPos : openFracPairNeg;
560 :
561 : // Answer.
562 0 : return sigma;
563 :
564 : }
565 :
566 : //--------------------------------------------------------------------------
567 :
568 : // Select identity, colour and anticolour.
569 :
570 : void Sigma2ffbar2HW::setIdColAcol() {
571 :
572 : // Sign of outgoing W.
573 0 : int sign = 1 - 2 * (abs(id1)%2);
574 0 : if (id1 < 0) sign = -sign;
575 0 : setId( id1, id2, idRes, 24 * sign);
576 :
577 : // Colour flow topologies. Swap when antiquarks.
578 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
579 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
580 0 : if (id1 < 0) swapColAcol();
581 :
582 0 : }
583 :
584 : //--------------------------------------------------------------------------
585 :
586 : // Evaluate weight for decay angles.
587 :
588 : double Sigma2ffbar2HW::weightDecay( Event& process, int iResBeg,
589 : int iResEnd) {
590 :
591 : // Identity of mother of decaying reseonance(s).
592 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
593 :
594 : // For Higgs decay hand over to standard routine.
595 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
596 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
597 :
598 : // For top decay hand over to standard routine.
599 0 : if (idMother == 6)
600 0 : return weightTopDecay( process, iResBeg, iResEnd);
601 :
602 : // If not decay of W+- created along with Higgs then done.
603 0 : if (iResBeg != 5 || iResEnd != 6) return 1.;
604 :
605 : // Order so that fbar(1) f(2) -> H() f'(3) fbar'(4).
606 0 : int i1 = (process[3].id() < 0) ? 3 : 4;
607 0 : int i2 = 7 - i1;
608 0 : int i3 = process[6].daughter1();
609 0 : int i4 = process[6].daughter2();
610 0 : if (process[i3].id() < 0) swap( i3, i4);
611 :
612 : // Evaluate relevant four-products.
613 0 : double pp13 = process[i1].p() * process[i3].p();
614 0 : double pp14 = process[i1].p() * process[i4].p();
615 0 : double pp23 = process[i2].p() * process[i3].p();
616 0 : double pp24 = process[i2].p() * process[i4].p();
617 :
618 : // Weight and maximum.
619 0 : double wt = pp13 * pp24;
620 0 : double wtMax = (pp13 + pp14) * (pp23 + pp24);
621 :
622 : // Done.
623 0 : return wt / wtMax;
624 :
625 0 : }
626 :
627 : //==========================================================================
628 :
629 : // Sigma3ff2HfftZZ class.
630 : // Cross section for f f' -> H f f' (Z0 Z0 fusion of SM or BSM Higgs).
631 : // (H can be H0 SM or H1, H2, A3 from BSM).
632 :
633 : //--------------------------------------------------------------------------
634 :
635 : // Initialize process.
636 :
637 : void Sigma3ff2HfftZZ::initProc() {
638 :
639 : // Properties specific to Higgs state.
640 0 : if (higgsType == 0) {
641 0 : nameSave = "f f' -> H0 f f'(Z0 Z0 fusion) (SM)";
642 0 : codeSave = 906;
643 0 : idRes = 25;
644 0 : coup2Z = 1.;
645 0 : }
646 0 : else if (higgsType == 1) {
647 0 : nameSave = "f f' -> h0(H1) f f' (Z0 Z0 fusion)";
648 0 : codeSave = 1006;
649 0 : idRes = 25;
650 0 : coup2Z = settingsPtr->parm("HiggsH1:coup2Z");
651 0 : }
652 0 : else if (higgsType == 2) {
653 0 : nameSave = "f f' -> H0(H2) f f' (Z0 Z0 fusion)";
654 0 : codeSave = 1026;
655 0 : idRes = 35;
656 0 : coup2Z = settingsPtr->parm("HiggsH2:coup2Z");
657 0 : }
658 0 : else if (higgsType == 3) {
659 0 : nameSave = "f f' -> A0(A3) f f' (Z0 Z0 fusion)";
660 0 : codeSave = 1046;
661 0 : idRes = 36;
662 0 : coup2Z = settingsPtr->parm("HiggsA3:coup2Z");
663 0 : }
664 :
665 : // Common fixed mass and coupling factor.
666 0 : mZS = pow2( particleDataPtr->m0(23) );
667 0 : prefac = 0.25 * mZS * pow3( 4. * M_PI / (couplingsPtr->sin2thetaW()
668 0 : * couplingsPtr->cos2thetaW()) );
669 :
670 : // Secondary open width fraction.
671 0 : openFrac = particleDataPtr->resOpenFrac(idRes);
672 :
673 0 : }
674 :
675 : //--------------------------------------------------------------------------
676 :
677 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
678 :
679 : void Sigma3ff2HfftZZ::sigmaKin() {
680 :
681 : // Required four-vector products.
682 0 : double pp12 = 0.5 * sH;
683 0 : double pp14 = 0.5 * mH * p4cm.pNeg();
684 0 : double pp15 = 0.5 * mH * p5cm.pNeg();
685 0 : double pp24 = 0.5 * mH * p4cm.pPos();
686 0 : double pp25 = 0.5 * mH * p5cm.pPos();
687 0 : double pp45 = p4cm * p5cm;
688 :
689 : // Propagator factors and two possible numerators.
690 0 : double prop = pow2( (2. * pp14 + mZS) * (2. * pp25 + mZS) );
691 0 : sigma1 = prefac * pp12 * pp45 / prop;
692 0 : sigma2 = prefac * pp15 * pp24 / prop;
693 :
694 0 : }
695 :
696 : //--------------------------------------------------------------------------
697 :
698 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
699 :
700 : double Sigma3ff2HfftZZ::sigmaHat() {
701 :
702 : // Flavour-dependent coupling factors for two incoming flavours.
703 0 : int id1Abs = abs(id1);
704 0 : int id2Abs = abs(id2);
705 0 : double lf1S = pow2( couplingsPtr->lf(id1Abs) );
706 0 : double rf1S = pow2( couplingsPtr->rf(id1Abs) );
707 0 : double lf2S = pow2( couplingsPtr->lf(id2Abs) );
708 0 : double rf2S = pow2( couplingsPtr->rf(id2Abs) );
709 0 : double c1 = lf1S * lf2S + rf1S * rf2S;
710 0 : double c2 = lf1S * rf2S + rf1S * lf2S;
711 :
712 : // Combine couplings and kinematics factors.
713 0 : double sigma = pow3(alpEM) * (c1 * sigma1 + c2 * sigma2) * pow2(coup2Z);
714 :
715 : // Secondary width for H0, H1, H2 or A3.
716 0 : sigma *= openFrac;
717 :
718 : // Answer.
719 0 : return sigma;
720 :
721 : }
722 :
723 : //--------------------------------------------------------------------------
724 :
725 : // Select identity, colour and anticolour.
726 :
727 : void Sigma3ff2HfftZZ::setIdColAcol() {
728 :
729 : // Trivial flavours: out = in.
730 0 : setId( id1, id2, idRes, id1, id2);
731 :
732 : // Colour flow topologies. Swap when antiquarks.
733 0 : if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
734 0 : setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
735 0 : else if (abs(id1) < 9 && abs(id2) < 9)
736 0 : setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
737 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
738 0 : else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0);
739 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
740 0 : if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
741 0 : swapColAcol();
742 :
743 0 : }
744 :
745 : //--------------------------------------------------------------------------
746 :
747 : // Evaluate weight for decay angles.
748 :
749 : double Sigma3ff2HfftZZ::weightDecay( Event& process, int iResBeg,
750 : int iResEnd) {
751 :
752 : // Identity of mother of decaying reseonance(s).
753 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
754 :
755 : // For Higgs decay hand over to standard routine.
756 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
757 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
758 :
759 : // For top decay hand over to standard routine.
760 0 : if (idMother == 6)
761 0 : return weightTopDecay( process, iResBeg, iResEnd);
762 :
763 : // Else done.
764 0 : return 1.;
765 :
766 0 : }
767 :
768 : //==========================================================================
769 :
770 : // Sigma3ff2HfftWW class.
771 : // Cross section for f_1 f_2 -> H0 f_3 f_4 (W+ W- fusion of SM or BSM Higgs).
772 :
773 : //--------------------------------------------------------------------------
774 :
775 : // Initialize process.
776 :
777 : void Sigma3ff2HfftWW::initProc() {
778 :
779 : // Properties specific to Higgs state.
780 0 : if (higgsType == 0) {
781 0 : nameSave = "f_1 f_2 -> H0 f_3 f_4 (W+ W- fusion) (SM)";
782 0 : codeSave = 907;
783 0 : idRes = 25;
784 0 : coup2W = 1.;
785 0 : }
786 0 : else if (higgsType == 1) {
787 0 : nameSave = "f_1 f_2 -> h0(H1) f_3 f_4 (W+ W- fusion)";
788 0 : codeSave = 1007;
789 0 : idRes = 25;
790 0 : coup2W = settingsPtr->parm("HiggsH1:coup2W");
791 0 : }
792 0 : else if (higgsType == 2) {
793 0 : nameSave = "f_1 f_2 -> H0(H2) f_3 f_4 (W+ W- fusion)";
794 0 : codeSave = 1027;
795 0 : idRes = 35;
796 0 : coup2W = settingsPtr->parm("HiggsH2:coup2W");
797 0 : }
798 0 : else if (higgsType == 3) {
799 0 : nameSave = "f_1 f_2 -> A0(A3) f_3 f_4 (W+ W- fusion)";
800 0 : codeSave = 1047;
801 0 : idRes = 36;
802 0 : coup2W = settingsPtr->parm("HiggsA3:coup2W");
803 0 : }
804 :
805 : // Common fixed mass and coupling factor.
806 0 : mWS = pow2( particleDataPtr->m0(24) );
807 0 : prefac = mWS * pow3( 4. * M_PI / couplingsPtr->sin2thetaW() );
808 :
809 : // Secondary open width fraction.
810 0 : openFrac = particleDataPtr->resOpenFrac(idRes);
811 :
812 0 : }
813 :
814 : //--------------------------------------------------------------------------
815 :
816 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
817 :
818 : void Sigma3ff2HfftWW::sigmaKin() {
819 :
820 : // Required four-vector products.
821 0 : double pp12 = 0.5 * sH;
822 0 : double pp14 = 0.5 * mH * p4cm.pNeg();
823 0 : double pp25 = 0.5 * mH * p5cm.pPos();
824 0 : double pp45 = p4cm * p5cm;
825 :
826 : // Cross section: kinematics part. Combine with couplings.
827 0 : double prop = pow2( (2. * pp14 + mWS) * (2. * pp25 + mWS) );
828 0 : sigma0 = prefac * pp12 * pp45 * pow2(coup2W) / prop;
829 :
830 0 : }
831 :
832 : //--------------------------------------------------------------------------
833 :
834 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
835 :
836 : double Sigma3ff2HfftWW::sigmaHat() {
837 :
838 : // Some flavour combinations not possible.
839 0 : int id1Abs = abs(id1);
840 0 : int id2Abs = abs(id2);
841 0 : if ( (id1Abs%2 == id2Abs%2 && id1 * id2 > 0)
842 0 : || (id1Abs%2 != id2Abs%2 && id1 * id2 < 0) ) return 0.;
843 :
844 : // Basic cross section. CKM factors for final states.
845 0 : double sigma = sigma0 * pow3(alpEM) * couplingsPtr->V2CKMsum(id1Abs)
846 0 : * couplingsPtr->V2CKMsum(id2Abs);
847 :
848 : // Secondary width for H0, H1, H2 or A3.
849 0 : sigma *= openFrac;
850 :
851 : // Spin-state extra factor 2 per incoming neutrino.
852 0 : if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
853 0 : if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
854 :
855 : // Answer.
856 : return sigma;
857 :
858 0 : }
859 :
860 : //--------------------------------------------------------------------------
861 :
862 : // Select identity, colour and anticolour.
863 :
864 : void Sigma3ff2HfftWW::setIdColAcol() {
865 :
866 : // Pick out-flavours by relative CKM weights.
867 0 : id4 = couplingsPtr->V2CKMpick(id1);
868 0 : id5 = couplingsPtr->V2CKMpick(id2);
869 0 : setId( id1, id2, idRes, id4, id5);
870 :
871 : // Colour flow topologies. Swap when antiquarks.
872 0 : if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
873 0 : setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
874 0 : else if (abs(id1) < 9 && abs(id2) < 9)
875 0 : setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
876 0 : else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
877 0 : else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0);
878 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
879 0 : if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
880 0 : swapColAcol();
881 :
882 0 : }
883 :
884 : //--------------------------------------------------------------------------
885 :
886 : // Evaluate weight for decay angles.
887 :
888 : double Sigma3ff2HfftWW::weightDecay( Event& process, int iResBeg,
889 : int iResEnd) {
890 :
891 : // Identity of mother of decaying reseonance(s).
892 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
893 :
894 : // For Higgs decay hand over to standard routine.
895 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
896 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
897 :
898 : // For top decay hand over to standard routine.
899 0 : if (idMother == 6)
900 0 : return weightTopDecay( process, iResBeg, iResEnd);
901 :
902 : // Else done.
903 0 : return 1.;
904 :
905 0 : }
906 :
907 : //==========================================================================
908 :
909 : // Sigma3gg2HQQbar class.
910 : // Cross section for g g -> H0 Q Qbar (Q Qbar fusion of SM or BSM Higgs).
911 :
912 : //--------------------------------------------------------------------------
913 :
914 : // Initialize process.
915 :
916 : void Sigma3gg2HQQbar::initProc() {
917 :
918 : // Properties specific to Higgs state for the "g g -> H ttbar" process.
919 : // (H can be H0 SM or H1, H2, A3 from BSM).
920 0 : if (higgsType == 0 && idNew == 6) {
921 0 : nameSave = "g g -> H t tbar (SM)";
922 0 : codeSave = 908;
923 0 : idRes = 25;
924 0 : coup2Q = 1.;
925 0 : }
926 0 : else if (higgsType == 1 && idNew == 6) {
927 0 : nameSave = "g g -> h0(H1) t tbar";
928 0 : codeSave = 1008;
929 0 : idRes = 25;
930 0 : coup2Q = settingsPtr->parm("HiggsH1:coup2u");
931 0 : }
932 0 : else if (higgsType == 2 && idNew == 6) {
933 0 : nameSave = "g g -> H0(H2) t tbar";
934 0 : codeSave = 1028;
935 0 : idRes = 35;
936 0 : coup2Q = settingsPtr->parm("HiggsH2:coup2u");
937 0 : }
938 0 : else if (higgsType == 3 && idNew == 6) {
939 0 : nameSave = "g g -> A0(A3) t tbar";
940 0 : codeSave = 1048;
941 0 : idRes = 36;
942 0 : coup2Q = settingsPtr->parm("HiggsA3:coup2u");
943 0 : }
944 :
945 : // Properties specific to Higgs state for the "g g -> H b bbar" process.
946 : // (H can be H0 SM or H1, H2, A3 from BSM).
947 0 : if (higgsType == 0 && idNew == 5) {
948 0 : nameSave = "g g -> H b bbar (SM)";
949 0 : codeSave = 912;
950 0 : idRes = 25;
951 0 : coup2Q = 1.;
952 0 : }
953 0 : else if (higgsType == 1 && idNew == 5) {
954 0 : nameSave = "g g -> h0(H1) b bbar";
955 0 : codeSave = 1012;
956 0 : idRes = 25;
957 0 : coup2Q = settingsPtr->parm("HiggsH1:coup2d");
958 0 : }
959 0 : else if (higgsType == 2 && idNew == 5) {
960 0 : nameSave = "g g -> H0(H2) b bbar";
961 0 : codeSave = 1032;
962 0 : idRes = 35;
963 0 : coup2Q = settingsPtr->parm("HiggsH2:coup2d");
964 0 : }
965 0 : else if (higgsType == 3 && idNew == 5) {
966 0 : nameSave = "g g -> A0(A3) b bbar";
967 0 : codeSave = 1052;
968 0 : idRes = 36;
969 0 : coup2Q = settingsPtr->parm("HiggsA3:coup2d");
970 0 : }
971 :
972 : // Common mass and coupling factors.
973 0 : double mWS = pow2(particleDataPtr->m0(24));
974 0 : prefac = (4. * M_PI / couplingsPtr->sin2thetaW()) * pow2(4. * M_PI)
975 0 : * 0.25 / mWS;
976 :
977 : // Secondary open width fraction.
978 0 : openFracTriplet = particleDataPtr->resOpenFrac(idRes, idNew, -idNew);
979 :
980 0 : }
981 :
982 : //--------------------------------------------------------------------------
983 :
984 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
985 :
986 : void Sigma3gg2HQQbar::sigmaKin() {
987 :
988 : // Running mass of heavy quark.
989 0 : double mQ2run = pow2( particleDataPtr->mRun(idNew, mH) );
990 :
991 : // Linear combination of p_Q and p_Qbar to ensure common mass.
992 0 : double mQ2 = m4 * m5;
993 : double epsi = 0.;
994 0 : if (m4 != m5) {
995 0 : double s45 = (p4cm + p5cm).m2Calc();
996 0 : mQ2 = 0.5 * (s4 + s5) - 0.25 * pow2(s4 - s5) / s45;
997 0 : epsi = 0.5 * (s5 - s4) / s45;
998 0 : }
999 :
1000 : // Set up kinematics: g(4) g(5) -> H(3) Q(1) Qbar(2) in outgoing sense.
1001 0 : Vec4 pTemp[6];
1002 0 : pTemp[4] = Vec4( 0., 0., -0.5* mH, -0.5* mH);
1003 0 : pTemp[5] = Vec4( 0., 0., 0.5* mH, -0.5* mH);
1004 0 : pTemp[1] = p4cm + epsi * (p4cm + p5cm);
1005 0 : pTemp[2] = p5cm - epsi * (p4cm + p5cm);
1006 0 : pTemp[3] = p3cm;
1007 :
1008 : // Four-product combinations.
1009 0 : double z1 = pTemp[1] * pTemp[2];
1010 0 : double z2 = pTemp[1] * pTemp[3];
1011 0 : double z3 = pTemp[1] * pTemp[4];
1012 0 : double z4 = pTemp[1] * pTemp[5];
1013 0 : double z5 = pTemp[2] * pTemp[3];
1014 0 : double z6 = pTemp[2] * pTemp[4];
1015 0 : double z7 = pTemp[2] * pTemp[5];
1016 0 : double z8 = pTemp[3] * pTemp[4];
1017 0 : double z9 = pTemp[3] * pTemp[5];
1018 0 : double z10 = pTemp[4] * pTemp[5];
1019 :
1020 : // Powers required as shorthand in matriz elements.
1021 0 : double mQ4 = mQ2 * mQ2;
1022 0 : double mQ6 = mQ2 * mQ4;
1023 0 : double z1S = z1 * z1;
1024 0 : double z2S = z2 * z2;
1025 0 : double z3S = z3 * z3;
1026 0 : double z4S = z4 * z4;
1027 0 : double z5S = z5 * z5;
1028 0 : double z6S = z6 * z6;
1029 0 : double z7S = z7 * z7;
1030 0 : double z8S = z8 * z8;
1031 0 : double z9S = z9 * z9;
1032 0 : double z10S = z10 * z10;
1033 :
1034 : // Evaluate matriz elements for g + g -> Q + Qbar + H.
1035 : // (H can be H0 SM or H1, H2, A3 from BSM).
1036 0 : double fm[9][9];
1037 0 : fm[1][1] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+2*z2+z4+z9+2*
1038 0 : z7+z5)+8*mQ2*s3*(-z1-z4+2*z7)+16*mQ2*(z2*z9+4*z2*
1039 0 : z7+z2*z5-2*z4*z7-2*z9*z7)+8*s3*z4*z7-16*z2*z9*z7;
1040 0 : fm[1][2] = 16*mQ6+8*mQ4*(-2*z1+z2-2*z3-2*z4-4*z10+z9-z8+2
1041 0 : *z7-4*z6+z5)+8*mQ2*(-2*z1*z2-2*z2*z4-2*z2*z10+z2*z7-2*
1042 0 : z2*z6-2*z3*z7+2*z4*z7+4*z10*z7-z9*z7-z8*z7)+16*z2*z7*(z4+
1043 : z10);
1044 0 : fm[1][3] = 16*mQ6-4*mQ4*s3+8*mQ4*(-2*z1+2*z2-2*z3-4*
1045 0 : z4-8*z10+z9+z8-2*z7-4*z6+2*z5)-(4*mQ2*s3)*(z1+z4+z10
1046 0 : +z6)+8*mQ2*(-2*z1*z2-2*z1*z10+z1*z9+z1*z8-2*z1*z5+z2S
1047 0 : -4*z2*z4-5*z2*z10+z2*z8-z2*z7-3*z2*z6+z2*z5+z3*z9+2*z3*z7
1048 0 : -z3*z5+z4*z8+2*z4*z6-3*z4*z5-5*z10*z5+z9*z8+z9*z6+z9*z5+
1049 0 : z8*z7-4*z6*z5+z5S)-(16*z2*z5)*(z1+z4+z10+z6);
1050 0 : fm[1][4] = 16*mQ6+4*mQ4*s3+16*mQ4*(-z1+z2-z3-z4+z10-
1051 0 : z9-z8+2*z7+2*z6-z5)+4*mQ2*s3*(z1+z3+z4+z10+2*z7+2*z6
1052 0 : )+8*mQ2*(4*z1*z10+4*z1*z7+4*z1*z6+2*z2*z10-z2*z9-z2*z8+
1053 0 : 4*z2*z7+4*z2*z6-z2*z5+4*z10*z5+4*z7*z5+4*z6*z5)-(8*s3*
1054 0 : z1)*(z10+z7+z6)+16*z2*z5*(z10+z7+z6);
1055 0 : fm[1][5] = 8*mQ4*(-2*z1-2*z4+z10-z9)+4*mQ2*(4*z1S-2*z1*
1056 0 : z2+8*z1*z3+6*z1*z10-2*z1*z9+4*z1*z8+4*z1*z7+4*z1*z6+2*z1*
1057 0 : z5+z2*z10+4*z3*z4-z3*z9+2*z3*z7+3*z4*z8-2*z4*z6+2*z4*z5-4
1058 0 : *z10*z7+3*z10*z5-3*z9*z6+3*z8*z7-4*z7S+4*z7*z5)+8*(z1S
1059 0 : *z9-z1S*z8-z1*z2*z7+z1*z2*z6+z1*z3*z9+z1*z3*z5-z1*z4*
1060 0 : z8-z1*z4*z5+z1*z10*z9+z1*z9*z7+z1*z9*z6-z1*z8*z7-z2*z3*z7
1061 0 : +z2*z4*z6-z2*z10*z7-z2*z7S+z3*z7*z5-z4*z10*z5-z4*z7*z5-
1062 0 : z4*z6*z5);
1063 0 : fm[1][6] = 16*mQ4*(-4*z1-z4+z9-z7)+4*mQ2*s3*(-2*z1-z4-
1064 0 : z7)+16*mQ2*(-2*z1S-3*z1*z2-2*z1*z4-3*z1*z9-2*z1*z7-3*
1065 0 : z1*z5-2*z2*z4-2*z7*z5)-8*s3*z4*z7+8*(-z1*z2*z9-2*z1*z2
1066 0 : *z5-z1*z9S-z1*z9*z5+z2S*z7-z2*z4*z5+z2*z9*z7-z2*z7*z5
1067 0 : +z4*z9*z5+z4*z5S);
1068 0 : fm[1][7] = 8*mQ4*(2*z3+z4+3*z10+z9+2*z8+3*z7+6*z6)+2*mQ2*
1069 0 : s3*(-2*z3-z4+3*z10+3*z7+6*z6)+4*mQ2*(4*z1*z10+4*z1*
1070 0 : z7+8*z1*z6+6*z2*z10+z2*z9+2*z2*z8+6*z2*z7+12*z2*z6-8*z3*
1071 0 : z7+4*z4*z7+4*z4*z6+4*z10*z5+4*z9*z7+4*z9*z6-8*z8*z7+4*z7*
1072 0 : z5+8*z6*z5)+4*s3*(-z1*z10-z1*z7-2*z1*z6+2*z3*z7-z4*z7-
1073 0 : z4*z6)+8*z2*(z10*z5+z9*z7+z9*z6-2*z8*z7+z7*z5+2*z6*z5);
1074 0 : fm[1][8] = 8*mQ4*(2*z3+z4+3*z10+2*z9+z8+3*z7+6*z6)+2*mQ2*
1075 0 : s3*(-2*z3-z4+2*z10+z7+2*z6)+4*mQ2*(4*z1*z10-2*z1*z9+
1076 0 : 2*z1*z8+4*z1*z7+8*z1*z6+5*z2*z10+2*z2*z9+z2*z8+4*z2*z7+8*
1077 0 : z2*z6-z3*z9-8*z3*z7+2*z3*z5+2*z4*z9-z4*z8+4*z4*z7+4*z4*z6
1078 0 : +4*z4*z5+5*z10*z5+z9S-z9*z8+2*z9*z7+5*z9*z6+z9*z5-7*z8*
1079 0 : z7+2*z8*z5+2*z7*z5+10*z6*z5)+2*s3*(-z1*z10+z3*z7-2*z4*
1080 0 : z7+z4*z6)+4*(-z1*z9S+z1*z9*z8-2*z1*z9*z5-z1*z8*z5+2*z2*
1081 0 : z10*z5+z2*z9*z7+z2*z9*z6-2*z2*z8*z7+3*z2*z6*z5+z3*z9*z5+
1082 0 : z3*z5S+z4*z9*z5-2*z4*z8*z5+2*z4*z5S);
1083 0 : fm[2][2] = 16*mQ6+16*mQ4*(-z1+z3-z4-z10+z7-z6)+16*mQ2*(
1084 0 : z3*z10+z3*z7+z3*z6+z4*z7+z10*z7)-16*z3*z10*z7;
1085 0 : fm[2][3] = 16*mQ6+8*mQ4*(-2*z1+z2+2*z3-4*z4-4*z10-z9+z8-2
1086 0 : *z7-2*z6+z5)+8*mQ2*(-2*z1*z5+4*z3*z10-z3*z9-z3*z8-2*z3*
1087 0 : z7+2*z3*z6+z3*z5-2*z4*z5-2*z10*z5-2*z6*z5)+16*z3*z5*(z10+
1088 : z6);
1089 0 : fm[2][4] = 8*mQ4*(-2*z1-2*z3+z10-z8)+4*mQ2*(4*z1S-2*z1*
1090 0 : z2+8*z1*z4+6*z1*z10+4*z1*z9-2*z1*z8+4*z1*z7+4*z1*z6+2*z1*
1091 0 : z5+z2*z10+4*z3*z4+3*z3*z9-2*z3*z7+2*z3*z5-z4*z8+2*z4*z6-4
1092 0 : *z10*z6+3*z10*z5+3*z9*z6-3*z8*z7-4*z6S+4*z6*z5)+8*(-z1S
1093 0 : *z9+z1S*z8+z1*z2*z7-z1*z2*z6-z1*z3*z9-z1*z3*z5+z1*z4
1094 0 : *z8+z1*z4*z5+z1*z10*z8-z1*z9*z6+z1*z8*z7+z1*z8*z6+z2*z3*
1095 0 : z7-z2*z4*z6-z2*z10*z6-z2*z6S-z3*z10*z5-z3*z7*z5-z3*z6*
1096 0 : z5+z4*z6*z5);
1097 0 : fm[2][5] = 16*mQ4*z10+8*mQ2*(2*z1S+2*z1*z3+2*z1*z4+2*z1
1098 0 : *z10+2*z1*z7+2*z1*z6+z3*z7+z4*z6)+8*(-2*pow3(z1)-2*z1S*z3-
1099 0 : 2*z1S*z4-2*z1S*z10-2*z1S*z7-2*z1S*z6-2*z1*z3*z4-
1100 0 : z1*z3*z10-2*z1*z3*z6-z1*z4*z10-2*z1*z4*z7-z1*z10S-z1*
1101 0 : z10*z7-z1*z10*z6-2*z1*z7*z6+z3S*z7-z3*z4*z7-z3*z4*z6+z3
1102 0 : *z10*z7+z3*z7S-z3*z7*z6+z4S*z6+z4*z10*z6-z4*z7*z6+z4*
1103 : z6S);
1104 0 : fm[2][6] = 8*mQ4*(-2*z1+z10-z9-2*z7)+4*mQ2*(4*z1S+2*z1*
1105 0 : z2+4*z1*z3+4*z1*z4+6*z1*z10-2*z1*z9+4*z1*z8+8*z1*z6-2*z1*
1106 0 : z5+4*z2*z4+3*z2*z10+2*z2*z7-3*z3*z9-2*z3*z7-4*z4S-4*z4*
1107 0 : z10+3*z4*z8+2*z4*z6+z10*z5-z9*z6+3*z8*z7+4*z7*z6)+8*(z1S
1108 0 : *z9-z1S*z8-z1*z2*z7+z1*z2*z6+z1*z3*z9+z1*z3*z5+z1*z4*
1109 0 : z9-z1*z4*z8-z1*z4*z5+z1*z10*z9+z1*z9*z6-z1*z8*z7-z2*z3*z7
1110 0 : -z2*z4*z7+z2*z4*z6-z2*z10*z7+z3*z7*z5-z4S*z5-z4*z10*z5-
1111 : z4*z6*z5);
1112 0 : fm[2][7] = 8*mQ4*(z3+2*z4+3*z10+z7+2*z6)+4*mQ2*(-4*z1*z3-
1113 0 : 2*z1*z4-2*z1*z10+z1*z9-z1*z8-4*z1*z7-2*z1*z6+z2*z3+2*z2*
1114 0 : z4+3*z2*z10+z2*z7+2*z2*z6-6*z3*z4-6*z3*z10-2*z3*z9-2*z3*
1115 0 : z7-4*z3*z6-z3*z5-6*z4S-6*z4*z10-3*z4*z9-z4*z8-4*z4*z7-2
1116 0 : *z4*z6-2*z4*z5-3*z10*z9-3*z10*z8-6*z10*z7-6*z10*z6+z10*z5
1117 0 : +z9*z7-2*z8*z7-2*z8*z6-6*z7*z6+z7*z5-6*z6S+2*z6*z5)+4*(
1118 0 : -z1S*z9+z1S*z8-2*z1*z2*z10-3*z1*z2*z7-3*z1*z2*z6+z1*
1119 0 : z3*z9-z1*z3*z5+z1*z4*z9+z1*z4*z8+z1*z4*z5+z1*z10*z9+z1*
1120 0 : z10*z8-z1*z9*z6+z1*z8*z6+z2*z3*z7-3*z2*z4*z7-z2*z4*z6-3*
1121 0 : z2*z10*z7-3*z2*z10*z6-3*z2*z7*z6-3*z2*z6S-2*z3*z4*z5-z3
1122 0 : *z10*z5-z3*z6*z5-z4S*z5-z4*z10*z5+z4*z6*z5);
1123 0 : fm[2][8] = 8*mQ4*(z3+2*z4+3*z10+z7+2*z6)+4*mQ2*(-4*z1*z3-
1124 0 : 2*z1*z4-2*z1*z10-z1*z9+z1*z8-4*z1*z7-2*z1*z6+z2*z3+2*z2*
1125 0 : z4+z2*z10-z2*z7-2*z2*z6-6*z3*z4-6*z3*z10-2*z3*z9+z3*z8-2*
1126 0 : z3*z7-4*z3*z6+z3*z5-6*z4S-6*z4*z10-2*z4*z9-4*z4*z7-2*z4
1127 0 : *z6+2*z4*z5-3*z10*z9-3*z10*z8-6*z10*z7-6*z10*z6+3*z10*z5-
1128 0 : z9*z6-2*z8*z7-3*z8*z6-6*z7*z6+z7*z5-6*z6S+2*z6*z5)+4*(
1129 0 : z1S*z9-z1S*z8-z1*z2*z7+z1*z2*z6-3*z1*z3*z5+z1*z4*z9-
1130 0 : z1*z4*z8-3*z1*z4*z5+z1*z10*z9+z1*z10*z8-2*z1*z10*z5+z1*z9
1131 0 : *z6+z1*z8*z7+z1*z8*z6-z2*z4*z7+z2*z4*z6-z2*z10*z7-z2*z10*
1132 0 : z6-2*z2*z7*z6-z2*z6S-3*z3*z4*z5-3*z3*z10*z5+z3*z7*z5-3*
1133 0 : z3*z6*z5-3*z4S*z5-3*z4*z10*z5-z4*z6*z5);
1134 0 : fm[3][3] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+z2+2*z3+z8+z6
1135 0 : +2*z5)+8*mQ2*s3*(-z1+2*z3-z6)+16*mQ2*(z2*z5-2*z3*
1136 0 : z8-2*z3*z6+4*z3*z5+z8*z5)+8*s3*z3*z6-16*z3*z8*z5;
1137 0 : fm[3][4] = 16*mQ4*(-4*z1-z3+z8-z6)+4*mQ2*s3*(-2*z1-z3-
1138 0 : z6)+16*mQ2*(-2*z1S-3*z1*z2-2*z1*z3-3*z1*z8-2*z1*z6-3*
1139 0 : z1*z5-2*z2*z3-2*z6*z5)-8*s3*z3*z6+8*(-z1*z2*z8-2*z1*z2
1140 0 : *z5-z1*z8S-z1*z8*z5+z2S*z6-z2*z3*z5+z2*z8*z6-z2*z6*z5
1141 0 : +z3*z8*z5+z3*z5S);
1142 0 : fm[3][5] = 8*mQ4*(-2*z1+z10-z8-2*z6)+4*mQ2*(4*z1S+2*z1*
1143 0 : z2+4*z1*z3+4*z1*z4+6*z1*z10+4*z1*z9-2*z1*z8+8*z1*z7-2*z1*
1144 0 : z5+4*z2*z3+3*z2*z10+2*z2*z6-4*z3S-4*z3*z10+3*z3*z9+2*z3
1145 0 : *z7-3*z4*z8-2*z4*z6+z10*z5+3*z9*z6-z8*z7+4*z7*z6)+8*(-z1S
1146 0 : *z9+z1S*z8+z1*z2*z7-z1*z2*z6-z1*z3*z9+z1*z3*z8-z1*z3
1147 0 : *z5+z1*z4*z8+z1*z4*z5+z1*z10*z8-z1*z9*z6+z1*z8*z7+z2*z3*
1148 0 : z7-z2*z3*z6-z2*z4*z6-z2*z10*z6-z3S*z5-z3*z10*z5-z3*z7*
1149 0 : z5+z4*z6*z5);
1150 0 : fm[3][6] = 16*mQ6+4*mQ4*s3+16*mQ4*(-z1-z2+2*z3+2*z4+
1151 0 : z10-z9-z8-z7-z6+z5)+4*mQ2*s3*(z1+2*z3+2*z4+z10+z7+z6
1152 0 : )+8*mQ2*(4*z1*z3+4*z1*z4+4*z1*z10+4*z2*z3+4*z2*z4+4*z2*
1153 0 : z10-z2*z5+4*z3*z5+4*z4*z5+2*z10*z5-z9*z5-z8*z5)-(8*s3*
1154 0 : z1)*(z3+z4+z10)+16*z2*z5*(z3+z4+z10);
1155 0 : fm[3][7] = 8*mQ4*(3*z3+6*z4+3*z10+z9+2*z8+2*z7+z6)+2*mQ2*
1156 0 : s3*(z3+2*z4+2*z10-2*z7-z6)+4*mQ2*(4*z1*z3+8*z1*z4+4*
1157 0 : z1*z10+2*z1*z9-2*z1*z8+2*z2*z3+10*z2*z4+5*z2*z10+2*z2*z9+
1158 0 : z2*z8+2*z2*z7+4*z2*z6-7*z3*z9+2*z3*z8-8*z3*z7+4*z3*z6+4*
1159 0 : z3*z5+5*z4*z8+4*z4*z6+8*z4*z5+5*z10*z5-z9*z8-z9*z6+z9*z5+
1160 0 : z8S-z8*z7+2*z8*z6+2*z8*z5)+2*s3*(-z1*z10+z3*z7-2*z3*
1161 0 : z6+z4*z6)+4*(-z1*z2*z9-2*z1*z2*z8+z1*z9*z8-z1*z8S+z2S
1162 0 : *z7+2*z2S*z6+3*z2*z4*z5+2*z2*z10*z5-2*z2*z9*z6+z2*z8*z7
1163 0 : +z2*z8*z6-2*z3*z9*z5+z3*z8*z5+z4*z8*z5);
1164 0 : fm[3][8] = 8*mQ4*(3*z3+6*z4+3*z10+2*z9+z8+2*z7+z6)+2*mQ2*
1165 0 : s3*(3*z3+6*z4+3*z10-2*z7-z6)+4*mQ2*(4*z1*z3+8*z1*z4+
1166 0 : 4*z1*z10+4*z2*z3+8*z2*z4+4*z2*z10-8*z3*z9+4*z3*z8-8*z3*z7
1167 0 : +4*z3*z6+6*z3*z5+4*z4*z8+4*z4*z6+12*z4*z5+6*z10*z5+2*z9*
1168 0 : z5+z8*z5)+4*s3*(-z1*z3-2*z1*z4-z1*z10+2*z3*z7-z3*z6-z4
1169 0 : *z6)+8*z5*(z2*z3+2*z2*z4+z2*z10-2*z3*z9+z3*z8+z4*z8);
1170 0 : fm[4][4] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+2*z2+z3+z8+2*
1171 0 : z6+z5)+8*mQ2*s3*(-z1-z3+2*z6)+16*mQ2*(z2*z8+4*z2*
1172 0 : z6+z2*z5-2*z3*z6-2*z8*z6)+8*s3*z3*z6-16*z2*z8*z6;
1173 0 : fm[4][5] = 16*mQ6+8*mQ4*(-2*z1+z2-2*z3-2*z4-4*z10-z9+z8-4
1174 0 : *z7+2*z6+z5)+8*mQ2*(-2*z1*z2-2*z2*z3-2*z2*z10-2*z2*z7+
1175 0 : z2*z6+2*z3*z6-2*z4*z6+4*z10*z6-z9*z6-z8*z6)+16*z2*z6*(z3+
1176 : z10);
1177 0 : fm[4][6] = 16*mQ6-4*mQ4*s3+8*mQ4*(-2*z1+2*z2-4*z3-2*
1178 0 : z4-8*z10+z9+z8-4*z7-2*z6+2*z5)-(4*mQ2*s3)*(z1+z3+z10
1179 0 : +z7)+8*mQ2*(-2*z1*z2-2*z1*z10+z1*z9+z1*z8-2*z1*z5+z2S
1180 0 : -4*z2*z3-5*z2*z10+z2*z9-3*z2*z7-z2*z6+z2*z5+z3*z9+2*z3*z7
1181 0 : -3*z3*z5+z4*z8+2*z4*z6-z4*z5-5*z10*z5+z9*z8+z9*z6+z8*z7+
1182 0 : z8*z5-4*z7*z5+z5S)-(16*z2*z5)*(z1+z3+z10+z7);
1183 0 : fm[4][7] = 8*mQ4*(-z3-2*z4-3*z10-2*z9-z8-6*z7-3*z6)+2*mQ2
1184 0 : *s3*(z3+2*z4-3*z10-6*z7-3*z6)+4*mQ2*(-4*z1*z10-8*z1*
1185 0 : z7-4*z1*z6-6*z2*z10-2*z2*z9-z2*z8-12*z2*z7-6*z2*z6-4*z3*
1186 0 : z7-4*z3*z6+8*z4*z6-4*z10*z5+8*z9*z6-4*z8*z7-4*z8*z6-8*z7*
1187 0 : z5-4*z6*z5)+4*s3*(z1*z10+2*z1*z7+z1*z6+z3*z7+z3*z6-2*
1188 0 : z4*z6)+8*z2*(-z10*z5+2*z9*z6-z8*z7-z8*z6-2*z7*z5-z6*z5);
1189 0 : fm[4][8] = 8*mQ4*(-z3-2*z4-3*z10-z9-2*z8-6*z7-3*z6)+2*mQ2
1190 0 : *s3*(z3+2*z4-2*z10-2*z7-z6)+4*mQ2*(-4*z1*z10-2*z1*z9
1191 0 : +2*z1*z8-8*z1*z7-4*z1*z6-5*z2*z10-z2*z9-2*z2*z8-8*z2*z7-4
1192 0 : *z2*z6+z3*z9-2*z3*z8-4*z3*z7-4*z3*z6-4*z3*z5+z4*z8+8*z4*
1193 0 : z6-2*z4*z5-5*z10*z5+z9*z8+7*z9*z6-2*z9*z5-z8S-5*z8*z7-2
1194 0 : *z8*z6-z8*z5-10*z7*z5-2*z6*z5)+2*s3*(z1*z10-z3*z7+2*z3
1195 0 : *z6-z4*z6)+4*(-z1*z9*z8+z1*z9*z5+z1*z8S+2*z1*z8*z5-2*z2
1196 0 : *z10*z5+2*z2*z9*z6-z2*z8*z7-z2*z8*z6-3*z2*z7*z5+2*z3*z9*
1197 0 : z5-z3*z8*z5-2*z3*z5S-z4*z8*z5-z4*z5S);
1198 0 : fm[5][5] = 16*mQ6+16*mQ4*(-z1-z3+z4-z10-z7+z6)+16*mQ2*(
1199 0 : z3*z6+z4*z10+z4*z7+z4*z6+z10*z6)-16*z4*z10*z6;
1200 0 : fm[5][6] = 16*mQ6+8*mQ4*(-2*z1+z2-4*z3+2*z4-4*z10+z9-z8-2
1201 0 : *z7-2*z6+z5)+8*mQ2*(-2*z1*z5-2*z3*z5+4*z4*z10-z4*z9-z4*
1202 0 : z8+2*z4*z7-2*z4*z6+z4*z5-2*z10*z5-2*z7*z5)+16*z4*z5*(z10+
1203 : z7);
1204 0 : fm[5][7] = 8*mQ4*(-2*z3-z4-3*z10-2*z7-z6)+4*mQ2*(2*z1*z3+
1205 0 : 4*z1*z4+2*z1*z10+z1*z9-z1*z8+2*z1*z7+4*z1*z6-2*z2*z3-z2*
1206 0 : z4-3*z2*z10-2*z2*z7-z2*z6+6*z3S+6*z3*z4+6*z3*z10+z3*z9+
1207 0 : 3*z3*z8+2*z3*z7+4*z3*z6+2*z3*z5+6*z4*z10+2*z4*z8+4*z4*z7+
1208 0 : 2*z4*z6+z4*z5+3*z10*z9+3*z10*z8+6*z10*z7+6*z10*z6-z10*z5+
1209 0 : 2*z9*z7+2*z9*z6-z8*z6+6*z7S+6*z7*z6-2*z7*z5-z6*z5)+4*(-
1210 0 : z1S*z9+z1S*z8+2*z1*z2*z10+3*z1*z2*z7+3*z1*z2*z6-z1*z3
1211 0 : *z9-z1*z3*z8-z1*z3*z5-z1*z4*z8+z1*z4*z5-z1*z10*z9-z1*z10*
1212 0 : z8-z1*z9*z7+z1*z8*z7+z2*z3*z7+3*z2*z3*z6-z2*z4*z6+3*z2*
1213 0 : z10*z7+3*z2*z10*z6+3*z2*z7S+3*z2*z7*z6+z3S*z5+2*z3*z4
1214 0 : *z5+z3*z10*z5-z3*z7*z5+z4*z10*z5+z4*z7*z5);
1215 0 : fm[5][8] = 8*mQ4*(-2*z3-z4-3*z10-2*z7-z6)+4*mQ2*(2*z1*z3+
1216 0 : 4*z1*z4+2*z1*z10-z1*z9+z1*z8+2*z1*z7+4*z1*z6-2*z2*z3-z2*
1217 0 : z4-z2*z10+2*z2*z7+z2*z6+6*z3S+6*z3*z4+6*z3*z10+2*z3*z8+
1218 0 : 2*z3*z7+4*z3*z6-2*z3*z5+6*z4*z10-z4*z9+2*z4*z8+4*z4*z7+2*
1219 0 : z4*z6-z4*z5+3*z10*z9+3*z10*z8+6*z10*z7+6*z10*z6-3*z10*z5+
1220 0 : 3*z9*z7+2*z9*z6+z8*z7+6*z7S+6*z7*z6-2*z7*z5-z6*z5)+4*(
1221 0 : z1S*z9-z1S*z8-z1*z2*z7+z1*z2*z6+z1*z3*z9-z1*z3*z8+3*
1222 0 : z1*z3*z5+3*z1*z4*z5-z1*z10*z9-z1*z10*z8+2*z1*z10*z5-z1*z9
1223 0 : *z7-z1*z9*z6-z1*z8*z7-z2*z3*z7+z2*z3*z6+z2*z10*z7+z2*z10*
1224 0 : z6+z2*z7S+2*z2*z7*z6+3*z3S*z5+3*z3*z4*z5+3*z3*z10*z5+
1225 0 : z3*z7*z5+3*z4*z10*z5+3*z4*z7*z5-z4*z6*z5);
1226 0 : fm[6][6] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+z2+2*z4+z9+z7
1227 0 : +2*z5)+8*mQ2*s3*(-z1+2*z4-z7)+16*mQ2*(z2*z5-2*z4*
1228 0 : z9-2*z4*z7+4*z4*z5+z9*z5)+8*s3*z4*z7-16*z4*z9*z5;
1229 0 : fm[6][7] = 8*mQ4*(-6*z3-3*z4-3*z10-2*z9-z8-z7-2*z6)+2*mQ2
1230 0 : *s3*(-2*z3-z4-2*z10+z7+2*z6)+4*mQ2*(-8*z1*z3-4*z1*z4
1231 0 : -4*z1*z10+2*z1*z9-2*z1*z8-10*z2*z3-2*z2*z4-5*z2*z10-z2*z9
1232 0 : -2*z2*z8-4*z2*z7-2*z2*z6-5*z3*z9-4*z3*z7-8*z3*z5-2*z4*z9+
1233 0 : 7*z4*z8-4*z4*z7+8*z4*z6-4*z4*z5-5*z10*z5-z9S+z9*z8-2*z9
1234 0 : *z7+z9*z6-2*z9*z5+z8*z7-z8*z5)+2*s3*(z1*z10-z3*z7+2*z4
1235 0 : *z7-z4*z6)+4*(2*z1*z2*z9+z1*z2*z8+z1*z9S-z1*z9*z8-2*
1236 0 : z2S*z7-z2S*z6-3*z2*z3*z5-2*z2*z10*z5-z2*z9*z7-z2*z9*z6+
1237 0 : 2*z2*z8*z7-z3*z9*z5-z4*z9*z5+2*z4*z8*z5);
1238 0 : fm[6][8] = 8*mQ4*(-6*z3-3*z4-3*z10-z9-2*z8-z7-2*z6)+2*mQ2
1239 0 : *s3*(-6*z3-3*z4-3*z10+z7+2*z6)+4*mQ2*(-8*z1*z3-4*z1*
1240 0 : z4-4*z1*z10-8*z2*z3-4*z2*z4-4*z2*z10-4*z3*z9-4*z3*z7-12*
1241 0 : z3*z5-4*z4*z9+8*z4*z8-4*z4*z7+8*z4*z6-6*z4*z5-6*z10*z5-z9
1242 0 : *z5-2*z8*z5)+4*s3*(2*z1*z3+z1*z4+z1*z10+z3*z7+z4*z7-2*
1243 0 : z4*z6)+8*z5*(-2*z2*z3-z2*z4-z2*z10-z3*z9-z4*z9+2*z4*z8);
1244 0 : fm[7][7] = 72*mQ4*z10+18*mQ2*s3*z10+8*mQ2*(z1*z10+9*
1245 0 : z2*z10+7*z3*z7+2*z3*z6+2*z4*z7+7*z4*z6+z10*z5+2*z9*z7+7*
1246 0 : z9*z6+7*z8*z7+2*z8*z6)+2*s3*(-z1*z10-7*z3*z7-2*z3*z6-2
1247 0 : *z4*z7-7*z4*z6)+4*z2*(z10*z5+2*z9*z7+7*z9*z6+7*z8*z7+2*z8
1248 : *z6);
1249 0 : fm[7][8] = 72*mQ4*z10+2*mQ2*s3*z10+4*mQ2*(2*z1*z10+
1250 0 : 10*z2*z10+7*z3*z9+2*z3*z8+14*z3*z7+4*z3*z6+2*z4*z9+7*z4*
1251 0 : z8+4*z4*z7+14*z4*z6+10*z10*z5+z9S+7*z9*z8+2*z9*z7+7*z9*
1252 0 : z6+z8S+7*z8*z7+2*z8*z6)+2*s3*(7*z1*z10-7*z3*z7-2*z3*
1253 0 : z6-2*z4*z7-7*z4*z6)+2*(-2*z1*z9S-14*z1*z9*z8-2*z1*z8S
1254 0 : +2*z2*z10*z5+2*z2*z9*z7+7*z2*z9*z6+7*z2*z8*z7+2*z2*z8*z6+
1255 0 : 7*z3*z9*z5+2*z3*z8*z5+2*z4*z9*z5+7*z4*z8*z5);
1256 0 : fm[8][8] = 72*mQ4*z10+18*mQ2*s3*z10+8*mQ2*(z1*z10+z2
1257 0 : *z10+7*z3*z9+2*z3*z8+7*z3*z7+2*z3*z6+2*z4*z9+7*z4*z8+2*z4
1258 0 : *z7+7*z4*z6+9*z10*z5)+2*s3*(-z1*z10-7*z3*z7-2*z3*z6-2*
1259 0 : z4*z7-7*z4*z6)+4*z5*(z2*z10+7*z3*z9+2*z3*z8+2*z4*z9+7*z4*
1260 : z8);
1261 0 : double fm99 = -4*mQ4*z10-mQ2*s3*z10+4*mQ2*(-z1*z10-z2*z10+
1262 0 : z3*z7+z4*z6-z10*z5+z9*z6+z8*z7)+s3*(z1*z10-z3*z7-z4*z6
1263 0 : )+2*z2*(-z10*z5+z9*z6+z8*z7);
1264 0 : double fm910 = -4*mQ4*z10-mQ2*s3*z10+2*mQ2*(-2*z1*z10-2*z2*
1265 0 : z10+2*z3*z9+2*z3*z7+2*z4*z6-2*z10*z5+z9*z8+2*z8*z7)+s3
1266 0 : *(z1*z10-z3*z7-z4*z6)+2*(-z1*z9*z8-z2*z10*z5+z2*z8*z7+z3*
1267 : z9*z5);
1268 0 : double fmxx = -4*mQ4*z10-mQ2*s3*z10+2*mQ2*(-2*z1*z10-2*z2*
1269 0 : z10+2*z4*z8+2*z4*z6+2*z3*z7-2*z10*z5+z9*z8+2*z9*z6)+s3
1270 0 : *(z1*z10-z3*z7-z4*z6)+2*(-z1*z9*z8-z2*z10*z5+z2*z9*z6+z4*
1271 : z8*z5);
1272 0 : fm910 = 0.5*(fmxx+fm910);
1273 0 : double fm1010 = -4*mQ4*z10-mQ2*s3*z10+4*mQ2*(-z1*z10-z2*z10+
1274 0 : z3*z7+z4*z6-z10*z5+z9*z3+z8*z4)+s3*(z1*z10-z3*z7-z4*z6
1275 0 : )+2*z5*(-z10*z2+z9*z3+z8*z4);
1276 0 : fm[7][7] -= 2. * fm99;
1277 0 : fm[7][8] -= 2. * fm910;
1278 0 : fm[8][8] -= 2. * fm1010;
1279 :
1280 : // Propagators.
1281 0 : double ss1 = (pTemp[1] + pTemp[3]).m2Calc() - mQ2;
1282 0 : double ss2 = (pTemp[1] + pTemp[4]).m2Calc() - mQ2;
1283 0 : double ss3 = (pTemp[1] + pTemp[5]).m2Calc() - mQ2;
1284 0 : double ss4 = (pTemp[2] + pTemp[3]).m2Calc() - mQ2;
1285 0 : double ss5 = (pTemp[2] + pTemp[4]).m2Calc() - mQ2;
1286 0 : double ss6 = (pTemp[2] + pTemp[5]).m2Calc() - mQ2;
1287 0 : double ss7 = sH;
1288 :
1289 : // Propagator combinations.
1290 0 : double dz[9];
1291 0 : dz[1] = ss1 * ss6;
1292 0 : dz[2] = ss2 * ss6;
1293 0 : dz[3] = ss2 * ss4;
1294 0 : dz[4] = ss1 * ss5;
1295 0 : dz[5] = ss3 * ss5;
1296 0 : dz[6] = ss3 * ss4;
1297 0 : dz[7] = ss7 * ss1;
1298 0 : dz[8] = ss7 * ss4;
1299 :
1300 : // Colour factors.
1301 0 : double clr[9][9];
1302 0 : for (int i = 1; i < 4; ++i)
1303 0 : for (int j = 1; j < 4; ++j) {
1304 0 : clr[i][j] = 16. / 3.;
1305 0 : clr[i][j+3] = -2. / 3.;
1306 0 : clr[i+3][j] = -2. / 3.;
1307 0 : clr[i+3][j+3] = 16. / 3.;
1308 : }
1309 0 : for (int i = 1; i < 4; ++i)
1310 0 : for (int j = 1; j < 3; ++j) {
1311 0 : clr[i][j+6] = -6.;
1312 0 : clr[i+3][j+6] = 6.;
1313 0 : clr[j+6][i] = -6.;
1314 0 : clr[j+6][i+3] = 6.;
1315 : }
1316 0 : for (int i = 1; i < 3; ++i)
1317 0 : for (int j = 1; j < 3; ++j)
1318 0 : clr[i+6][j+6] = 12.;
1319 :
1320 : // Produce final result: matrix elements * colours * propagators.
1321 : double wtSum = 0.;
1322 0 : for (int i = 1; i < 9; ++i)
1323 0 : for (int j = i; j < 9; ++j) {
1324 0 : double fac = (j == i) ? 4. : 8.;
1325 0 : wtSum += fm[i][j] * fac * clr[i][j] / (dz[i] * dz[j]);
1326 : }
1327 0 : wtSum *= -1./256.;
1328 :
1329 : // Combine factors.
1330 0 : sigma = prefac * alpEM * pow2(alpS) * mQ2run * wtSum *pow2(coup2Q);
1331 :
1332 : // Secondary width for H, Q and Qbar (latter for top only).
1333 : // (H can be H0 SM or H1, H2, A3 from BSM).
1334 0 : sigma *= openFracTriplet;
1335 :
1336 0 : }
1337 :
1338 : //--------------------------------------------------------------------------
1339 :
1340 : // Select identity, colour and anticolour.
1341 :
1342 : void Sigma3gg2HQQbar::setIdColAcol() {
1343 :
1344 : // Pick out-flavours by relative CKM weights.
1345 0 : setId( id1, id2, idRes, idNew, -idNew);
1346 :
1347 : // Colour flow topologies.
1348 0 : if (rndmPtr->flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 0, 0, 3);
1349 0 : else setColAcol( 1, 2, 3, 1, 0, 0, 3, 0, 0, 2);
1350 :
1351 0 : }
1352 :
1353 : //--------------------------------------------------------------------------
1354 :
1355 : // Evaluate weight for decay angles.
1356 :
1357 : double Sigma3gg2HQQbar::weightDecay( Event& process, int iResBeg,
1358 : int iResEnd) {
1359 :
1360 : // Identity of mother of decaying reseonance(s).
1361 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
1362 :
1363 : // For Higgs decay hand over to standard routine.
1364 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
1365 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
1366 :
1367 : // For top decay hand over to standard routine.
1368 0 : if (idMother == 6)
1369 0 : return weightTopDecay( process, iResBeg, iResEnd);
1370 :
1371 : // Else done.
1372 0 : return 1.;
1373 :
1374 0 : }
1375 :
1376 : //==========================================================================
1377 :
1378 : // Sigma3qqbar2HQQbar class.
1379 : // Cross section for q qbar -> H0 Q Qbar (Q Qbar fusion of SM Higgs).
1380 : // REDUCE output and part of the rest courtesy Z. Kunszt,
1381 : // see Z. Kunszt, Nucl. Phys. B247 (1984) 339.
1382 :
1383 : //--------------------------------------------------------------------------
1384 :
1385 : // Initialize process.
1386 :
1387 : void Sigma3qqbar2HQQbar::initProc() {
1388 :
1389 : // Properties specific to Higgs state for the "q qbar -> H ttbar" process.
1390 : // (H can be H0 SM or H1, H2, A3 from BSM).
1391 :
1392 0 : if (higgsType == 0 && idNew == 6) {
1393 0 : nameSave = "q qbar -> H t tbar (SM)";
1394 0 : codeSave = 909;
1395 0 : idRes = 25;
1396 0 : coup2Q = 1.;
1397 0 : }
1398 0 : else if (higgsType == 1 && idNew == 6) {
1399 0 : nameSave = "q qbar -> h0(H1) t tbar";
1400 0 : codeSave = 1009;
1401 0 : idRes = 25;
1402 0 : coup2Q = settingsPtr->parm("HiggsH1:coup2u");
1403 0 : }
1404 0 : else if (higgsType == 2 && idNew == 6) {
1405 0 : nameSave = "q qbar -> H0(H2) t tbar";
1406 0 : codeSave = 1029;
1407 0 : idRes = 35;
1408 0 : coup2Q = settingsPtr->parm("HiggsH2:coup2u");
1409 0 : }
1410 0 : else if (higgsType == 3 && idNew == 6) {
1411 0 : nameSave = "q qbar -> A0(A3) t tbar";
1412 0 : codeSave = 1049;
1413 0 : idRes = 36;
1414 0 : coup2Q = settingsPtr->parm("HiggsA3:coup2u");
1415 0 : }
1416 :
1417 : // Properties specific to Higgs state for the "q qbar -> H b bbar" process.
1418 : // (H can be H0 SM or H1, H2, A3 from BSM).
1419 0 : if (higgsType == 0 && idNew == 5) {
1420 0 : nameSave = "q qbar -> H b bbar (SM)";
1421 0 : codeSave = 913;
1422 0 : idRes = 25;
1423 0 : coup2Q = 1.;
1424 0 : }
1425 0 : else if (higgsType == 1 && idNew == 5) {
1426 0 : nameSave = "q qbar -> h0(H1) b bbar";
1427 0 : codeSave = 1013;
1428 0 : idRes = 25;
1429 0 : coup2Q = settingsPtr->parm("HiggsH1:coup2d");
1430 0 : }
1431 0 : else if (higgsType == 2 && idNew == 5) {
1432 0 : nameSave = "q qbar -> H0(H2) b bbar";
1433 0 : codeSave = 1033;
1434 0 : idRes = 35;
1435 0 : coup2Q = settingsPtr->parm("HiggsH2:coup2d");
1436 0 : }
1437 0 : else if (higgsType == 3 && idNew == 5) {
1438 0 : nameSave = "q qbar -> A0(A3) b bbar";
1439 0 : codeSave = 1053;
1440 0 : idRes = 36;
1441 0 : coup2Q = settingsPtr->parm("HiggsA3:coup2d");
1442 0 : }
1443 :
1444 : // Common mass and coupling factors.
1445 0 : double mWS = pow2(particleDataPtr->m0(24));
1446 0 : prefac = (4. * M_PI / couplingsPtr->sin2thetaW()) * pow2(4. * M_PI)
1447 0 : * 0.25 / mWS;
1448 :
1449 : // Secondary open width fraction.
1450 0 : openFracTriplet = particleDataPtr->resOpenFrac(idRes, idNew, -idNew);
1451 :
1452 0 : }
1453 :
1454 : //--------------------------------------------------------------------------
1455 :
1456 : // Evaluate sigma(sHat), part independent of incoming flavour.
1457 :
1458 : void Sigma3qqbar2HQQbar::sigmaKin() {
1459 :
1460 : // Running mass of heavy quark.
1461 0 : double mQ2run = pow2( particleDataPtr->mRun(idNew, mH) );
1462 :
1463 : // Linear combination of p_Q and p_Qbar to ensure common mass.
1464 0 : double mQ2 = m4 * m5;
1465 : double epsi = 0.;
1466 0 : if (m4 != m5) {
1467 0 : double s45 = (p4cm + p5cm).m2Calc();
1468 0 : mQ2 = 0.5 * (s4 + s5) - 0.25 * pow2(s4 - s5) / s45;
1469 0 : epsi = 0.5 * (s5 - s4) / s45;
1470 0 : }
1471 :
1472 : // Set up kinematics: q(4) qbar(5) -> H(3) Q(1) Qbar(2) in outgoing sense.
1473 0 : Vec4 pTemp[6];
1474 0 : pTemp[4] = Vec4( 0., 0., -0.5* mH, -0.5* mH);
1475 0 : pTemp[5] = Vec4( 0., 0., 0.5* mH, -0.5* mH);
1476 0 : pTemp[1] = p4cm + epsi * (p4cm + p5cm);
1477 0 : pTemp[2] = p5cm - epsi * (p4cm + p5cm);
1478 0 : pTemp[3] = p3cm;
1479 :
1480 : // Four-product combinations.
1481 0 : double z1 = pTemp[1] * pTemp[2];
1482 0 : double z2 = pTemp[1] * pTemp[3];
1483 0 : double z3 = pTemp[1] * pTemp[4];
1484 0 : double z4 = pTemp[1] * pTemp[5];
1485 0 : double z5 = pTemp[2] * pTemp[3];
1486 0 : double z6 = pTemp[2] * pTemp[4];
1487 0 : double z7 = pTemp[2] * pTemp[5];
1488 0 : double z8 = pTemp[3] * pTemp[4];
1489 0 : double z9 = pTemp[3] * pTemp[5];
1490 0 : double z10 = pTemp[4] * pTemp[5];
1491 :
1492 : // Powers required as shorthand in matriz elements.
1493 0 : double mQ4 = mQ2 * mQ2;
1494 :
1495 : // Evaluate matrix elements for q + qbar -> Q + Qbar + H.
1496 : // (H can be H0 SM or H1, H2, A3 from BSM).
1497 0 : double a11 = -8.*mQ4*z10-2.*mQ2*s3*z10-(8.*mQ2)*(z2*z10+z3
1498 0 : *z7+z4*z6+z9*z6+z8*z7)+2.*s3*(z3*z7+z4*z6)-(4.*z2)*(z9
1499 0 : *z6+z8*z7);
1500 0 : double a12 = -8.*mQ4*z10+4.*mQ2*(-z2*z10-z3*z9-2.*z3*z7-z4*z8-
1501 0 : 2.*z4*z6-z10*z5-z9*z8-z9*z6-z8*z7)+2.*s3*(-z1*z10+z3*z7
1502 0 : +z4*z6)+2.*(2.*z1*z9*z8-z2*z9*z6-z2*z8*z7-z3*z9*z5-z4*z8*
1503 : z5);
1504 0 : double a22 = -8.*mQ4*z10-2.*mQ2*s3*z10-(8.*mQ2)*(z3*z9+z3*
1505 0 : z7+z4*z8+z4*z6+z10*z5)+2.*s3*(z3*z7+z4*z6)-(4.*z5)*(z3
1506 0 : *z9+z4*z8);
1507 :
1508 : // Propagators and propagator combinations.
1509 0 : double ss1 = (pTemp[1] + pTemp[3]).m2Calc() - mQ2;
1510 0 : double ss4 = (pTemp[2] + pTemp[3]).m2Calc() - mQ2;
1511 0 : double ss7 = sH;
1512 0 : double dz7 = ss7 * ss1;
1513 0 : double dz8 = ss7 * ss4;
1514 :
1515 : // Produce final result: matrix elements * propagators.
1516 0 : a11 /= (dz7 * dz7);
1517 0 : a12 /= (dz7 * dz8);
1518 0 : a22 /= (dz8 * dz8);
1519 0 : double wtSum = -(a11 + a22 + 2.*a12) * (8./9.);
1520 :
1521 : // Combine factors.
1522 0 : sigma = prefac * alpEM * pow2(alpS) * mQ2run * wtSum * pow2(coup2Q);
1523 :
1524 : // Secondary width for H, Q and Qbar (latter for top only).
1525 : // (H can be H0 SM or H1, H2, A3 from BSM).
1526 0 : sigma *= openFracTriplet;
1527 :
1528 0 : }
1529 :
1530 : //--------------------------------------------------------------------------
1531 :
1532 : // Select identity, colour and anticolour.
1533 :
1534 : void Sigma3qqbar2HQQbar::setIdColAcol() {
1535 :
1536 : // Pick out-flavours by relative CKM weights.
1537 0 : setId( id1, id2, idRes, idNew, -idNew);
1538 :
1539 : // Colour flow topologies.
1540 0 : if (id1 > 0) setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
1541 0 : else setColAcol( 0, 1, 2, 0, 0, 0, 2, 0, 0, 1);
1542 :
1543 0 : }
1544 :
1545 : //--------------------------------------------------------------------------
1546 :
1547 : // Evaluate weight for decay angles.
1548 :
1549 : double Sigma3qqbar2HQQbar::weightDecay( Event& process, int iResBeg,
1550 : int iResEnd) {
1551 :
1552 : // Identity of mother of decaying reseonance(s).
1553 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
1554 :
1555 : // For Higgs decay hand over to standard routine.
1556 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
1557 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
1558 :
1559 : // For top decay hand over to standard routine.
1560 0 : if (idMother == 6)
1561 0 : return weightTopDecay( process, iResBeg, iResEnd);
1562 :
1563 : // Else done.
1564 0 : return 1.;
1565 :
1566 0 : }
1567 :
1568 : //==========================================================================
1569 :
1570 : // Sigma2qg2Hq class.
1571 : // Cross section for q g -> H q.
1572 : // (H can be H0 SM or H1, H2, A3 from BSM).
1573 :
1574 : //--------------------------------------------------------------------------
1575 :
1576 : // Initialize process.
1577 :
1578 : void Sigma2qg2Hq::initProc() {
1579 :
1580 : // Properties specific to Higgs state for the "c g -> H c" process.
1581 : // (H can be H0 SM or H1, H2, A3 from BSM).
1582 0 : if (higgsType == 0 && idNew == 4) {
1583 0 : nameSave = "c g -> H c (SM)";
1584 0 : codeSave = 911;
1585 0 : idRes = 25;
1586 0 : }
1587 0 : else if (higgsType == 1 && idNew == 4) {
1588 0 : nameSave = "c g -> h0(H1) c";
1589 0 : codeSave = 1011;
1590 0 : idRes = 25;
1591 0 : }
1592 0 : else if (higgsType == 2 && idNew == 4) {
1593 0 : nameSave = "c g -> H0(H2) c";
1594 0 : codeSave = 1031;
1595 0 : idRes = 35;
1596 0 : }
1597 0 : else if (higgsType == 3 && idNew == 4) {
1598 0 : nameSave = "c g -> A0(A3) c";
1599 0 : codeSave = 1051;
1600 0 : idRes = 36;
1601 0 : }
1602 :
1603 : // Properties specific to Higgs state for the "b g -> H b" process.
1604 : // (H can be H0 SM or H1, H2, A3 from BSM).
1605 0 : if (higgsType == 0 && idNew == 5) {
1606 0 : nameSave = "b g -> H b (SM)";
1607 0 : codeSave = 911;
1608 0 : idRes = 25;
1609 0 : }
1610 0 : else if (higgsType == 1 && idNew == 5) {
1611 0 : nameSave = "b g -> h0(H1) b";
1612 0 : codeSave = 1011;
1613 0 : idRes = 25;
1614 0 : }
1615 0 : else if (higgsType == 2 && idNew == 5) {
1616 0 : nameSave = "b g -> H0(H2) b";
1617 0 : codeSave = 1031;
1618 0 : idRes = 35;
1619 0 : }
1620 0 : else if (higgsType == 3 && idNew == 5) {
1621 0 : nameSave = "b g -> A0(A3) b";
1622 0 : codeSave = 1051;
1623 0 : idRes = 36;
1624 0 : }
1625 :
1626 : // Standard parameters.
1627 0 : m2W = pow2( particleDataPtr->m0(24) );
1628 0 : thetaWRat = 1. / (24. * couplingsPtr->sin2thetaW());
1629 :
1630 : // Secondary open width fraction.
1631 0 : openFrac = particleDataPtr->resOpenFrac(idRes);
1632 :
1633 :
1634 0 : }
1635 :
1636 : //--------------------------------------------------------------------------
1637 :
1638 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1639 :
1640 : void Sigma2qg2Hq::sigmaKin() {
1641 :
1642 : // Running mass provides coupling.
1643 0 : double m2Run = pow2( particleDataPtr->mRun(idNew, mH) );
1644 :
1645 : // Cross section, including couplings and kinematics.
1646 0 : sigma = (M_PI / sH2) * alpS * alpEM * thetaWRat * (m2Run/m2W)
1647 0 : * (sH / (s4 - uH) + 2. * s4 * (s3 - uH) / pow2(s4 - uH)
1648 0 : + (s4 - uH) / sH - 2. * s4 / (s4 - uH)
1649 0 : + 2. * (s3 - uH) * (s3 - s4 - sH) / ((s4 - uH) * sH) );
1650 :
1651 : // Include secondary width for H0, H1, H2 or A3. Done.
1652 0 : sigma *= openFrac;
1653 :
1654 0 : }
1655 :
1656 : //--------------------------------------------------------------------------
1657 :
1658 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
1659 :
1660 : double Sigma2qg2Hq::sigmaHat() {
1661 :
1662 : // Check that specified flavour present.
1663 0 : if (abs(id1) != idNew && abs(id2) != idNew) return 0.;
1664 :
1665 : // Answer.
1666 0 : return sigma;
1667 :
1668 0 : }
1669 :
1670 :
1671 : //--------------------------------------------------------------------------
1672 :
1673 : // Select identity, colour and anticolour.
1674 :
1675 : void Sigma2qg2Hq::setIdColAcol() {
1676 :
1677 : // Flavour set up for q g -> H0 q.
1678 0 : int idq = (id2 == 21) ? id1 : id2;
1679 0 : setId( id1, id2, idRes, idq);
1680 :
1681 : // tH defined between f and f': must swap tHat <-> uHat if q g in.
1682 0 : swapTU = (id2 == 21);
1683 :
1684 : // Colour flow topologies. Swap when antiquarks.
1685 0 : if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
1686 0 : else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
1687 0 : if (idq < 0) swapColAcol();
1688 :
1689 0 : }
1690 :
1691 : //--------------------------------------------------------------------------
1692 :
1693 : // Evaluate weight for decay angles.
1694 :
1695 : double Sigma2qg2Hq::weightDecay( Event& process, int iResBeg,
1696 : int iResEnd) {
1697 :
1698 : // Identity of mother of decaying reseonance(s).
1699 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
1700 :
1701 : // For Higgs decay hand over to standard routine.
1702 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
1703 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
1704 :
1705 : // For top decay hand over to standard routine.
1706 0 : if (idMother == 6)
1707 0 : return weightTopDecay( process, iResBeg, iResEnd);
1708 :
1709 : // Else done.
1710 0 : return 1.;
1711 :
1712 0 : }
1713 :
1714 : //==========================================================================
1715 :
1716 : // Sigma2gg2Hglt class.
1717 : // Cross section for g g -> H g (H SM Higgs or BSM Higgs) via top loop.
1718 : // (H can be H0 SM or H1, H2, A3 from BSM).
1719 :
1720 : //--------------------------------------------------------------------------
1721 :
1722 : // Initialize process.
1723 :
1724 : void Sigma2gg2Hglt::initProc() {
1725 :
1726 : // Properties specific to Higgs state.
1727 0 : if (higgsType == 0) {
1728 0 : nameSave = "g g -> H g (SM; top loop)";
1729 0 : codeSave = 914;
1730 0 : idRes = 25;
1731 0 : }
1732 0 : else if (higgsType == 1) {
1733 0 : nameSave = "g g -> h0(H1) g (BSM; top loop)";
1734 0 : codeSave = 1014;
1735 0 : idRes = 25;
1736 0 : }
1737 0 : else if (higgsType == 2) {
1738 0 : nameSave = "g g -> H0(H2) g (BSM; top loop)";
1739 0 : codeSave = 1034;
1740 0 : idRes = 35;
1741 0 : }
1742 0 : else if (higgsType == 3) {
1743 0 : nameSave = "g g -> A0(A3) g (BSM; top loop)";
1744 0 : codeSave = 1054;
1745 0 : idRes = 36;
1746 0 : }
1747 :
1748 : // Normalization factor by g g -> H partial width.
1749 : // (H can be H0 SM or H1, H2, A3 from BSM).
1750 0 : double mHiggs = particleDataPtr->m0(idRes);
1751 0 : widHgg = particleDataPtr->resWidthChan(idRes, mHiggs, 21, 21);
1752 :
1753 : // Secondary open width fraction.
1754 0 : openFrac = particleDataPtr->resOpenFrac(idRes);
1755 :
1756 0 : }
1757 :
1758 : //--------------------------------------------------------------------------
1759 :
1760 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1761 :
1762 : void Sigma2gg2Hglt::sigmaKin() {
1763 :
1764 : // Evaluate cross section. Secondary width for H0, H1, H2 or A3.
1765 0 : sigma = (M_PI / sH2) * (3. / 16.) * alpS * (widHgg / m3)
1766 0 : * (sH2 * sH2 + tH2 * tH2 + uH2 * uH2 + pow4(s3))
1767 0 : / (sH * tH * uH * s3);
1768 0 : sigma *= openFrac;
1769 :
1770 0 : }
1771 :
1772 : //--------------------------------------------------------------------------
1773 :
1774 : // Select identity, colour and anticolour.
1775 :
1776 : void Sigma2gg2Hglt::setIdColAcol() {
1777 :
1778 : // Flavour set up for g g -> H g trivial.
1779 : // (H can be H0 SM or H1, H2, A3 from BSM).
1780 0 : setId( 21, 21, idRes, 21);
1781 :
1782 : // Colour flow topologies: random choice between two mirrors.
1783 0 : if (rndmPtr->flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
1784 0 : else setColAcol( 1, 2, 3, 1, 0, 0, 3, 2);
1785 :
1786 0 : }
1787 :
1788 : //--------------------------------------------------------------------------
1789 :
1790 : // Evaluate weight for decay angles.
1791 :
1792 : double Sigma2gg2Hglt::weightDecay( Event& process, int iResBeg,
1793 : int iResEnd) {
1794 :
1795 : // Identity of mother of decaying reseonance(s).
1796 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
1797 :
1798 : // For Higgs decay hand over to standard routine.
1799 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
1800 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
1801 :
1802 : // For top decay hand over to standard routine.
1803 0 : if (idMother == 6)
1804 0 : return weightTopDecay( process, iResBeg, iResEnd);
1805 :
1806 : // Else done.
1807 0 : return 1.;
1808 :
1809 0 : }
1810 :
1811 : //==========================================================================
1812 :
1813 : // Sigma2qg2Hqlt class.
1814 : // Cross section for q g -> H q (H SM or BSM Higgs) via top loop.
1815 : // (H can be H0 SM or H1, H2, A3 from BSM).
1816 :
1817 : //--------------------------------------------------------------------------
1818 :
1819 : // Initialize process.
1820 :
1821 : void Sigma2qg2Hqlt::initProc() {
1822 :
1823 : // Properties specific to Higgs state.
1824 0 : if (higgsType == 0) {
1825 0 : nameSave = "q g -> H q (SM; top loop)";
1826 0 : codeSave = 915;
1827 0 : idRes = 25;
1828 0 : }
1829 0 : else if (higgsType == 1) {
1830 0 : nameSave = "q g -> h0(H1) q (BSM; top loop)";
1831 0 : codeSave = 1015;
1832 0 : idRes = 25;
1833 0 : }
1834 0 : else if (higgsType == 2) {
1835 0 : nameSave = "q g -> H0(H2) q (BSM; top loop)";
1836 0 : codeSave = 1035;
1837 0 : idRes = 35;
1838 0 : }
1839 0 : else if (higgsType == 3) {
1840 0 : nameSave = "q g -> A0(A3) q (BSM; top loop)";
1841 0 : codeSave = 1055;
1842 0 : idRes = 36;
1843 0 : }
1844 :
1845 : // Normalization factor by g g -> H partial width.
1846 : // (H can be H0 SM or H1, H2, A3 from BSM).
1847 0 : double mHiggs = particleDataPtr->m0(idRes);
1848 0 : widHgg = particleDataPtr->resWidthChan(idRes, mHiggs, 21, 21);
1849 :
1850 : // Secondary open width fraction.
1851 0 : openFrac = particleDataPtr->resOpenFrac(idRes);
1852 :
1853 0 : }
1854 :
1855 : //--------------------------------------------------------------------------
1856 :
1857 : // Evaluate sigmaHat(sHat, part independent of incoming flavour).
1858 :
1859 : void Sigma2qg2Hqlt::sigmaKin() {
1860 :
1861 : // Evaluate cross section. Secondary width for H0, H1, H2 or A3.
1862 0 : sigma = (M_PI / sH2) * (1. / 12.) * alpS * (widHgg / m3)
1863 0 : * (sH2 + uH2) / (-tH * s3);
1864 0 : sigma *= openFrac;
1865 :
1866 0 : }
1867 :
1868 : //--------------------------------------------------------------------------
1869 :
1870 : // Select identity, colour and anticolour.
1871 :
1872 : void Sigma2qg2Hqlt::setIdColAcol() {
1873 :
1874 : // Flavour set up for q g -> H q.
1875 : // (H can be H0 SM or H1, H2, A3 from BSM).
1876 0 : int idq = (id2 == 21) ? id1 : id2;
1877 0 : setId( id1, id2, idRes, idq);
1878 :
1879 : // tH defined between f and f': must swap tHat <-> uHat if q g in.
1880 0 : swapTU = (id2 == 21);
1881 :
1882 : // Colour flow topologies. Swap when antiquarks.
1883 0 : if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
1884 0 : else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
1885 0 : if (idq < 0) swapColAcol();
1886 :
1887 0 : }
1888 :
1889 : //--------------------------------------------------------------------------
1890 :
1891 : // Evaluate weight for decay angles.
1892 :
1893 : double Sigma2qg2Hqlt::weightDecay( Event& process, int iResBeg,
1894 : int iResEnd) {
1895 :
1896 : // Identity of mother of decaying reseonance(s).
1897 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
1898 :
1899 : // For Higgs decay hand over to standard routine.
1900 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
1901 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
1902 :
1903 : // For top decay hand over to standard routine.
1904 0 : if (idMother == 6)
1905 0 : return weightTopDecay( process, iResBeg, iResEnd);
1906 :
1907 : // Else done.
1908 0 : return 1.;
1909 :
1910 0 : }
1911 :
1912 : //==========================================================================
1913 :
1914 : // Sigma2qqbar2Hglt class.
1915 : // Cross section for q qbar -> H g (H SM or BSM Higgs) via top loop.
1916 : // (H can be H0 SM or H1, H2, A3 from BSM).
1917 :
1918 : //--------------------------------------------------------------------------
1919 :
1920 : // Initialize process.
1921 :
1922 : void Sigma2qqbar2Hglt::initProc() {
1923 :
1924 : // Properties specific to Higgs state.
1925 0 : if (higgsType == 0) {
1926 0 : nameSave = "q qbar -> H g (SM; top loop)";
1927 0 : codeSave = 916;
1928 0 : idRes = 25;
1929 0 : }
1930 0 : else if (higgsType == 1) {
1931 0 : nameSave = "q qbar -> h0(H1) g (BSM; top loop)";
1932 0 : codeSave = 1016;
1933 0 : idRes = 25;
1934 0 : }
1935 0 : else if (higgsType == 2) {
1936 0 : nameSave = "q qbar -> H0(H2) g (BSM; top loop)";
1937 0 : codeSave = 1036;
1938 0 : idRes = 35;
1939 0 : }
1940 0 : else if (higgsType == 3) {
1941 0 : nameSave = "q qbar -> A0(A3) g (BSM; top loop)";
1942 0 : codeSave = 1056;
1943 0 : idRes = 36;
1944 0 : }
1945 :
1946 : // Normalization factor by g g -> H partial width.
1947 : // (H can be H0 SM or H1, H2, A3 from BSM).
1948 0 : double mHiggs = particleDataPtr->m0(idRes);
1949 0 : widHgg = particleDataPtr->resWidthChan(idRes, mHiggs, 21, 21);
1950 :
1951 : // Secondary open width fraction.
1952 0 : openFrac = particleDataPtr->resOpenFrac(idRes);
1953 :
1954 :
1955 0 : }
1956 :
1957 : //--------------------------------------------------------------------------
1958 :
1959 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1960 :
1961 : void Sigma2qqbar2Hglt::sigmaKin() {
1962 :
1963 : // Evaluate cross section. Secondary width for H0, H1, H2 or A3.
1964 0 : sigma = (M_PI / sH2) * (2. / 9.) * alpS * (widHgg / m3)
1965 0 : * (tH2 + uH2) / (sH * s3);
1966 0 : sigma *= openFrac;
1967 :
1968 0 : }
1969 :
1970 : //--------------------------------------------------------------------------
1971 :
1972 : // Select identity, colour and anticolour.
1973 :
1974 : void Sigma2qqbar2Hglt::setIdColAcol() {
1975 :
1976 : // Flavours trivial.
1977 0 : setId( id1, id2, idRes, 21);
1978 :
1979 : // Colour flow topologies. Swap when antiquarks.
1980 0 : setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
1981 0 : if (id1 < 0) swapColAcol();
1982 :
1983 0 : }
1984 :
1985 : //--------------------------------------------------------------------------
1986 :
1987 : // Evaluate weight for decay angles.
1988 :
1989 : double Sigma2qqbar2Hglt::weightDecay( Event& process, int iResBeg,
1990 : int iResEnd) {
1991 :
1992 : // Identity of mother of decaying reseonance(s).
1993 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
1994 :
1995 : // For Higgs decay hand over to standard routine.
1996 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
1997 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
1998 :
1999 : // For top decay hand over to standard routine.
2000 0 : if (idMother == 6)
2001 0 : return weightTopDecay( process, iResBeg, iResEnd);
2002 :
2003 : // Else done.
2004 0 : return 1.;
2005 :
2006 0 : }
2007 :
2008 :
2009 : //==========================================================================
2010 :
2011 : // Sigma1ffbar2Hchg class.
2012 : // Cross section for f fbar -> H+- (f is quark or lepton).
2013 :
2014 : //--------------------------------------------------------------------------
2015 :
2016 : // Initialize process.
2017 :
2018 : void Sigma1ffbar2Hchg::initProc() {
2019 :
2020 : // Find pointer to H+-.
2021 0 : HResPtr = particleDataPtr->particleDataEntryPtr(37);
2022 :
2023 : // Store H+- mass and width for propagator.
2024 0 : mRes = HResPtr->m0();
2025 0 : GammaRes = HResPtr->mWidth();
2026 0 : m2Res = mRes*mRes;
2027 0 : GamMRat = GammaRes / mRes;
2028 :
2029 : // Couplings.
2030 0 : m2W = pow2(particleDataPtr->m0(24));
2031 0 : thetaWRat = 1. / (8. * couplingsPtr->sin2thetaW());
2032 0 : tan2Beta = pow2(settingsPtr->parm("HiggsHchg:tanBeta"));
2033 :
2034 0 : }
2035 :
2036 : //--------------------------------------------------------------------------
2037 :
2038 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2039 :
2040 : void Sigma1ffbar2Hchg::sigmaKin() {
2041 :
2042 : // Set up Breit-Wigner. Width out only includes open channels.
2043 0 : sigBW = 4. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
2044 0 : widthOutPos = HResPtr->resWidthOpen( 37, mH);
2045 0 : widthOutNeg = HResPtr->resWidthOpen(-37, mH);
2046 :
2047 0 : }
2048 :
2049 : //--------------------------------------------------------------------------
2050 :
2051 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2052 :
2053 : double Sigma1ffbar2Hchg::sigmaHat() {
2054 :
2055 : // Only allow generation-diagonal states.
2056 0 : int id1Abs = abs(id1);
2057 0 : int id2Abs = abs(id2);
2058 0 : int idUp = max(id1Abs, id2Abs);
2059 0 : int idDn = min(id1Abs, id2Abs);
2060 0 : if (idUp%2 != 0 || idUp - idDn != 1) return 0.;
2061 :
2062 : // Calculate mass-dependent incoming width. Total cross section.
2063 0 : double m2RunUp = pow2(particleDataPtr->mRun(idUp, mH));
2064 0 : double m2RunDn = pow2(particleDataPtr->mRun(idDn, mH));
2065 0 : double widthIn = alpEM * thetaWRat * (mH/m2W)
2066 0 : * (m2RunDn * tan2Beta + m2RunUp / tan2Beta);
2067 0 : int idUpChg = (id1Abs%2 == 0) ? id1 : id2;
2068 0 : double sigma = (idUpChg > 0) ? widthIn * sigBW * widthOutPos
2069 0 : : widthIn * sigBW * widthOutNeg;
2070 :
2071 : // Colour factor. Answer.
2072 0 : if (idUp < 9) sigma /= 3.;
2073 : return sigma;
2074 :
2075 0 : }
2076 :
2077 : //--------------------------------------------------------------------------
2078 :
2079 : // Select identity, colour and anticolour.
2080 :
2081 : void Sigma1ffbar2Hchg::setIdColAcol() {
2082 :
2083 : // Charge of Higgs. Fill flavours.
2084 0 : int idUpChg = (abs(id1)%2 == 0) ? id1 : id2;
2085 0 : int idHchg = (idUpChg > 0) ? 37 : -37;
2086 0 : setId( id1, id2, idHchg);
2087 :
2088 : // Colour flow topologies. Swap when antiquarks.
2089 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2090 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
2091 0 : if (id1 < 0) swapColAcol();
2092 :
2093 0 : }
2094 :
2095 : //--------------------------------------------------------------------------
2096 :
2097 : // Evaluate weight for decay angles.
2098 :
2099 : double Sigma1ffbar2Hchg::weightDecay( Event& process, int iResBeg,
2100 : int iResEnd) {
2101 :
2102 : // Identity of mother of decaying reseonance(s).
2103 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
2104 :
2105 : // For Higgs decay hand over to standard routine.
2106 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
2107 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
2108 :
2109 : // For top decay hand over to standard routine.
2110 0 : if (idMother == 6)
2111 0 : return weightTopDecay( process, iResBeg, iResEnd);
2112 :
2113 : // Else done.
2114 0 : return 1.;
2115 :
2116 0 : }
2117 :
2118 : //==========================================================================
2119 :
2120 : // Sigma2qg2Hq class.
2121 : // Cross section for q g -> H+- q'.
2122 :
2123 : //--------------------------------------------------------------------------
2124 :
2125 : // Initialize process.
2126 :
2127 : void Sigma2qg2Hchgq::initProc() {
2128 :
2129 : // Standard parameters.
2130 0 : m2W = pow2( particleDataPtr->m0(24) );
2131 0 : thetaWRat = 1. / (24. * couplingsPtr->sin2thetaW());
2132 0 : tan2Beta = pow2(settingsPtr->parm("HiggsHchg:tanBeta"));
2133 :
2134 : // Incoming flavour within same doublet. Uptype and downtype flavours.
2135 0 : idOld = (idNew%2 == 0) ? idNew - 1 : idNew + 1;
2136 0 : idUp = max(idOld, idNew);
2137 0 : idDn = min(idOld, idNew);
2138 :
2139 : // Secondary open width fraction.
2140 0 : openFracPos = (idOld%2 == 0) ? particleDataPtr->resOpenFrac( 37, idNew)
2141 0 : : particleDataPtr->resOpenFrac(-37, idNew);
2142 0 : openFracNeg = (idOld%2 == 0) ? particleDataPtr->resOpenFrac(-37, -idNew)
2143 0 : : particleDataPtr->resOpenFrac( 37, -idNew);
2144 :
2145 0 : }
2146 :
2147 : //--------------------------------------------------------------------------
2148 :
2149 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2150 :
2151 : void Sigma2qg2Hchgq::sigmaKin() {
2152 :
2153 : // Running masses provides coupling.
2154 0 : double m2RunUp = pow2(particleDataPtr->mRun(idUp, mH));
2155 0 : double m2RunDn = pow2(particleDataPtr->mRun(idDn, mH));
2156 :
2157 : // Cross section, including couplings and kinematics.
2158 0 : sigma = (M_PI / sH2) * alpS * alpEM * thetaWRat
2159 0 : * (m2RunDn * tan2Beta + m2RunUp / tan2Beta) / m2W
2160 0 : * (sH / (s4 - uH) + 2. * s4 * (s3 - uH) / pow2(s4 - uH)
2161 0 : + (s4 - uH) / sH - 2. * s4 / (s4 - uH)
2162 0 : + 2. * (s3 - uH) * (s3 - s4 - sH) / ((s4 - uH) * sH) );
2163 :
2164 0 : }
2165 :
2166 : //--------------------------------------------------------------------------
2167 :
2168 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2169 :
2170 : double Sigma2qg2Hchgq::sigmaHat() {
2171 :
2172 : // Check that specified flavour present.
2173 0 : if (abs(id1) != idOld && abs(id2) != idOld) return 0.;
2174 :
2175 : // Answer.
2176 0 : return (id1 == idOld || id2 == idOld) ? sigma * openFracPos
2177 0 : : sigma * openFracNeg;
2178 :
2179 0 : }
2180 :
2181 : //--------------------------------------------------------------------------
2182 :
2183 : // Select identity, colour and anticolour.
2184 :
2185 : void Sigma2qg2Hchgq::setIdColAcol() {
2186 :
2187 : // Flavour set up for q g -> H+- q'.
2188 0 : int idq = (id2 == 21) ? id1 : id2;
2189 0 : id3 = ( (idq > 0 && idOld%2 == 0) || (idq < 0 && idOld%2 != 0) )
2190 : ? 37 : -37;
2191 0 : id4 = (idq > 0) ? idNew : -idNew;
2192 0 : setId( id1, id2, id3, id4);
2193 :
2194 : // tH defined between f and f': must swap tHat <-> uHat if q g in.
2195 0 : swapTU = (id2 == 21);
2196 :
2197 : // Colour flow topologies. Swap when antiquarks.
2198 0 : if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2199 0 : else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2200 0 : if (idq < 0) swapColAcol();
2201 :
2202 0 : }
2203 :
2204 : //--------------------------------------------------------------------------
2205 :
2206 : // Evaluate weight for decay angles.
2207 :
2208 : double Sigma2qg2Hchgq::weightDecay( Event& process, int iResBeg,
2209 : int iResEnd) {
2210 :
2211 : // Identity of mother of decaying reseonance(s).
2212 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
2213 :
2214 : // For Higgs decay hand over to standard routine.
2215 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
2216 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
2217 :
2218 : // For top decay hand over to standard routine.
2219 0 : if (idMother == 6)
2220 0 : return weightTopDecay( process, iResBeg, iResEnd);
2221 :
2222 : // Else done.
2223 0 : return 1.;
2224 :
2225 0 : }
2226 :
2227 : //==========================================================================
2228 :
2229 : // Sigma2ffbar2A3H12 class.
2230 : // Cross section for f fbar -> A0(H_3) h0(H_1) or A0(H_3) H0(H_2).
2231 :
2232 : //--------------------------------------------------------------------------
2233 :
2234 : // Initialize process.
2235 :
2236 : void Sigma2ffbar2A3H12::initProc() {
2237 :
2238 : // Set up whether h0(H_1) or H0(H_2).
2239 0 : higgs12 = (higgsType == 1) ? 25 : 35;
2240 0 : codeSave = (higgsType == 1) ? 1081 : 1082;
2241 0 : nameSave = (higgsType == 1) ? "f fbar -> A0(H3) h0(H1)"
2242 : : "f fbar -> A0(H3) H0(H2)";
2243 0 : coupZA3H12 = (higgsType == 1) ? settingsPtr->parm("HiggsA3:coup2H1Z")
2244 0 : : settingsPtr->parm("HiggsA3:coup2H2Z");
2245 :
2246 : // Standard parameters.
2247 0 : double mZ = particleDataPtr->m0(23);
2248 0 : double GammaZ = particleDataPtr->mWidth(23);
2249 0 : m2Z = mZ * mZ;
2250 0 : mGammaZ = mZ * GammaZ;
2251 0 : thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW()
2252 0 : * couplingsPtr->cos2thetaW());
2253 :
2254 : // Secondary open width fraction.
2255 0 : openFrac = particleDataPtr->resOpenFrac(36, higgs12);
2256 :
2257 0 : }
2258 :
2259 : //--------------------------------------------------------------------------
2260 :
2261 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2262 :
2263 : void Sigma2ffbar2A3H12::sigmaKin() {
2264 :
2265 : // Common kinematics factora.
2266 0 : sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat * coupZA3H12)
2267 0 : * (uH * tH - s3 * s4) / ( pow2(sH - m2Z) + pow2(mGammaZ) );
2268 :
2269 0 : }
2270 :
2271 : //--------------------------------------------------------------------------
2272 :
2273 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2274 :
2275 : double Sigma2ffbar2A3H12::sigmaHat() {
2276 :
2277 : // Couplings for incoming flavour.
2278 0 : int idAbs = abs(id1);
2279 0 : double lIn = couplingsPtr->lf(idAbs);
2280 0 : double rIn = couplingsPtr->rf(idAbs);
2281 :
2282 : // Combine to total cross section. Colour factor.
2283 0 : double sigma = (pow2(lIn) + pow2(rIn)) * sigma0 * openFrac;
2284 0 : if (idAbs < 9) sigma /= 3.;
2285 0 : return sigma;
2286 :
2287 0 : }
2288 :
2289 : //--------------------------------------------------------------------------
2290 :
2291 : // Select identity, colour and anticolour.
2292 :
2293 : void Sigma2ffbar2A3H12::setIdColAcol() {
2294 :
2295 : // Flavours trivial
2296 0 : setId( id1, id2, 36, higgs12);
2297 :
2298 : // Colour flow topologies. Swap when antiquarks.
2299 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2300 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
2301 0 : if (id1 < 0) swapColAcol();
2302 :
2303 0 : }
2304 :
2305 : //--------------------------------------------------------------------------
2306 :
2307 : // Evaluate weight for decay angles.
2308 :
2309 : double Sigma2ffbar2A3H12::weightDecay( Event& process, int iResBeg,
2310 : int iResEnd) {
2311 :
2312 : // Identity of mother of decaying reseonance(s).
2313 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
2314 :
2315 : // For Higgs decay hand over to standard routine.
2316 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
2317 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
2318 :
2319 : // For top decay hand over to standard routine.
2320 0 : if (idMother == 6)
2321 0 : return weightTopDecay( process, iResBeg, iResEnd);
2322 :
2323 : // Else done.
2324 0 : return 1.;
2325 :
2326 0 : }
2327 :
2328 : //==========================================================================
2329 :
2330 : // Sigma2ffbar2HchgH12 class.
2331 : // Cross section for f fbar -> H+- h0(H_1) or H+- H0(H_2).
2332 :
2333 : //--------------------------------------------------------------------------
2334 :
2335 : // Initialize process.
2336 :
2337 : void Sigma2ffbar2HchgH12::initProc() {
2338 :
2339 : // Set up whether h0(H_1) or H0(H_2).
2340 0 : higgs12 = (higgsType == 1) ? 25 : 35;
2341 0 : codeSave = (higgsType == 1) ? 1083 : 1084;
2342 0 : nameSave = (higgsType == 1) ? "f fbar' -> H+- h0(H1)"
2343 : : "f fbar' -> H+- H0(H2)";
2344 0 : coupWHchgH12 = (higgsType == 1) ? settingsPtr->parm("HiggsHchg:coup2H1W")
2345 0 : : settingsPtr->parm("HiggsHchg:coup2H2W");
2346 :
2347 : // Standard parameters.
2348 0 : double mW = particleDataPtr->m0(24);
2349 0 : double GammaW = particleDataPtr->mWidth(24);
2350 0 : m2W = mW * mW;
2351 0 : mGammaW = mW * GammaW;
2352 0 : thetaWRat = 1. / (2. * couplingsPtr->sin2thetaW());
2353 :
2354 : // Secondary open width fraction.
2355 0 : openFracPos = particleDataPtr->resOpenFrac( 37, higgs12);
2356 0 : openFracNeg = particleDataPtr->resOpenFrac(-37, higgs12);
2357 :
2358 0 : }
2359 :
2360 : //--------------------------------------------------------------------------
2361 :
2362 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2363 :
2364 : void Sigma2ffbar2HchgH12::sigmaKin() {
2365 :
2366 : // Common kinematics factora.
2367 0 : sigma0 = 0.5 * (M_PI / sH2) * pow2(alpEM * thetaWRat * coupWHchgH12)
2368 0 : * (uH * tH - s3 * s4) / ( pow2(sH - m2W) + pow2(mGammaW) );
2369 :
2370 0 : }
2371 :
2372 : //--------------------------------------------------------------------------
2373 :
2374 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2375 :
2376 : double Sigma2ffbar2HchgH12::sigmaHat() {
2377 :
2378 : // Combine to total cross section. CKM and colour factor.
2379 0 : int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2380 0 : double sigma = (idUp > 0) ? sigma0 * openFracPos : sigma0 * openFracNeg;
2381 0 : if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
2382 0 : return sigma;
2383 :
2384 : }
2385 :
2386 : //--------------------------------------------------------------------------
2387 :
2388 : // Select identity, colour and anticolour.
2389 :
2390 : void Sigma2ffbar2HchgH12::setIdColAcol() {
2391 :
2392 : // Charge of Higgs. Fill flavours.
2393 0 : int idUpChg = (abs(id1)%2 == 0) ? id1 : id2;
2394 0 : int idHchg = (idUpChg > 0) ? 37 : -37;
2395 0 : setId( id1, id2, idHchg, higgs12);
2396 :
2397 : // Colour flow topologies. Swap when antiquarks.
2398 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2399 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
2400 0 : if (id1 < 0) swapColAcol();
2401 :
2402 0 : }
2403 :
2404 : //--------------------------------------------------------------------------
2405 :
2406 : // Evaluate weight for decay angles.
2407 :
2408 : double Sigma2ffbar2HchgH12::weightDecay( Event& process, int iResBeg,
2409 : int iResEnd) {
2410 :
2411 : // Identity of mother of decaying reseonance(s).
2412 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
2413 :
2414 : // For Higgs decay hand over to standard routine.
2415 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
2416 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
2417 :
2418 : // For top decay hand over to standard routine.
2419 0 : if (idMother == 6)
2420 0 : return weightTopDecay( process, iResBeg, iResEnd);
2421 :
2422 : // Else done.
2423 0 : return 1.;
2424 :
2425 0 : }
2426 :
2427 : //==========================================================================
2428 :
2429 : // Sigma2ffbar2HposHneg class.
2430 : // Cross section for q g -> H+- q'.
2431 :
2432 : //--------------------------------------------------------------------------
2433 :
2434 : // Initialize process.
2435 :
2436 : void Sigma2ffbar2HposHneg::initProc() {
2437 :
2438 : // Standard parameters.
2439 0 : double mZ = particleDataPtr->m0(23);
2440 0 : double GammaZ = particleDataPtr->mWidth(23);
2441 0 : m2Z = mZ * mZ;
2442 0 : mGammaZ = mZ * GammaZ;
2443 0 : thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW()
2444 0 : * couplingsPtr->cos2thetaW());
2445 :
2446 : // Charged Higgs coupling to gamma and Z0.
2447 0 : eH = -1.;
2448 0 : lH = -1. + 2. * couplingsPtr->sin2thetaW();
2449 :
2450 : // Secondary open width fraction.
2451 0 : openFrac = particleDataPtr->resOpenFrac(37, -37);
2452 :
2453 0 : }
2454 :
2455 : //--------------------------------------------------------------------------
2456 :
2457 : // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2458 :
2459 : void Sigma2ffbar2HposHneg::sigmaKin() {
2460 :
2461 : // Common kinematics factora.
2462 0 : double preFac = M_PI * pow2(alpEM) * ((uH * tH - s3 * s4) / sH2);
2463 0 : double propZ = 1. / ( pow2(sH - m2Z) + pow2(mGammaZ) );
2464 :
2465 : // Separate parts for gamma*, interference and Z0.
2466 0 : gamSig = preFac * 2. * pow2(eH) / sH2;
2467 0 : intSig = preFac * 2. * eH * lH * thetaWRat * propZ * (sH - m2Z) / sH;
2468 0 : resSig = preFac * pow2(lH * thetaWRat) * propZ;
2469 :
2470 0 : }
2471 :
2472 : //--------------------------------------------------------------------------
2473 :
2474 : // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2475 :
2476 : double Sigma2ffbar2HposHneg::sigmaHat() {
2477 :
2478 : // Couplings for incoming flavour.
2479 0 : int idAbs = abs(id1);
2480 0 : double eIn = couplingsPtr->ef(idAbs);
2481 0 : double lIn = couplingsPtr->lf(idAbs);
2482 0 : double rIn = couplingsPtr->rf(idAbs);
2483 :
2484 : // Combine to total cross section. Colour factor.
2485 0 : double sigma = (pow2(eIn) * gamSig + eIn * (lIn + rIn) * intSig
2486 0 : + (pow2(lIn) + pow2(rIn)) * resSig) * openFrac;
2487 0 : if (idAbs < 9) sigma /= 3.;
2488 0 : return sigma;
2489 :
2490 0 : }
2491 :
2492 : //--------------------------------------------------------------------------
2493 :
2494 : // Select identity, colour and anticolour.
2495 :
2496 : void Sigma2ffbar2HposHneg::setIdColAcol() {
2497 :
2498 : // Flavours trivial
2499 0 : setId( id1, id2, 37, -37);
2500 :
2501 : // Colour flow topologies. Swap when antiquarks.
2502 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2503 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
2504 0 : if (id1 < 0) swapColAcol();
2505 :
2506 0 : }
2507 :
2508 : //--------------------------------------------------------------------------
2509 :
2510 : // Evaluate weight for decay angles.
2511 :
2512 : double Sigma2ffbar2HposHneg::weightDecay( Event& process, int iResBeg,
2513 : int iResEnd) {
2514 :
2515 : // Identity of mother of decaying reseonance(s).
2516 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
2517 :
2518 : // For Higgs decay hand over to standard routine.
2519 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
2520 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
2521 :
2522 : // For top decay hand over to standard routine.
2523 0 : if (idMother == 6)
2524 0 : return weightTopDecay( process, iResBeg, iResEnd);
2525 :
2526 : // Else done.
2527 0 : return 1.;
2528 :
2529 0 : }
2530 :
2531 : //==========================================================================
2532 :
2533 : } // end namespace Pythia8
|
