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