StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaHiggs.cc
1 // SigmaHiggs.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 Torbjorn Sjostrand.
3 // Part of code written by Marc Montull, CERN summer student 2007.
4 // PYTHIA is licenced under the GNU GPL v2 or later, 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  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 = particleDataPtr->particleDataEntryPtr(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 
63 void 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 
78 double 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 
94 void 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 
110 double 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 
138 void 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 = particleDataPtr->particleDataEntryPtr(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 
177 void 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 
198 void 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 
212 double 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 
241 void 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 = particleDataPtr->particleDataEntryPtr(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 
280 void 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 
301 void 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 
315 double 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 
344 void 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 = settingsPtr->parm("HiggsH1:coup2Z");
358  }
359  else if (higgsType == 2) {
360  nameSave = "f fbar -> H0(H2) Z0";
361  codeSave = 1024;
362  idRes = 35;
363  coup2Z = settingsPtr->parm("HiggsH2:coup2Z");
364  }
365  else if (higgsType == 3) {
366  nameSave = "f fbar -> A0(A3) ZO";
367  codeSave = 1044;
368  idRes = 36;
369  coup2Z = settingsPtr->parm("HiggsA3:coup2Z");
370  }
371 
372  // Store Z0 mass and width for propagator. Common coupling factor.
373  mZ = particleDataPtr->m0(23);
374  widZ = particleDataPtr->mWidth(23);
375  mZS = mZ*mZ;
376  mwZS = pow2(mZ * widZ);
377  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
378  * couplingsPtr->cos2thetaW());
379 
380  // Secondary open width fraction.
381  openFracPair = particleDataPtr->resOpenFrac(idRes, 23);
382 
383 }
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  sigma0 = (M_PI / sH2) * 8. * pow2(alpEM * thetaWRat * coup2Z)
393  * (tH * uH - s3 * s4 + 2. * sH * s4) / (pow2(sH - mZS) + mwZS);
394 
395 }
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  int idAbs = abs(id1);
405  double sigma = sigma0 * couplingsPtr->vf2af2(idAbs);
406  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  sigma *= openFracPair;
410 
411  // Answer.
412  return sigma;
413 
414 }
415 
416 //--------------------------------------------------------------------------
417 
418 // Select identity, colour and anticolour.
419 
420 void Sigma2ffbar2HZ::setIdColAcol() {
421 
422  // Flavours trivial.
423  setId( id1, id2, idRes, 23);
424 
425  // Colour flow topologies. Swap when antiquarks.
426  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
427  else setColAcol( 0, 0, 0, 0, 0, 0);
428  if (id1 < 0) swapColAcol();
429 
430 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
441 
442  // For Higgs decay hand over to standard routine.
443  if (idMother == 25 || idMother == 35 || idMother == 36)
444  return weightHiggsDecay( process, iResBeg, iResEnd);
445 
446  // For top decay hand over to standard routine.
447  if (idMother == 6)
448  return weightTopDecay( process, iResBeg, iResEnd);
449 
450  // If not decay of Z0 created along with Higgs then done.
451  if (iResBeg != 5 || iResEnd != 6) return 1.;
452 
453  // Order so that fbar(1) f(2) -> H() f'(3) fbar'(4).
454  int i1 = (process[3].id() < 0) ? 3 : 4;
455  int i2 = 7 - i1;
456  int i3 = process[6].daughter1();
457  int i4 = process[6].daughter2();
458  if (process[i3].id() < 0) swap( i3, i4);
459 
460  // Find left- and righthanded couplings of fermion pairs.
461  int idAbs = process[i1].idAbs();
462  double liS = pow2( couplingsPtr->lf(idAbs) );
463  double riS = pow2( couplingsPtr->rf(idAbs) );
464  idAbs = process[i3].idAbs();
465  double lfS = pow2( couplingsPtr->lf(idAbs) );
466  double rfS = pow2( couplingsPtr->rf(idAbs) );
467 
468  // Evaluate relevant four-products.
469  double pp13 = process[i1].p() * process[i3].p();
470  double pp14 = process[i1].p() * process[i4].p();
471  double pp23 = process[i2].p() * process[i3].p();
472  double pp24 = process[i2].p() * process[i4].p();
473 
474  // Weight and maximum.
475  double wt = (liS * lfS + riS * rfS) * pp13 * pp24
476  + (liS * rfS + riS * lfS) * pp14 * pp23;
477  double wtMax = (liS + riS) * (lfS + rfS) * (pp13 + pp14) * (pp23 + pp24);
478 
479  // Done.
480  return wt / wtMax;
481 
482 }
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  if (higgsType == 0) {
498  nameSave = "f fbar -> H0 W+- (SM)";
499  codeSave = 905;
500  idRes = 25;
501  coup2W = 1.;
502  }
503  else if (higgsType == 1) {
504  nameSave = "f fbar -> h0(H1) W+-";
505  codeSave = 1005;
506  idRes = 25;
507  coup2W = settingsPtr->parm("HiggsH1:coup2W");
508  }
509  else if (higgsType == 2) {
510  nameSave = "f fbar -> H0(H2) W+-";
511  codeSave = 1025;
512  idRes = 35;
513  coup2W = settingsPtr->parm("HiggsH2:coup2W");
514  }
515  else if (higgsType == 3) {
516  nameSave = "f fbar -> A0(A3) W+-";
517  codeSave = 1045;
518  idRes = 36;
519  coup2W = settingsPtr->parm("HiggsA3:coup2W");
520  }
521 
522  // Store W+- mass and width for propagator. Common coupling factor.
523  mW = particleDataPtr->m0(24);
524  widW = particleDataPtr->mWidth(24);
525  mWS = mW*mW;
526  mwWS = pow2(mW * widW);
527  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW());
528 
529  // Secondary open width fractions.
530  openFracPairPos = particleDataPtr->resOpenFrac(idRes, 24);
531  openFracPairNeg = particleDataPtr->resOpenFrac(idRes, -24);
532 
533 }
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  sigma0 = (M_PI / sH2) * 2. * pow2(alpEM * thetaWRat * coup2W)
543  * (tH * uH - s3 * s4 + 2. * sH * s4) / (pow2(sH - mWS) + mwWS);
544 
545 }
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  double sigma = sigma0;
555  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
556 
557  // Secondary width for H0 and W+-.
558  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
559  sigma *= (idUp > 0) ? openFracPairPos : openFracPairNeg;
560 
561  // Answer.
562  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  int sign = 1 - 2 * (abs(id1)%2);
574  if (id1 < 0) sign = -sign;
575  setId( id1, id2, idRes, 24 * sign);
576 
577  // Colour flow topologies. Swap when antiquarks.
578  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
579  else setColAcol( 0, 0, 0, 0, 0, 0);
580  if (id1 < 0) swapColAcol();
581 
582 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
593 
594  // For Higgs decay hand over to standard routine.
595  if (idMother == 25 || idMother == 35 || idMother == 36)
596  return weightHiggsDecay( process, iResBeg, iResEnd);
597 
598  // For top decay hand over to standard routine.
599  if (idMother == 6)
600  return weightTopDecay( process, iResBeg, iResEnd);
601 
602  // If not decay of W+- created along with Higgs then done.
603  if (iResBeg != 5 || iResEnd != 6) return 1.;
604 
605  // Order so that fbar(1) f(2) -> H() f'(3) fbar'(4).
606  int i1 = (process[3].id() < 0) ? 3 : 4;
607  int i2 = 7 - i1;
608  int i3 = process[6].daughter1();
609  int i4 = process[6].daughter2();
610  if (process[i3].id() < 0) swap( i3, i4);
611 
612  // Evaluate relevant four-products.
613  double pp13 = process[i1].p() * process[i3].p();
614  double pp14 = process[i1].p() * process[i4].p();
615  double pp23 = process[i2].p() * process[i3].p();
616  double pp24 = process[i2].p() * process[i4].p();
617 
618  // Weight and maximum.
619  double wt = pp13 * pp24;
620  double wtMax = (pp13 + pp14) * (pp23 + pp24);
621 
622  // Done.
623  return wt / wtMax;
624 
625 }
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  if (higgsType == 0) {
641  nameSave = "f f' -> H0 f f'(Z0 Z0 fusion) (SM)";
642  codeSave = 906;
643  idRes = 25;
644  coup2Z = 1.;
645  }
646  else if (higgsType == 1) {
647  nameSave = "f f' -> h0(H1) f f' (Z0 Z0 fusion)";
648  codeSave = 1006;
649  idRes = 25;
650  coup2Z = settingsPtr->parm("HiggsH1:coup2Z");
651  }
652  else if (higgsType == 2) {
653  nameSave = "f f' -> H0(H2) f f' (Z0 Z0 fusion)";
654  codeSave = 1026;
655  idRes = 35;
656  coup2Z = settingsPtr->parm("HiggsH2:coup2Z");
657  }
658  else if (higgsType == 3) {
659  nameSave = "f f' -> A0(A3) f f' (Z0 Z0 fusion)";
660  codeSave = 1046;
661  idRes = 36;
662  coup2Z = settingsPtr->parm("HiggsA3:coup2Z");
663  }
664 
665  // Common fixed mass and coupling factor.
666  mZS = pow2( particleDataPtr->m0(23) );
667  prefac = 0.25 * mZS * pow3( 4. * M_PI / (couplingsPtr->sin2thetaW()
668  * couplingsPtr->cos2thetaW()) );
669 
670  // Secondary open width fraction.
671  openFrac = particleDataPtr->resOpenFrac(idRes);
672 
673 }
674 
675 //--------------------------------------------------------------------------
676 
677 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
678 
679 void Sigma3ff2HfftZZ::sigmaKin() {
680 
681  // Required four-vector products.
682  double pp12 = 0.5 * sH;
683  double pp14 = 0.5 * mH * p4cm.pNeg();
684  double pp15 = 0.5 * mH * p5cm.pNeg();
685  double pp24 = 0.5 * mH * p4cm.pPos();
686  double pp25 = 0.5 * mH * p5cm.pPos();
687  double pp45 = p4cm * p5cm;
688 
689  // Propagator factors and two possible numerators.
690  double prop = pow2( (2. * pp14 + mZS) * (2. * pp25 + mZS) );
691  sigma1 = prefac * pp12 * pp45 / prop;
692  sigma2 = prefac * pp15 * pp24 / prop;
693 
694 }
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  int id1Abs = abs(id1);
704  int id2Abs = abs(id2);
705  double lf1S = pow2( couplingsPtr->lf(id1Abs) );
706  double rf1S = pow2( couplingsPtr->rf(id1Abs) );
707  double lf2S = pow2( couplingsPtr->lf(id2Abs) );
708  double rf2S = pow2( couplingsPtr->rf(id2Abs) );
709  double c1 = lf1S * lf2S + rf1S * rf2S;
710  double c2 = lf1S * rf2S + rf1S * lf2S;
711 
712  // Combine couplings and kinematics factors.
713  double sigma = pow3(alpEM) * (c1 * sigma1 + c2 * sigma2) * pow2(coup2Z);
714 
715  // Secondary width for H0, H1, H2 or A3.
716  sigma *= openFrac;
717 
718  // Answer.
719  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  setId( id1, id2, idRes, id1, id2);
731 
732  // Colour flow topologies. Swap when antiquarks.
733  if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
734  setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
735  else if (abs(id1) < 9 && abs(id2) < 9)
736  setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
737  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
738  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0);
739  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
740  if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
741  swapColAcol();
742 
743 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
754 
755  // For Higgs decay hand over to standard routine.
756  if (idMother == 25 || idMother == 35 || idMother == 36)
757  return weightHiggsDecay( process, iResBeg, iResEnd);
758 
759  // For top decay hand over to standard routine.
760  if (idMother == 6)
761  return weightTopDecay( process, iResBeg, iResEnd);
762 
763  // Else done.
764  return 1.;
765 
766 }
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  if (higgsType == 0) {
781  nameSave = "f_1 f_2 -> H0 f_3 f_4 (W+ W- fusion) (SM)";
782  codeSave = 907;
783  idRes = 25;
784  coup2W = 1.;
785  }
786  else if (higgsType == 1) {
787  nameSave = "f_1 f_2 -> h0(H1) f_3 f_4 (W+ W- fusion)";
788  codeSave = 1007;
789  idRes = 25;
790  coup2W = settingsPtr->parm("HiggsH1:coup2W");
791  }
792  else if (higgsType == 2) {
793  nameSave = "f_1 f_2 -> H0(H2) f_3 f_4 (W+ W- fusion)";
794  codeSave = 1027;
795  idRes = 35;
796  coup2W = settingsPtr->parm("HiggsH2:coup2W");
797  }
798  else if (higgsType == 3) {
799  nameSave = "f_1 f_2 -> A0(A3) f_3 f_4 (W+ W- fusion)";
800  codeSave = 1047;
801  idRes = 36;
802  coup2W = settingsPtr->parm("HiggsA3:coup2W");
803  }
804 
805  // Common fixed mass and coupling factor.
806  mWS = pow2( particleDataPtr->m0(24) );
807  prefac = mWS * pow3( 4. * M_PI / couplingsPtr->sin2thetaW() );
808 
809  // Secondary open width fraction.
810  openFrac = particleDataPtr->resOpenFrac(idRes);
811 
812 }
813 
814 //--------------------------------------------------------------------------
815 
816 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
817 
818 void Sigma3ff2HfftWW::sigmaKin() {
819 
820  // Required four-vector products.
821  double pp12 = 0.5 * sH;
822  double pp14 = 0.5 * mH * p4cm.pNeg();
823  double pp25 = 0.5 * mH * p5cm.pPos();
824  double pp45 = p4cm * p5cm;
825 
826  // Cross section: kinematics part. Combine with couplings.
827  double prop = pow2( (2. * pp14 + mWS) * (2. * pp25 + mWS) );
828  sigma0 = prefac * pp12 * pp45 * pow2(coup2W) / prop;
829 
830 }
831 
832 //--------------------------------------------------------------------------
833 
834 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
835 
836 double Sigma3ff2HfftWW::sigmaHat() {
837 
838  // Some flavour combinations not possible.
839  int id1Abs = abs(id1);
840  int id2Abs = abs(id2);
841  if ( (id1Abs%2 == id2Abs%2 && id1 * id2 > 0)
842  || (id1Abs%2 != id2Abs%2 && id1 * id2 < 0) ) return 0.;
843 
844  // Basic cross section. CKM factors for final states.
845  double sigma = sigma0 * pow3(alpEM) * couplingsPtr->V2CKMsum(id1Abs)
846  * couplingsPtr->V2CKMsum(id2Abs);
847 
848  // Secondary width for H0, H1, H2 or A3.
849  sigma *= openFrac;
850 
851  // Spin-state extra factor 2 per incoming neutrino.
852  if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
853  if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
854 
855  // Answer.
856  return sigma;
857 
858 }
859 
860 //--------------------------------------------------------------------------
861 
862 // Select identity, colour and anticolour.
863 
864 void Sigma3ff2HfftWW::setIdColAcol() {
865 
866  // Pick out-flavours by relative CKM weights.
867  id4 = couplingsPtr->V2CKMpick(id1);
868  id5 = couplingsPtr->V2CKMpick(id2);
869  setId( id1, id2, idRes, id4, id5);
870 
871  // Colour flow topologies. Swap when antiquarks.
872  if (abs(id1) < 9 && abs(id2) < 9 && id1*id2 > 0)
873  setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
874  else if (abs(id1) < 9 && abs(id2) < 9)
875  setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
876  else if (abs(id1) < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
877  else if (abs(id2) < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0);
878  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
879  if ( (abs(id1) < 9 && id1 < 0) || (abs(id1) > 10 && id2 < 0) )
880  swapColAcol();
881 
882 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
893 
894  // For Higgs decay hand over to standard routine.
895  if (idMother == 25 || idMother == 35 || idMother == 36)
896  return weightHiggsDecay( process, iResBeg, iResEnd);
897 
898  // For top decay hand over to standard routine.
899  if (idMother == 6)
900  return weightTopDecay( process, iResBeg, iResEnd);
901 
902  // Else done.
903  return 1.;
904 
905 }
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  if (higgsType == 0 && idNew == 6) {
921  nameSave = "g g -> H t tbar (SM)";
922  codeSave = 908;
923  idRes = 25;
924  coup2Q = 1.;
925  }
926  else if (higgsType == 1 && idNew == 6) {
927  nameSave = "g g -> h0(H1) t tbar";
928  codeSave = 1008;
929  idRes = 25;
930  coup2Q = settingsPtr->parm("HiggsH1:coup2u");
931  }
932  else if (higgsType == 2 && idNew == 6) {
933  nameSave = "g g -> H0(H2) t tbar";
934  codeSave = 1028;
935  idRes = 35;
936  coup2Q = settingsPtr->parm("HiggsH2:coup2u");
937  }
938  else if (higgsType == 3 && idNew == 6) {
939  nameSave = "g g -> A0(A3) t tbar";
940  codeSave = 1048;
941  idRes = 36;
942  coup2Q = settingsPtr->parm("HiggsA3:coup2u");
943  }
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  if (higgsType == 0 && idNew == 5) {
948  nameSave = "g g -> H b bbar (SM)";
949  codeSave = 912;
950  idRes = 25;
951  coup2Q = 1.;
952  }
953  else if (higgsType == 1 && idNew == 5) {
954  nameSave = "g g -> h0(H1) b bbar";
955  codeSave = 1012;
956  idRes = 25;
957  coup2Q = settingsPtr->parm("HiggsH1:coup2d");
958  }
959  else if (higgsType == 2 && idNew == 5) {
960  nameSave = "g g -> H0(H2) b bbar";
961  codeSave = 1032;
962  idRes = 35;
963  coup2Q = settingsPtr->parm("HiggsH2:coup2d");
964  }
965  else if (higgsType == 3 && idNew == 5) {
966  nameSave = "g g -> A0(A3) b bbar";
967  codeSave = 1052;
968  idRes = 36;
969  coup2Q = settingsPtr->parm("HiggsA3:coup2d");
970  }
971 
972  // Common mass and coupling factors.
973  double mWS = pow2(particleDataPtr->m0(24));
974  prefac = (4. * M_PI / couplingsPtr->sin2thetaW()) * pow2(4. * M_PI)
975  * 0.25 / mWS;
976 
977  // Secondary open width fraction.
978  openFracTriplet = particleDataPtr->resOpenFrac(idRes, idNew, -idNew);
979 
980 }
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  double mQ2run = pow2( particleDataPtr->mRun(idNew, mH) );
990 
991  // Linear combination of p_Q and p_Qbar to ensure common mass.
992  double mQ2 = m4 * m5;
993  double epsi = 0.;
994  if (m4 != m5) {
995  double s45 = (p4cm + p5cm).m2Calc();
996  mQ2 = 0.5 * (s4 + s5) - 0.25 * pow2(s4 - s5) / s45;
997  epsi = 0.5 * (s5 - s4) / s45;
998  }
999 
1000  // Set up kinematics: g(4) g(5) -> H(3) Q(1) Qbar(2) in outgoing sense.
1001  Vec4 pTemp[6];
1002  pTemp[4] = Vec4( 0., 0., -0.5* mH, -0.5* mH);
1003  pTemp[5] = Vec4( 0., 0., 0.5* mH, -0.5* mH);
1004  pTemp[1] = p4cm + epsi * (p4cm + p5cm);
1005  pTemp[2] = p5cm - epsi * (p4cm + p5cm);
1006  pTemp[3] = p3cm;
1007 
1008  // Four-product combinations.
1009  double z1 = pTemp[1] * pTemp[2];
1010  double z2 = pTemp[1] * pTemp[3];
1011  double z3 = pTemp[1] * pTemp[4];
1012  double z4 = pTemp[1] * pTemp[5];
1013  double z5 = pTemp[2] * pTemp[3];
1014  double z6 = pTemp[2] * pTemp[4];
1015  double z7 = pTemp[2] * pTemp[5];
1016  double z8 = pTemp[3] * pTemp[4];
1017  double z9 = pTemp[3] * pTemp[5];
1018  double z10 = pTemp[4] * pTemp[5];
1019 
1020  // Powers required as shorthand in matriz elements.
1021  double mQ4 = mQ2 * mQ2;
1022  double mQ6 = mQ2 * mQ4;
1023  double z1S = z1 * z1;
1024  double z2S = z2 * z2;
1025  double z3S = z3 * z3;
1026  double z4S = z4 * z4;
1027  double z5S = z5 * z5;
1028  double z6S = z6 * z6;
1029  double z7S = z7 * z7;
1030  double z8S = z8 * z8;
1031  double z9S = z9 * z9;
1032  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  double fm[9][9];
1037  fm[1][1] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+2*z2+z4+z9+2*
1038  z7+z5)+8*mQ2*s3*(-z1-z4+2*z7)+16*mQ2*(z2*z9+4*z2*
1039  z7+z2*z5-2*z4*z7-2*z9*z7)+8*s3*z4*z7-16*z2*z9*z7;
1040  fm[1][2] = 16*mQ6+8*mQ4*(-2*z1+z2-2*z3-2*z4-4*z10+z9-z8+2
1041  *z7-4*z6+z5)+8*mQ2*(-2*z1*z2-2*z2*z4-2*z2*z10+z2*z7-2*
1042  z2*z6-2*z3*z7+2*z4*z7+4*z10*z7-z9*z7-z8*z7)+16*z2*z7*(z4+
1043  z10);
1044  fm[1][3] = 16*mQ6-4*mQ4*s3+8*mQ4*(-2*z1+2*z2-2*z3-4*
1045  z4-8*z10+z9+z8-2*z7-4*z6+2*z5)-(4*mQ2*s3)*(z1+z4+z10
1046  +z6)+8*mQ2*(-2*z1*z2-2*z1*z10+z1*z9+z1*z8-2*z1*z5+z2S
1047  -4*z2*z4-5*z2*z10+z2*z8-z2*z7-3*z2*z6+z2*z5+z3*z9+2*z3*z7
1048  -z3*z5+z4*z8+2*z4*z6-3*z4*z5-5*z10*z5+z9*z8+z9*z6+z9*z5+
1049  z8*z7-4*z6*z5+z5S)-(16*z2*z5)*(z1+z4+z10+z6);
1050  fm[1][4] = 16*mQ6+4*mQ4*s3+16*mQ4*(-z1+z2-z3-z4+z10-
1051  z9-z8+2*z7+2*z6-z5)+4*mQ2*s3*(z1+z3+z4+z10+2*z7+2*z6
1052  )+8*mQ2*(4*z1*z10+4*z1*z7+4*z1*z6+2*z2*z10-z2*z9-z2*z8+
1053  4*z2*z7+4*z2*z6-z2*z5+4*z10*z5+4*z7*z5+4*z6*z5)-(8*s3*
1054  z1)*(z10+z7+z6)+16*z2*z5*(z10+z7+z6);
1055  fm[1][5] = 8*mQ4*(-2*z1-2*z4+z10-z9)+4*mQ2*(4*z1S-2*z1*
1056  z2+8*z1*z3+6*z1*z10-2*z1*z9+4*z1*z8+4*z1*z7+4*z1*z6+2*z1*
1057  z5+z2*z10+4*z3*z4-z3*z9+2*z3*z7+3*z4*z8-2*z4*z6+2*z4*z5-4
1058  *z10*z7+3*z10*z5-3*z9*z6+3*z8*z7-4*z7S+4*z7*z5)+8*(z1S
1059  *z9-z1S*z8-z1*z2*z7+z1*z2*z6+z1*z3*z9+z1*z3*z5-z1*z4*
1060  z8-z1*z4*z5+z1*z10*z9+z1*z9*z7+z1*z9*z6-z1*z8*z7-z2*z3*z7
1061  +z2*z4*z6-z2*z10*z7-z2*z7S+z3*z7*z5-z4*z10*z5-z4*z7*z5-
1062  z4*z6*z5);
1063  fm[1][6] = 16*mQ4*(-4*z1-z4+z9-z7)+4*mQ2*s3*(-2*z1-z4-
1064  z7)+16*mQ2*(-2*z1S-3*z1*z2-2*z1*z4-3*z1*z9-2*z1*z7-3*
1065  z1*z5-2*z2*z4-2*z7*z5)-8*s3*z4*z7+8*(-z1*z2*z9-2*z1*z2
1066  *z5-z1*z9S-z1*z9*z5+z2S*z7-z2*z4*z5+z2*z9*z7-z2*z7*z5
1067  +z4*z9*z5+z4*z5S);
1068  fm[1][7] = 8*mQ4*(2*z3+z4+3*z10+z9+2*z8+3*z7+6*z6)+2*mQ2*
1069  s3*(-2*z3-z4+3*z10+3*z7+6*z6)+4*mQ2*(4*z1*z10+4*z1*
1070  z7+8*z1*z6+6*z2*z10+z2*z9+2*z2*z8+6*z2*z7+12*z2*z6-8*z3*
1071  z7+4*z4*z7+4*z4*z6+4*z10*z5+4*z9*z7+4*z9*z6-8*z8*z7+4*z7*
1072  z5+8*z6*z5)+4*s3*(-z1*z10-z1*z7-2*z1*z6+2*z3*z7-z4*z7-
1073  z4*z6)+8*z2*(z10*z5+z9*z7+z9*z6-2*z8*z7+z7*z5+2*z6*z5);
1074  fm[1][8] = 8*mQ4*(2*z3+z4+3*z10+2*z9+z8+3*z7+6*z6)+2*mQ2*
1075  s3*(-2*z3-z4+2*z10+z7+2*z6)+4*mQ2*(4*z1*z10-2*z1*z9+
1076  2*z1*z8+4*z1*z7+8*z1*z6+5*z2*z10+2*z2*z9+z2*z8+4*z2*z7+8*
1077  z2*z6-z3*z9-8*z3*z7+2*z3*z5+2*z4*z9-z4*z8+4*z4*z7+4*z4*z6
1078  +4*z4*z5+5*z10*z5+z9S-z9*z8+2*z9*z7+5*z9*z6+z9*z5-7*z8*
1079  z7+2*z8*z5+2*z7*z5+10*z6*z5)+2*s3*(-z1*z10+z3*z7-2*z4*
1080  z7+z4*z6)+4*(-z1*z9S+z1*z9*z8-2*z1*z9*z5-z1*z8*z5+2*z2*
1081  z10*z5+z2*z9*z7+z2*z9*z6-2*z2*z8*z7+3*z2*z6*z5+z3*z9*z5+
1082  z3*z5S+z4*z9*z5-2*z4*z8*z5+2*z4*z5S);
1083  fm[2][2] = 16*mQ6+16*mQ4*(-z1+z3-z4-z10+z7-z6)+16*mQ2*(
1084  z3*z10+z3*z7+z3*z6+z4*z7+z10*z7)-16*z3*z10*z7;
1085  fm[2][3] = 16*mQ6+8*mQ4*(-2*z1+z2+2*z3-4*z4-4*z10-z9+z8-2
1086  *z7-2*z6+z5)+8*mQ2*(-2*z1*z5+4*z3*z10-z3*z9-z3*z8-2*z3*
1087  z7+2*z3*z6+z3*z5-2*z4*z5-2*z10*z5-2*z6*z5)+16*z3*z5*(z10+
1088  z6);
1089  fm[2][4] = 8*mQ4*(-2*z1-2*z3+z10-z8)+4*mQ2*(4*z1S-2*z1*
1090  z2+8*z1*z4+6*z1*z10+4*z1*z9-2*z1*z8+4*z1*z7+4*z1*z6+2*z1*
1091  z5+z2*z10+4*z3*z4+3*z3*z9-2*z3*z7+2*z3*z5-z4*z8+2*z4*z6-4
1092  *z10*z6+3*z10*z5+3*z9*z6-3*z8*z7-4*z6S+4*z6*z5)+8*(-z1S
1093  *z9+z1S*z8+z1*z2*z7-z1*z2*z6-z1*z3*z9-z1*z3*z5+z1*z4
1094  *z8+z1*z4*z5+z1*z10*z8-z1*z9*z6+z1*z8*z7+z1*z8*z6+z2*z3*
1095  z7-z2*z4*z6-z2*z10*z6-z2*z6S-z3*z10*z5-z3*z7*z5-z3*z6*
1096  z5+z4*z6*z5);
1097  fm[2][5] = 16*mQ4*z10+8*mQ2*(2*z1S+2*z1*z3+2*z1*z4+2*z1
1098  *z10+2*z1*z7+2*z1*z6+z3*z7+z4*z6)+8*(-2*pow3(z1)-2*z1S*z3-
1099  2*z1S*z4-2*z1S*z10-2*z1S*z7-2*z1S*z6-2*z1*z3*z4-
1100  z1*z3*z10-2*z1*z3*z6-z1*z4*z10-2*z1*z4*z7-z1*z10S-z1*
1101  z10*z7-z1*z10*z6-2*z1*z7*z6+z3S*z7-z3*z4*z7-z3*z4*z6+z3
1102  *z10*z7+z3*z7S-z3*z7*z6+z4S*z6+z4*z10*z6-z4*z7*z6+z4*
1103  z6S);
1104  fm[2][6] = 8*mQ4*(-2*z1+z10-z9-2*z7)+4*mQ2*(4*z1S+2*z1*
1105  z2+4*z1*z3+4*z1*z4+6*z1*z10-2*z1*z9+4*z1*z8+8*z1*z6-2*z1*
1106  z5+4*z2*z4+3*z2*z10+2*z2*z7-3*z3*z9-2*z3*z7-4*z4S-4*z4*
1107  z10+3*z4*z8+2*z4*z6+z10*z5-z9*z6+3*z8*z7+4*z7*z6)+8*(z1S
1108  *z9-z1S*z8-z1*z2*z7+z1*z2*z6+z1*z3*z9+z1*z3*z5+z1*z4*
1109  z9-z1*z4*z8-z1*z4*z5+z1*z10*z9+z1*z9*z6-z1*z8*z7-z2*z3*z7
1110  -z2*z4*z7+z2*z4*z6-z2*z10*z7+z3*z7*z5-z4S*z5-z4*z10*z5-
1111  z4*z6*z5);
1112  fm[2][7] = 8*mQ4*(z3+2*z4+3*z10+z7+2*z6)+4*mQ2*(-4*z1*z3-
1113  2*z1*z4-2*z1*z10+z1*z9-z1*z8-4*z1*z7-2*z1*z6+z2*z3+2*z2*
1114  z4+3*z2*z10+z2*z7+2*z2*z6-6*z3*z4-6*z3*z10-2*z3*z9-2*z3*
1115  z7-4*z3*z6-z3*z5-6*z4S-6*z4*z10-3*z4*z9-z4*z8-4*z4*z7-2
1116  *z4*z6-2*z4*z5-3*z10*z9-3*z10*z8-6*z10*z7-6*z10*z6+z10*z5
1117  +z9*z7-2*z8*z7-2*z8*z6-6*z7*z6+z7*z5-6*z6S+2*z6*z5)+4*(
1118  -z1S*z9+z1S*z8-2*z1*z2*z10-3*z1*z2*z7-3*z1*z2*z6+z1*
1119  z3*z9-z1*z3*z5+z1*z4*z9+z1*z4*z8+z1*z4*z5+z1*z10*z9+z1*
1120  z10*z8-z1*z9*z6+z1*z8*z6+z2*z3*z7-3*z2*z4*z7-z2*z4*z6-3*
1121  z2*z10*z7-3*z2*z10*z6-3*z2*z7*z6-3*z2*z6S-2*z3*z4*z5-z3
1122  *z10*z5-z3*z6*z5-z4S*z5-z4*z10*z5+z4*z6*z5);
1123  fm[2][8] = 8*mQ4*(z3+2*z4+3*z10+z7+2*z6)+4*mQ2*(-4*z1*z3-
1124  2*z1*z4-2*z1*z10-z1*z9+z1*z8-4*z1*z7-2*z1*z6+z2*z3+2*z2*
1125  z4+z2*z10-z2*z7-2*z2*z6-6*z3*z4-6*z3*z10-2*z3*z9+z3*z8-2*
1126  z3*z7-4*z3*z6+z3*z5-6*z4S-6*z4*z10-2*z4*z9-4*z4*z7-2*z4
1127  *z6+2*z4*z5-3*z10*z9-3*z10*z8-6*z10*z7-6*z10*z6+3*z10*z5-
1128  z9*z6-2*z8*z7-3*z8*z6-6*z7*z6+z7*z5-6*z6S+2*z6*z5)+4*(
1129  z1S*z9-z1S*z8-z1*z2*z7+z1*z2*z6-3*z1*z3*z5+z1*z4*z9-
1130  z1*z4*z8-3*z1*z4*z5+z1*z10*z9+z1*z10*z8-2*z1*z10*z5+z1*z9
1131  *z6+z1*z8*z7+z1*z8*z6-z2*z4*z7+z2*z4*z6-z2*z10*z7-z2*z10*
1132  z6-2*z2*z7*z6-z2*z6S-3*z3*z4*z5-3*z3*z10*z5+z3*z7*z5-3*
1133  z3*z6*z5-3*z4S*z5-3*z4*z10*z5-z4*z6*z5);
1134  fm[3][3] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+z2+2*z3+z8+z6
1135  +2*z5)+8*mQ2*s3*(-z1+2*z3-z6)+16*mQ2*(z2*z5-2*z3*
1136  z8-2*z3*z6+4*z3*z5+z8*z5)+8*s3*z3*z6-16*z3*z8*z5;
1137  fm[3][4] = 16*mQ4*(-4*z1-z3+z8-z6)+4*mQ2*s3*(-2*z1-z3-
1138  z6)+16*mQ2*(-2*z1S-3*z1*z2-2*z1*z3-3*z1*z8-2*z1*z6-3*
1139  z1*z5-2*z2*z3-2*z6*z5)-8*s3*z3*z6+8*(-z1*z2*z8-2*z1*z2
1140  *z5-z1*z8S-z1*z8*z5+z2S*z6-z2*z3*z5+z2*z8*z6-z2*z6*z5
1141  +z3*z8*z5+z3*z5S);
1142  fm[3][5] = 8*mQ4*(-2*z1+z10-z8-2*z6)+4*mQ2*(4*z1S+2*z1*
1143  z2+4*z1*z3+4*z1*z4+6*z1*z10+4*z1*z9-2*z1*z8+8*z1*z7-2*z1*
1144  z5+4*z2*z3+3*z2*z10+2*z2*z6-4*z3S-4*z3*z10+3*z3*z9+2*z3
1145  *z7-3*z4*z8-2*z4*z6+z10*z5+3*z9*z6-z8*z7+4*z7*z6)+8*(-z1S
1146  *z9+z1S*z8+z1*z2*z7-z1*z2*z6-z1*z3*z9+z1*z3*z8-z1*z3
1147  *z5+z1*z4*z8+z1*z4*z5+z1*z10*z8-z1*z9*z6+z1*z8*z7+z2*z3*
1148  z7-z2*z3*z6-z2*z4*z6-z2*z10*z6-z3S*z5-z3*z10*z5-z3*z7*
1149  z5+z4*z6*z5);
1150  fm[3][6] = 16*mQ6+4*mQ4*s3+16*mQ4*(-z1-z2+2*z3+2*z4+
1151  z10-z9-z8-z7-z6+z5)+4*mQ2*s3*(z1+2*z3+2*z4+z10+z7+z6
1152  )+8*mQ2*(4*z1*z3+4*z1*z4+4*z1*z10+4*z2*z3+4*z2*z4+4*z2*
1153  z10-z2*z5+4*z3*z5+4*z4*z5+2*z10*z5-z9*z5-z8*z5)-(8*s3*
1154  z1)*(z3+z4+z10)+16*z2*z5*(z3+z4+z10);
1155  fm[3][7] = 8*mQ4*(3*z3+6*z4+3*z10+z9+2*z8+2*z7+z6)+2*mQ2*
1156  s3*(z3+2*z4+2*z10-2*z7-z6)+4*mQ2*(4*z1*z3+8*z1*z4+4*
1157  z1*z10+2*z1*z9-2*z1*z8+2*z2*z3+10*z2*z4+5*z2*z10+2*z2*z9+
1158  z2*z8+2*z2*z7+4*z2*z6-7*z3*z9+2*z3*z8-8*z3*z7+4*z3*z6+4*
1159  z3*z5+5*z4*z8+4*z4*z6+8*z4*z5+5*z10*z5-z9*z8-z9*z6+z9*z5+
1160  z8S-z8*z7+2*z8*z6+2*z8*z5)+2*s3*(-z1*z10+z3*z7-2*z3*
1161  z6+z4*z6)+4*(-z1*z2*z9-2*z1*z2*z8+z1*z9*z8-z1*z8S+z2S
1162  *z7+2*z2S*z6+3*z2*z4*z5+2*z2*z10*z5-2*z2*z9*z6+z2*z8*z7
1163  +z2*z8*z6-2*z3*z9*z5+z3*z8*z5+z4*z8*z5);
1164  fm[3][8] = 8*mQ4*(3*z3+6*z4+3*z10+2*z9+z8+2*z7+z6)+2*mQ2*
1165  s3*(3*z3+6*z4+3*z10-2*z7-z6)+4*mQ2*(4*z1*z3+8*z1*z4+
1166  4*z1*z10+4*z2*z3+8*z2*z4+4*z2*z10-8*z3*z9+4*z3*z8-8*z3*z7
1167  +4*z3*z6+6*z3*z5+4*z4*z8+4*z4*z6+12*z4*z5+6*z10*z5+2*z9*
1168  z5+z8*z5)+4*s3*(-z1*z3-2*z1*z4-z1*z10+2*z3*z7-z3*z6-z4
1169  *z6)+8*z5*(z2*z3+2*z2*z4+z2*z10-2*z3*z9+z3*z8+z4*z8);
1170  fm[4][4] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+2*z2+z3+z8+2*
1171  z6+z5)+8*mQ2*s3*(-z1-z3+2*z6)+16*mQ2*(z2*z8+4*z2*
1172  z6+z2*z5-2*z3*z6-2*z8*z6)+8*s3*z3*z6-16*z2*z8*z6;
1173  fm[4][5] = 16*mQ6+8*mQ4*(-2*z1+z2-2*z3-2*z4-4*z10-z9+z8-4
1174  *z7+2*z6+z5)+8*mQ2*(-2*z1*z2-2*z2*z3-2*z2*z10-2*z2*z7+
1175  z2*z6+2*z3*z6-2*z4*z6+4*z10*z6-z9*z6-z8*z6)+16*z2*z6*(z3+
1176  z10);
1177  fm[4][6] = 16*mQ6-4*mQ4*s3+8*mQ4*(-2*z1+2*z2-4*z3-2*
1178  z4-8*z10+z9+z8-4*z7-2*z6+2*z5)-(4*mQ2*s3)*(z1+z3+z10
1179  +z7)+8*mQ2*(-2*z1*z2-2*z1*z10+z1*z9+z1*z8-2*z1*z5+z2S
1180  -4*z2*z3-5*z2*z10+z2*z9-3*z2*z7-z2*z6+z2*z5+z3*z9+2*z3*z7
1181  -3*z3*z5+z4*z8+2*z4*z6-z4*z5-5*z10*z5+z9*z8+z9*z6+z8*z7+
1182  z8*z5-4*z7*z5+z5S)-(16*z2*z5)*(z1+z3+z10+z7);
1183  fm[4][7] = 8*mQ4*(-z3-2*z4-3*z10-2*z9-z8-6*z7-3*z6)+2*mQ2
1184  *s3*(z3+2*z4-3*z10-6*z7-3*z6)+4*mQ2*(-4*z1*z10-8*z1*
1185  z7-4*z1*z6-6*z2*z10-2*z2*z9-z2*z8-12*z2*z7-6*z2*z6-4*z3*
1186  z7-4*z3*z6+8*z4*z6-4*z10*z5+8*z9*z6-4*z8*z7-4*z8*z6-8*z7*
1187  z5-4*z6*z5)+4*s3*(z1*z10+2*z1*z7+z1*z6+z3*z7+z3*z6-2*
1188  z4*z6)+8*z2*(-z10*z5+2*z9*z6-z8*z7-z8*z6-2*z7*z5-z6*z5);
1189  fm[4][8] = 8*mQ4*(-z3-2*z4-3*z10-z9-2*z8-6*z7-3*z6)+2*mQ2
1190  *s3*(z3+2*z4-2*z10-2*z7-z6)+4*mQ2*(-4*z1*z10-2*z1*z9
1191  +2*z1*z8-8*z1*z7-4*z1*z6-5*z2*z10-z2*z9-2*z2*z8-8*z2*z7-4
1192  *z2*z6+z3*z9-2*z3*z8-4*z3*z7-4*z3*z6-4*z3*z5+z4*z8+8*z4*
1193  z6-2*z4*z5-5*z10*z5+z9*z8+7*z9*z6-2*z9*z5-z8S-5*z8*z7-2
1194  *z8*z6-z8*z5-10*z7*z5-2*z6*z5)+2*s3*(z1*z10-z3*z7+2*z3
1195  *z6-z4*z6)+4*(-z1*z9*z8+z1*z9*z5+z1*z8S+2*z1*z8*z5-2*z2
1196  *z10*z5+2*z2*z9*z6-z2*z8*z7-z2*z8*z6-3*z2*z7*z5+2*z3*z9*
1197  z5-z3*z8*z5-2*z3*z5S-z4*z8*z5-z4*z5S);
1198  fm[5][5] = 16*mQ6+16*mQ4*(-z1-z3+z4-z10-z7+z6)+16*mQ2*(
1199  z3*z6+z4*z10+z4*z7+z4*z6+z10*z6)-16*z4*z10*z6;
1200  fm[5][6] = 16*mQ6+8*mQ4*(-2*z1+z2-4*z3+2*z4-4*z10+z9-z8-2
1201  *z7-2*z6+z5)+8*mQ2*(-2*z1*z5-2*z3*z5+4*z4*z10-z4*z9-z4*
1202  z8+2*z4*z7-2*z4*z6+z4*z5-2*z10*z5-2*z7*z5)+16*z4*z5*(z10+
1203  z7);
1204  fm[5][7] = 8*mQ4*(-2*z3-z4-3*z10-2*z7-z6)+4*mQ2*(2*z1*z3+
1205  4*z1*z4+2*z1*z10+z1*z9-z1*z8+2*z1*z7+4*z1*z6-2*z2*z3-z2*
1206  z4-3*z2*z10-2*z2*z7-z2*z6+6*z3S+6*z3*z4+6*z3*z10+z3*z9+
1207  3*z3*z8+2*z3*z7+4*z3*z6+2*z3*z5+6*z4*z10+2*z4*z8+4*z4*z7+
1208  2*z4*z6+z4*z5+3*z10*z9+3*z10*z8+6*z10*z7+6*z10*z6-z10*z5+
1209  2*z9*z7+2*z9*z6-z8*z6+6*z7S+6*z7*z6-2*z7*z5-z6*z5)+4*(-
1210  z1S*z9+z1S*z8+2*z1*z2*z10+3*z1*z2*z7+3*z1*z2*z6-z1*z3
1211  *z9-z1*z3*z8-z1*z3*z5-z1*z4*z8+z1*z4*z5-z1*z10*z9-z1*z10*
1212  z8-z1*z9*z7+z1*z8*z7+z2*z3*z7+3*z2*z3*z6-z2*z4*z6+3*z2*
1213  z10*z7+3*z2*z10*z6+3*z2*z7S+3*z2*z7*z6+z3S*z5+2*z3*z4
1214  *z5+z3*z10*z5-z3*z7*z5+z4*z10*z5+z4*z7*z5);
1215  fm[5][8] = 8*mQ4*(-2*z3-z4-3*z10-2*z7-z6)+4*mQ2*(2*z1*z3+
1216  4*z1*z4+2*z1*z10-z1*z9+z1*z8+2*z1*z7+4*z1*z6-2*z2*z3-z2*
1217  z4-z2*z10+2*z2*z7+z2*z6+6*z3S+6*z3*z4+6*z3*z10+2*z3*z8+
1218  2*z3*z7+4*z3*z6-2*z3*z5+6*z4*z10-z4*z9+2*z4*z8+4*z4*z7+2*
1219  z4*z6-z4*z5+3*z10*z9+3*z10*z8+6*z10*z7+6*z10*z6-3*z10*z5+
1220  3*z9*z7+2*z9*z6+z8*z7+6*z7S+6*z7*z6-2*z7*z5-z6*z5)+4*(
1221  z1S*z9-z1S*z8-z1*z2*z7+z1*z2*z6+z1*z3*z9-z1*z3*z8+3*
1222  z1*z3*z5+3*z1*z4*z5-z1*z10*z9-z1*z10*z8+2*z1*z10*z5-z1*z9
1223  *z7-z1*z9*z6-z1*z8*z7-z2*z3*z7+z2*z3*z6+z2*z10*z7+z2*z10*
1224  z6+z2*z7S+2*z2*z7*z6+3*z3S*z5+3*z3*z4*z5+3*z3*z10*z5+
1225  z3*z7*z5+3*z4*z10*z5+3*z4*z7*z5-z4*z6*z5);
1226  fm[6][6] = 64*mQ6+16*mQ4*s3+32*mQ4*(z1+z2+2*z4+z9+z7
1227  +2*z5)+8*mQ2*s3*(-z1+2*z4-z7)+16*mQ2*(z2*z5-2*z4*
1228  z9-2*z4*z7+4*z4*z5+z9*z5)+8*s3*z4*z7-16*z4*z9*z5;
1229  fm[6][7] = 8*mQ4*(-6*z3-3*z4-3*z10-2*z9-z8-z7-2*z6)+2*mQ2
1230  *s3*(-2*z3-z4-2*z10+z7+2*z6)+4*mQ2*(-8*z1*z3-4*z1*z4
1231  -4*z1*z10+2*z1*z9-2*z1*z8-10*z2*z3-2*z2*z4-5*z2*z10-z2*z9
1232  -2*z2*z8-4*z2*z7-2*z2*z6-5*z3*z9-4*z3*z7-8*z3*z5-2*z4*z9+
1233  7*z4*z8-4*z4*z7+8*z4*z6-4*z4*z5-5*z10*z5-z9S+z9*z8-2*z9
1234  *z7+z9*z6-2*z9*z5+z8*z7-z8*z5)+2*s3*(z1*z10-z3*z7+2*z4
1235  *z7-z4*z6)+4*(2*z1*z2*z9+z1*z2*z8+z1*z9S-z1*z9*z8-2*
1236  z2S*z7-z2S*z6-3*z2*z3*z5-2*z2*z10*z5-z2*z9*z7-z2*z9*z6+
1237  2*z2*z8*z7-z3*z9*z5-z4*z9*z5+2*z4*z8*z5);
1238  fm[6][8] = 8*mQ4*(-6*z3-3*z4-3*z10-z9-2*z8-z7-2*z6)+2*mQ2
1239  *s3*(-6*z3-3*z4-3*z10+z7+2*z6)+4*mQ2*(-8*z1*z3-4*z1*
1240  z4-4*z1*z10-8*z2*z3-4*z2*z4-4*z2*z10-4*z3*z9-4*z3*z7-12*
1241  z3*z5-4*z4*z9+8*z4*z8-4*z4*z7+8*z4*z6-6*z4*z5-6*z10*z5-z9
1242  *z5-2*z8*z5)+4*s3*(2*z1*z3+z1*z4+z1*z10+z3*z7+z4*z7-2*
1243  z4*z6)+8*z5*(-2*z2*z3-z2*z4-z2*z10-z3*z9-z4*z9+2*z4*z8);
1244  fm[7][7] = 72*mQ4*z10+18*mQ2*s3*z10+8*mQ2*(z1*z10+9*
1245  z2*z10+7*z3*z7+2*z3*z6+2*z4*z7+7*z4*z6+z10*z5+2*z9*z7+7*
1246  z9*z6+7*z8*z7+2*z8*z6)+2*s3*(-z1*z10-7*z3*z7-2*z3*z6-2
1247  *z4*z7-7*z4*z6)+4*z2*(z10*z5+2*z9*z7+7*z9*z6+7*z8*z7+2*z8
1248  *z6);
1249  fm[7][8] = 72*mQ4*z10+2*mQ2*s3*z10+4*mQ2*(2*z1*z10+
1250  10*z2*z10+7*z3*z9+2*z3*z8+14*z3*z7+4*z3*z6+2*z4*z9+7*z4*
1251  z8+4*z4*z7+14*z4*z6+10*z10*z5+z9S+7*z9*z8+2*z9*z7+7*z9*
1252  z6+z8S+7*z8*z7+2*z8*z6)+2*s3*(7*z1*z10-7*z3*z7-2*z3*
1253  z6-2*z4*z7-7*z4*z6)+2*(-2*z1*z9S-14*z1*z9*z8-2*z1*z8S
1254  +2*z2*z10*z5+2*z2*z9*z7+7*z2*z9*z6+7*z2*z8*z7+2*z2*z8*z6+
1255  7*z3*z9*z5+2*z3*z8*z5+2*z4*z9*z5+7*z4*z8*z5);
1256  fm[8][8] = 72*mQ4*z10+18*mQ2*s3*z10+8*mQ2*(z1*z10+z2
1257  *z10+7*z3*z9+2*z3*z8+7*z3*z7+2*z3*z6+2*z4*z9+7*z4*z8+2*z4
1258  *z7+7*z4*z6+9*z10*z5)+2*s3*(-z1*z10-7*z3*z7-2*z3*z6-2*
1259  z4*z7-7*z4*z6)+4*z5*(z2*z10+7*z3*z9+2*z3*z8+2*z4*z9+7*z4*
1260  z8);
1261  double fm99 = -4*mQ4*z10-mQ2*s3*z10+4*mQ2*(-z1*z10-z2*z10+
1262  z3*z7+z4*z6-z10*z5+z9*z6+z8*z7)+s3*(z1*z10-z3*z7-z4*z6
1263  )+2*z2*(-z10*z5+z9*z6+z8*z7);
1264  double fm910 = -4*mQ4*z10-mQ2*s3*z10+2*mQ2*(-2*z1*z10-2*z2*
1265  z10+2*z3*z9+2*z3*z7+2*z4*z6-2*z10*z5+z9*z8+2*z8*z7)+s3
1266  *(z1*z10-z3*z7-z4*z6)+2*(-z1*z9*z8-z2*z10*z5+z2*z8*z7+z3*
1267  z9*z5);
1268  double fmxx = -4*mQ4*z10-mQ2*s3*z10+2*mQ2*(-2*z1*z10-2*z2*
1269  z10+2*z4*z8+2*z4*z6+2*z3*z7-2*z10*z5+z9*z8+2*z9*z6)+s3
1270  *(z1*z10-z3*z7-z4*z6)+2*(-z1*z9*z8-z2*z10*z5+z2*z9*z6+z4*
1271  z8*z5);
1272  fm910 = 0.5*(fmxx+fm910);
1273  double fm1010 = -4*mQ4*z10-mQ2*s3*z10+4*mQ2*(-z1*z10-z2*z10+
1274  z3*z7+z4*z6-z10*z5+z9*z3+z8*z4)+s3*(z1*z10-z3*z7-z4*z6
1275  )+2*z5*(-z10*z2+z9*z3+z8*z4);
1276  fm[7][7] -= 2. * fm99;
1277  fm[7][8] -= 2. * fm910;
1278  fm[8][8] -= 2. * fm1010;
1279 
1280  // Propagators.
1281  double ss1 = (pTemp[1] + pTemp[3]).m2Calc() - mQ2;
1282  double ss2 = (pTemp[1] + pTemp[4]).m2Calc() - mQ2;
1283  double ss3 = (pTemp[1] + pTemp[5]).m2Calc() - mQ2;
1284  double ss4 = (pTemp[2] + pTemp[3]).m2Calc() - mQ2;
1285  double ss5 = (pTemp[2] + pTemp[4]).m2Calc() - mQ2;
1286  double ss6 = (pTemp[2] + pTemp[5]).m2Calc() - mQ2;
1287  double ss7 = sH;
1288 
1289  // Propagator combinations.
1290  double dz[9];
1291  dz[1] = ss1 * ss6;
1292  dz[2] = ss2 * ss6;
1293  dz[3] = ss2 * ss4;
1294  dz[4] = ss1 * ss5;
1295  dz[5] = ss3 * ss5;
1296  dz[6] = ss3 * ss4;
1297  dz[7] = ss7 * ss1;
1298  dz[8] = ss7 * ss4;
1299 
1300  // Colour factors.
1301  double clr[9][9];
1302  for (int i = 1; i < 4; ++i)
1303  for (int j = 1; j < 4; ++j) {
1304  clr[i][j] = 16. / 3.;
1305  clr[i][j+3] = -2. / 3.;
1306  clr[i+3][j] = -2. / 3.;
1307  clr[i+3][j+3] = 16. / 3.;
1308  }
1309  for (int i = 1; i < 4; ++i)
1310  for (int j = 1; j < 3; ++j) {
1311  clr[i][j+6] = -6.;
1312  clr[i+3][j+6] = 6.;
1313  clr[j+6][i] = -6.;
1314  clr[j+6][i+3] = 6.;
1315  }
1316  for (int i = 1; i < 3; ++i)
1317  for (int j = 1; j < 3; ++j)
1318  clr[i+6][j+6] = 12.;
1319 
1320  // Produce final result: matrix elements * colours * propagators.
1321  double wtSum = 0.;
1322  for (int i = 1; i < 9; ++i)
1323  for (int j = i; j < 9; ++j) {
1324  double fac = (j == i) ? 4. : 8.;
1325  wtSum += fm[i][j] * fac * clr[i][j] / (dz[i] * dz[j]);
1326  }
1327  wtSum *= -1./256.;
1328 
1329  // Combine factors.
1330  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  sigma *= openFracTriplet;
1335 
1336 }
1337 
1338 //--------------------------------------------------------------------------
1339 
1340 // Select identity, colour and anticolour.
1341 
1342 void Sigma3gg2HQQbar::setIdColAcol() {
1343 
1344  // Pick out-flavours by relative CKM weights.
1345  setId( id1, id2, idRes, idNew, -idNew);
1346 
1347  // Colour flow topologies.
1348  if (rndmPtr->flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 0, 0, 3);
1349  else setColAcol( 1, 2, 3, 1, 0, 0, 3, 0, 0, 2);
1350 
1351 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
1362 
1363  // For Higgs decay hand over to standard routine.
1364  if (idMother == 25 || idMother == 35 || idMother == 36)
1365  return weightHiggsDecay( process, iResBeg, iResEnd);
1366 
1367  // For top decay hand over to standard routine.
1368  if (idMother == 6)
1369  return weightTopDecay( process, iResBeg, iResEnd);
1370 
1371  // Else done.
1372  return 1.;
1373 
1374 }
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  if (higgsType == 0 && idNew == 6) {
1393  nameSave = "q qbar -> H t tbar (SM)";
1394  codeSave = 909;
1395  idRes = 25;
1396  coup2Q = 1.;
1397  }
1398  else if (higgsType == 1 && idNew == 6) {
1399  nameSave = "q qbar -> h0(H1) t tbar";
1400  codeSave = 1009;
1401  idRes = 25;
1402  coup2Q = settingsPtr->parm("HiggsH1:coup2u");
1403  }
1404  else if (higgsType == 2 && idNew == 6) {
1405  nameSave = "q qbar -> H0(H2) t tbar";
1406  codeSave = 1029;
1407  idRes = 35;
1408  coup2Q = settingsPtr->parm("HiggsH2:coup2u");
1409  }
1410  else if (higgsType == 3 && idNew == 6) {
1411  nameSave = "q qbar -> A0(A3) t tbar";
1412  codeSave = 1049;
1413  idRes = 36;
1414  coup2Q = settingsPtr->parm("HiggsA3:coup2u");
1415  }
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  if (higgsType == 0 && idNew == 5) {
1420  nameSave = "q qbar -> H b bbar (SM)";
1421  codeSave = 913;
1422  idRes = 25;
1423  coup2Q = 1.;
1424  }
1425  else if (higgsType == 1 && idNew == 5) {
1426  nameSave = "q qbar -> h0(H1) b bbar";
1427  codeSave = 1013;
1428  idRes = 25;
1429  coup2Q = settingsPtr->parm("HiggsH1:coup2d");
1430  }
1431  else if (higgsType == 2 && idNew == 5) {
1432  nameSave = "q qbar -> H0(H2) b bbar";
1433  codeSave = 1033;
1434  idRes = 35;
1435  coup2Q = settingsPtr->parm("HiggsH2:coup2d");
1436  }
1437  else if (higgsType == 3 && idNew == 5) {
1438  nameSave = "q qbar -> A0(A3) b bbar";
1439  codeSave = 1053;
1440  idRes = 36;
1441  coup2Q = settingsPtr->parm("HiggsA3:coup2d");
1442  }
1443 
1444  // Common mass and coupling factors.
1445  double mWS = pow2(particleDataPtr->m0(24));
1446  prefac = (4. * M_PI / couplingsPtr->sin2thetaW()) * pow2(4. * M_PI)
1447  * 0.25 / mWS;
1448 
1449  // Secondary open width fraction.
1450  openFracTriplet = particleDataPtr->resOpenFrac(idRes, idNew, -idNew);
1451 
1452 }
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  double mQ2run = pow2( particleDataPtr->mRun(idNew, mH) );
1462 
1463  // Linear combination of p_Q and p_Qbar to ensure common mass.
1464  double mQ2 = m4 * m5;
1465  double epsi = 0.;
1466  if (m4 != m5) {
1467  double s45 = (p4cm + p5cm).m2Calc();
1468  mQ2 = 0.5 * (s4 + s5) - 0.25 * pow2(s4 - s5) / s45;
1469  epsi = 0.5 * (s5 - s4) / s45;
1470  }
1471 
1472  // Set up kinematics: q(4) qbar(5) -> H(3) Q(1) Qbar(2) in outgoing sense.
1473  Vec4 pTemp[6];
1474  pTemp[4] = Vec4( 0., 0., -0.5* mH, -0.5* mH);
1475  pTemp[5] = Vec4( 0., 0., 0.5* mH, -0.5* mH);
1476  pTemp[1] = p4cm + epsi * (p4cm + p5cm);
1477  pTemp[2] = p5cm - epsi * (p4cm + p5cm);
1478  pTemp[3] = p3cm;
1479 
1480  // Four-product combinations.
1481  double z1 = pTemp[1] * pTemp[2];
1482  double z2 = pTemp[1] * pTemp[3];
1483  double z3 = pTemp[1] * pTemp[4];
1484  double z4 = pTemp[1] * pTemp[5];
1485  double z5 = pTemp[2] * pTemp[3];
1486  double z6 = pTemp[2] * pTemp[4];
1487  double z7 = pTemp[2] * pTemp[5];
1488  double z8 = pTemp[3] * pTemp[4];
1489  double z9 = pTemp[3] * pTemp[5];
1490  double z10 = pTemp[4] * pTemp[5];
1491 
1492  // Powers required as shorthand in matriz elements.
1493  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  double a11 = -8.*mQ4*z10-2.*mQ2*s3*z10-(8.*mQ2)*(z2*z10+z3
1498  *z7+z4*z6+z9*z6+z8*z7)+2.*s3*(z3*z7+z4*z6)-(4.*z2)*(z9
1499  *z6+z8*z7);
1500  double a12 = -8.*mQ4*z10+4.*mQ2*(-z2*z10-z3*z9-2.*z3*z7-z4*z8-
1501  2.*z4*z6-z10*z5-z9*z8-z9*z6-z8*z7)+2.*s3*(-z1*z10+z3*z7
1502  +z4*z6)+2.*(2.*z1*z9*z8-z2*z9*z6-z2*z8*z7-z3*z9*z5-z4*z8*
1503  z5);
1504  double a22 = -8.*mQ4*z10-2.*mQ2*s3*z10-(8.*mQ2)*(z3*z9+z3*
1505  z7+z4*z8+z4*z6+z10*z5)+2.*s3*(z3*z7+z4*z6)-(4.*z5)*(z3
1506  *z9+z4*z8);
1507 
1508  // Propagators and propagator combinations.
1509  double ss1 = (pTemp[1] + pTemp[3]).m2Calc() - mQ2;
1510  double ss4 = (pTemp[2] + pTemp[3]).m2Calc() - mQ2;
1511  double ss7 = sH;
1512  double dz7 = ss7 * ss1;
1513  double dz8 = ss7 * ss4;
1514 
1515  // Produce final result: matrix elements * propagators.
1516  a11 /= (dz7 * dz7);
1517  a12 /= (dz7 * dz8);
1518  a22 /= (dz8 * dz8);
1519  double wtSum = -(a11 + a22 + 2.*a12) * (8./9.);
1520 
1521  // Combine factors.
1522  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  sigma *= openFracTriplet;
1527 
1528 }
1529 
1530 //--------------------------------------------------------------------------
1531 
1532 // Select identity, colour and anticolour.
1533 
1534 void Sigma3qqbar2HQQbar::setIdColAcol() {
1535 
1536  // Pick out-flavours by relative CKM weights.
1537  setId( id1, id2, idRes, idNew, -idNew);
1538 
1539  // Colour flow topologies.
1540  if (id1 > 0) setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
1541  else setColAcol( 0, 1, 2, 0, 0, 0, 2, 0, 0, 1);
1542 
1543 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
1554 
1555  // For Higgs decay hand over to standard routine.
1556  if (idMother == 25 || idMother == 35 || idMother == 36)
1557  return weightHiggsDecay( process, iResBeg, iResEnd);
1558 
1559  // For top decay hand over to standard routine.
1560  if (idMother == 6)
1561  return weightTopDecay( process, iResBeg, iResEnd);
1562 
1563  // Else done.
1564  return 1.;
1565 
1566 }
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  if (higgsType == 0 && idNew == 4) {
1583  nameSave = "c g -> H c (SM)";
1584  codeSave = 911;
1585  idRes = 25;
1586  }
1587  else if (higgsType == 1 && idNew == 4) {
1588  nameSave = "c g -> h0(H1) c";
1589  codeSave = 1011;
1590  idRes = 25;
1591  }
1592  else if (higgsType == 2 && idNew == 4) {
1593  nameSave = "c g -> H0(H2) c";
1594  codeSave = 1031;
1595  idRes = 35;
1596  }
1597  else if (higgsType == 3 && idNew == 4) {
1598  nameSave = "c g -> A0(A3) c";
1599  codeSave = 1051;
1600  idRes = 36;
1601  }
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  if (higgsType == 0 && idNew == 5) {
1606  nameSave = "b g -> H b (SM)";
1607  codeSave = 911;
1608  idRes = 25;
1609  }
1610  else if (higgsType == 1 && idNew == 5) {
1611  nameSave = "b g -> h0(H1) b";
1612  codeSave = 1011;
1613  idRes = 25;
1614  }
1615  else if (higgsType == 2 && idNew == 5) {
1616  nameSave = "b g -> H0(H2) b";
1617  codeSave = 1031;
1618  idRes = 35;
1619  }
1620  else if (higgsType == 3 && idNew == 5) {
1621  nameSave = "b g -> A0(A3) b";
1622  codeSave = 1051;
1623  idRes = 36;
1624  }
1625 
1626  // Standard parameters.
1627  m2W = pow2( particleDataPtr->m0(24) );
1628  thetaWRat = 1. / (24. * couplingsPtr->sin2thetaW());
1629 
1630  // Secondary open width fraction.
1631  openFrac = particleDataPtr->resOpenFrac(idRes);
1632 
1633 
1634 }
1635 
1636 //--------------------------------------------------------------------------
1637 
1638 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
1639 
1640 void Sigma2qg2Hq::sigmaKin() {
1641 
1642  // Running mass provides coupling.
1643  double m2Run = pow2( particleDataPtr->mRun(idNew, mH) );
1644 
1645  // Cross section, including couplings and kinematics.
1646  sigma = (M_PI / sH2) * alpS * alpEM * thetaWRat * (m2Run/m2W)
1647  * (sH / (s4 - uH) + 2. * s4 * (s3 - uH) / pow2(s4 - uH)
1648  + (s4 - uH) / sH - 2. * s4 / (s4 - uH)
1649  + 2. * (s3 - uH) * (s3 - s4 - sH) / ((s4 - uH) * sH) );
1650 
1651  // Include secondary width for H0, H1, H2 or A3. Done.
1652  sigma *= openFrac;
1653 
1654 }
1655 
1656 //--------------------------------------------------------------------------
1657 
1658 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
1659 
1660 double Sigma2qg2Hq::sigmaHat() {
1661 
1662  // Check that specified flavour present.
1663  if (abs(id1) != idNew && abs(id2) != idNew) return 0.;
1664 
1665  // Answer.
1666  return sigma;
1667 
1668 }
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  int idq = (id2 == 21) ? id1 : id2;
1679  setId( id1, id2, idRes, idq);
1680 
1681  // tH defined between f and f': must swap tHat <-> uHat if q g in.
1682  swapTU = (id2 == 21);
1683 
1684  // Colour flow topologies. Swap when antiquarks.
1685  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
1686  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
1687  if (idq < 0) swapColAcol();
1688 
1689 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
1700 
1701  // For Higgs decay hand over to standard routine.
1702  if (idMother == 25 || idMother == 35 || idMother == 36)
1703  return weightHiggsDecay( process, iResBeg, iResEnd);
1704 
1705  // For top decay hand over to standard routine.
1706  if (idMother == 6)
1707  return weightTopDecay( process, iResBeg, iResEnd);
1708 
1709  // Else done.
1710  return 1.;
1711 
1712 }
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  if (higgsType == 0) {
1728  nameSave = "g g -> H g (SM; top loop)";
1729  codeSave = 914;
1730  idRes = 25;
1731  }
1732  else if (higgsType == 1) {
1733  nameSave = "g g -> h0(H1) g (BSM; top loop)";
1734  codeSave = 1014;
1735  idRes = 25;
1736  }
1737  else if (higgsType == 2) {
1738  nameSave = "g g -> H0(H2) g (BSM; top loop)";
1739  codeSave = 1034;
1740  idRes = 35;
1741  }
1742  else if (higgsType == 3) {
1743  nameSave = "g g -> A0(A3) g (BSM; top loop)";
1744  codeSave = 1054;
1745  idRes = 36;
1746  }
1747 
1748  // Normalization factor by g g -> H partial width.
1749  // (H can be H0 SM or H1, H2, A3 from BSM).
1750  double mHiggs = particleDataPtr->m0(idRes);
1751  widHgg = particleDataPtr->resWidthChan(idRes, mHiggs, 21, 21);
1752 
1753  // Secondary open width fraction.
1754  openFrac = particleDataPtr->resOpenFrac(idRes);
1755 
1756 }
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  sigma = (M_PI / sH2) * (3. / 16.) * alpS * (widHgg / m3)
1766  * (sH2 * sH2 + tH2 * tH2 + uH2 * uH2 + pow4(s3))
1767  / (sH * tH * uH * s3);
1768  sigma *= openFrac;
1769 
1770 }
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  setId( 21, 21, idRes, 21);
1781 
1782  // Colour flow topologies: random choice between two mirrors.
1783  if (rndmPtr->flat() < 0.5) setColAcol( 1, 2, 2, 3, 0, 0, 1, 3);
1784  else setColAcol( 1, 2, 3, 1, 0, 0, 3, 2);
1785 
1786 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
1797 
1798  // For Higgs decay hand over to standard routine.
1799  if (idMother == 25 || idMother == 35 || idMother == 36)
1800  return weightHiggsDecay( process, iResBeg, iResEnd);
1801 
1802  // For top decay hand over to standard routine.
1803  if (idMother == 6)
1804  return weightTopDecay( process, iResBeg, iResEnd);
1805 
1806  // Else done.
1807  return 1.;
1808 
1809 }
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  if (higgsType == 0) {
1825  nameSave = "q g -> H q (SM; top loop)";
1826  codeSave = 915;
1827  idRes = 25;
1828  }
1829  else if (higgsType == 1) {
1830  nameSave = "q g -> h0(H1) q (BSM; top loop)";
1831  codeSave = 1015;
1832  idRes = 25;
1833  }
1834  else if (higgsType == 2) {
1835  nameSave = "q g -> H0(H2) q (BSM; top loop)";
1836  codeSave = 1035;
1837  idRes = 35;
1838  }
1839  else if (higgsType == 3) {
1840  nameSave = "q g -> A0(A3) q (BSM; top loop)";
1841  codeSave = 1055;
1842  idRes = 36;
1843  }
1844 
1845  // Normalization factor by g g -> H partial width.
1846  // (H can be H0 SM or H1, H2, A3 from BSM).
1847  double mHiggs = particleDataPtr->m0(idRes);
1848  widHgg = particleDataPtr->resWidthChan(idRes, mHiggs, 21, 21);
1849 
1850  // Secondary open width fraction.
1851  openFrac = particleDataPtr->resOpenFrac(idRes);
1852 
1853 }
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  sigma = (M_PI / sH2) * (1. / 12.) * alpS * (widHgg / m3)
1863  * (sH2 + uH2) / (-tH * s3);
1864  sigma *= openFrac;
1865 
1866 }
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  int idq = (id2 == 21) ? id1 : id2;
1877  setId( id1, id2, idRes, idq);
1878 
1879  // tH defined between f and f': must swap tHat <-> uHat if q g in.
1880  swapTU = (id2 == 21);
1881 
1882  // Colour flow topologies. Swap when antiquarks.
1883  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
1884  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
1885  if (idq < 0) swapColAcol();
1886 
1887 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
1898 
1899  // For Higgs decay hand over to standard routine.
1900  if (idMother == 25 || idMother == 35 || idMother == 36)
1901  return weightHiggsDecay( process, iResBeg, iResEnd);
1902 
1903  // For top decay hand over to standard routine.
1904  if (idMother == 6)
1905  return weightTopDecay( process, iResBeg, iResEnd);
1906 
1907  // Else done.
1908  return 1.;
1909 
1910 }
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  if (higgsType == 0) {
1926  nameSave = "q qbar -> H g (SM; top loop)";
1927  codeSave = 916;
1928  idRes = 25;
1929  }
1930  else if (higgsType == 1) {
1931  nameSave = "q qbar -> h0(H1) g (BSM; top loop)";
1932  codeSave = 1016;
1933  idRes = 25;
1934  }
1935  else if (higgsType == 2) {
1936  nameSave = "q qbar -> H0(H2) g (BSM; top loop)";
1937  codeSave = 1036;
1938  idRes = 35;
1939  }
1940  else if (higgsType == 3) {
1941  nameSave = "q qbar -> A0(A3) g (BSM; top loop)";
1942  codeSave = 1056;
1943  idRes = 36;
1944  }
1945 
1946  // Normalization factor by g g -> H partial width.
1947  // (H can be H0 SM or H1, H2, A3 from BSM).
1948  double mHiggs = particleDataPtr->m0(idRes);
1949  widHgg = particleDataPtr->resWidthChan(idRes, mHiggs, 21, 21);
1950 
1951  // Secondary open width fraction.
1952  openFrac = particleDataPtr->resOpenFrac(idRes);
1953 
1954 
1955 }
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  sigma = (M_PI / sH2) * (2. / 9.) * alpS * (widHgg / m3)
1965  * (tH2 + uH2) / (sH * s3);
1966  sigma *= openFrac;
1967 
1968 }
1969 
1970 //--------------------------------------------------------------------------
1971 
1972 // Select identity, colour and anticolour.
1973 
1974 void Sigma2qqbar2Hglt::setIdColAcol() {
1975 
1976  // Flavours trivial.
1977  setId( id1, id2, idRes, 21);
1978 
1979  // Colour flow topologies. Swap when antiquarks.
1980  setColAcol( 1, 0, 0, 2, 0, 0, 1, 2);
1981  if (id1 < 0) swapColAcol();
1982 
1983 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
1994 
1995  // For Higgs decay hand over to standard routine.
1996  if (idMother == 25 || idMother == 35 || idMother == 36)
1997  return weightHiggsDecay( process, iResBeg, iResEnd);
1998 
1999  // For top decay hand over to standard routine.
2000  if (idMother == 6)
2001  return weightTopDecay( process, iResBeg, iResEnd);
2002 
2003  // Else done.
2004  return 1.;
2005 
2006 }
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  HResPtr = particleDataPtr->particleDataEntryPtr(37);
2022 
2023  // Store H+- mass and width for propagator.
2024  mRes = HResPtr->m0();
2025  GammaRes = HResPtr->mWidth();
2026  m2Res = mRes*mRes;
2027  GamMRat = GammaRes / mRes;
2028 
2029  // Couplings.
2030  m2W = pow2(particleDataPtr->m0(24));
2031  thetaWRat = 1. / (8. * couplingsPtr->sin2thetaW());
2032  tan2Beta = pow2(settingsPtr->parm("HiggsHchg:tanBeta"));
2033 
2034 }
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  sigBW = 4. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
2044  widthOutPos = HResPtr->resWidthOpen( 37, mH);
2045  widthOutNeg = HResPtr->resWidthOpen(-37, mH);
2046 
2047 }
2048 
2049 //--------------------------------------------------------------------------
2050 
2051 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2052 
2053 double Sigma1ffbar2Hchg::sigmaHat() {
2054 
2055  // Only allow generation-diagonal states.
2056  int id1Abs = abs(id1);
2057  int id2Abs = abs(id2);
2058  int idUp = max(id1Abs, id2Abs);
2059  int idDn = min(id1Abs, id2Abs);
2060  if (idUp%2 != 0 || idUp - idDn != 1) return 0.;
2061 
2062  // Calculate mass-dependent incoming width. Total cross section.
2063  double m2RunUp = pow2(particleDataPtr->mRun(idUp, mH));
2064  double m2RunDn = pow2(particleDataPtr->mRun(idDn, mH));
2065  double widthIn = alpEM * thetaWRat * (mH/m2W)
2066  * (m2RunDn * tan2Beta + m2RunUp / tan2Beta);
2067  int idUpChg = (id1Abs%2 == 0) ? id1 : id2;
2068  double sigma = (idUpChg > 0) ? widthIn * sigBW * widthOutPos
2069  : widthIn * sigBW * widthOutNeg;
2070 
2071  // Colour factor. Answer.
2072  if (idUp < 9) sigma /= 3.;
2073  return sigma;
2074 
2075 }
2076 
2077 //--------------------------------------------------------------------------
2078 
2079 // Select identity, colour and anticolour.
2080 
2081 void Sigma1ffbar2Hchg::setIdColAcol() {
2082 
2083  // Charge of Higgs. Fill flavours.
2084  int idUpChg = (abs(id1)%2 == 0) ? id1 : id2;
2085  int idHchg = (idUpChg > 0) ? 37 : -37;
2086  setId( id1, id2, idHchg);
2087 
2088  // Colour flow topologies. Swap when antiquarks.
2089  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2090  else setColAcol( 0, 0, 0, 0, 0, 0);
2091  if (id1 < 0) swapColAcol();
2092 
2093 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
2104 
2105  // For Higgs decay hand over to standard routine.
2106  if (idMother == 25 || idMother == 35 || idMother == 36)
2107  return weightHiggsDecay( process, iResBeg, iResEnd);
2108 
2109  // For top decay hand over to standard routine.
2110  if (idMother == 6)
2111  return weightTopDecay( process, iResBeg, iResEnd);
2112 
2113  // Else done.
2114  return 1.;
2115 
2116 }
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  m2W = pow2( particleDataPtr->m0(24) );
2131  thetaWRat = 1. / (24. * couplingsPtr->sin2thetaW());
2132  tan2Beta = pow2(settingsPtr->parm("HiggsHchg:tanBeta"));
2133 
2134  // Incoming flavour within same doublet. Uptype and downtype flavours.
2135  idOld = (idNew%2 == 0) ? idNew - 1 : idNew + 1;
2136  idUp = max(idOld, idNew);
2137  idDn = min(idOld, idNew);
2138 
2139  // Secondary open width fraction.
2140  openFracPos = (idOld%2 == 0) ? particleDataPtr->resOpenFrac( 37, idNew)
2141  : particleDataPtr->resOpenFrac(-37, idNew);
2142  openFracNeg = (idOld%2 == 0) ? particleDataPtr->resOpenFrac(-37, -idNew)
2143  : particleDataPtr->resOpenFrac( 37, -idNew);
2144 
2145 }
2146 
2147 //--------------------------------------------------------------------------
2148 
2149 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2150 
2151 void Sigma2qg2Hchgq::sigmaKin() {
2152 
2153  // Running masses provides coupling.
2154  double m2RunUp = pow2(particleDataPtr->mRun(idUp, mH));
2155  double m2RunDn = pow2(particleDataPtr->mRun(idDn, mH));
2156 
2157  // Cross section, including couplings and kinematics.
2158  sigma = (M_PI / sH2) * alpS * alpEM * thetaWRat
2159  * (m2RunDn * tan2Beta + m2RunUp / tan2Beta) / m2W
2160  * (sH / (s4 - uH) + 2. * s4 * (s3 - uH) / pow2(s4 - uH)
2161  + (s4 - uH) / sH - 2. * s4 / (s4 - uH)
2162  + 2. * (s3 - uH) * (s3 - s4 - sH) / ((s4 - uH) * sH) );
2163 
2164 }
2165 
2166 //--------------------------------------------------------------------------
2167 
2168 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2169 
2170 double Sigma2qg2Hchgq::sigmaHat() {
2171 
2172  // Check that specified flavour present.
2173  if (abs(id1) != idOld && abs(id2) != idOld) return 0.;
2174 
2175  // Answer.
2176  return (id1 == idOld || id2 == idOld) ? sigma * openFracPos
2177  : sigma * openFracNeg;
2178 
2179 }
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  int idq = (id2 == 21) ? id1 : id2;
2189  id3 = ( (idq > 0 && idOld%2 == 0) || (idq < 0 && idOld%2 != 0) )
2190  ? 37 : -37;
2191  id4 = (idq > 0) ? idNew : -idNew;
2192  setId( id1, id2, id3, id4);
2193 
2194  // tH defined between f and f': must swap tHat <-> uHat if q g in.
2195  swapTU = (id2 == 21);
2196 
2197  // Colour flow topologies. Swap when antiquarks.
2198  if (id2 == 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
2199  else setColAcol( 2, 1, 1, 0, 0, 0, 2, 0);
2200  if (idq < 0) swapColAcol();
2201 
2202 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
2213 
2214  // For Higgs decay hand over to standard routine.
2215  if (idMother == 25 || idMother == 35 || idMother == 36)
2216  return weightHiggsDecay( process, iResBeg, iResEnd);
2217 
2218  // For top decay hand over to standard routine.
2219  if (idMother == 6)
2220  return weightTopDecay( process, iResBeg, iResEnd);
2221 
2222  // Else done.
2223  return 1.;
2224 
2225 }
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  higgs12 = (higgsType == 1) ? 25 : 35;
2240  codeSave = (higgsType == 1) ? 1081 : 1082;
2241  nameSave = (higgsType == 1) ? "f fbar -> A0(H3) h0(H1)"
2242  : "f fbar -> A0(H3) H0(H2)";
2243  coupZA3H12 = (higgsType == 1) ? settingsPtr->parm("HiggsA3:coup2H1Z")
2244  : settingsPtr->parm("HiggsA3:coup2H2Z");
2245 
2246  // Standard parameters.
2247  double mZ = particleDataPtr->m0(23);
2248  double GammaZ = particleDataPtr->mWidth(23);
2249  m2Z = mZ * mZ;
2250  mGammaZ = mZ * GammaZ;
2251  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW()
2252  * couplingsPtr->cos2thetaW());
2253 
2254  // Secondary open width fraction.
2255  openFrac = particleDataPtr->resOpenFrac(36, higgs12);
2256 
2257 }
2258 
2259 //--------------------------------------------------------------------------
2260 
2261 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2262 
2263 void Sigma2ffbar2A3H12::sigmaKin() {
2264 
2265  // Common kinematics factora.
2266  sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat * coupZA3H12)
2267  * (uH * tH - s3 * s4) / ( pow2(sH - m2Z) + pow2(mGammaZ) );
2268 
2269 }
2270 
2271 //--------------------------------------------------------------------------
2272 
2273 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2274 
2275 double Sigma2ffbar2A3H12::sigmaHat() {
2276 
2277  // Couplings for incoming flavour.
2278  int idAbs = abs(id1);
2279  double lIn = couplingsPtr->lf(idAbs);
2280  double rIn = couplingsPtr->rf(idAbs);
2281 
2282  // Combine to total cross section. Colour factor.
2283  double sigma = (pow2(lIn) + pow2(rIn)) * sigma0 * openFrac;
2284  if (idAbs < 9) sigma /= 3.;
2285  return sigma;
2286 
2287 }
2288 
2289 //--------------------------------------------------------------------------
2290 
2291 // Select identity, colour and anticolour.
2292 
2293 void Sigma2ffbar2A3H12::setIdColAcol() {
2294 
2295  // Flavours trivial
2296  setId( id1, id2, 36, higgs12);
2297 
2298  // Colour flow topologies. Swap when antiquarks.
2299  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2300  else setColAcol( 0, 0, 0, 0, 0, 0);
2301  if (id1 < 0) swapColAcol();
2302 
2303 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
2314 
2315  // For Higgs decay hand over to standard routine.
2316  if (idMother == 25 || idMother == 35 || idMother == 36)
2317  return weightHiggsDecay( process, iResBeg, iResEnd);
2318 
2319  // For top decay hand over to standard routine.
2320  if (idMother == 6)
2321  return weightTopDecay( process, iResBeg, iResEnd);
2322 
2323  // Else done.
2324  return 1.;
2325 
2326 }
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  higgs12 = (higgsType == 1) ? 25 : 35;
2341  codeSave = (higgsType == 1) ? 1083 : 1084;
2342  nameSave = (higgsType == 1) ? "f fbar' -> H+- h0(H1)"
2343  : "f fbar' -> H+- H0(H2)";
2344  coupWHchgH12 = (higgsType == 1) ? settingsPtr->parm("HiggsHchg:coup2H1W")
2345  : settingsPtr->parm("HiggsHchg:coup2H2W");
2346 
2347  // Standard parameters.
2348  double mW = particleDataPtr->m0(24);
2349  double GammaW = particleDataPtr->mWidth(24);
2350  m2W = mW * mW;
2351  mGammaW = mW * GammaW;
2352  thetaWRat = 1. / (2. * couplingsPtr->sin2thetaW());
2353 
2354  // Secondary open width fraction.
2355  openFracPos = particleDataPtr->resOpenFrac( 37, higgs12);
2356  openFracNeg = particleDataPtr->resOpenFrac(-37, higgs12);
2357 
2358 }
2359 
2360 //--------------------------------------------------------------------------
2361 
2362 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2363 
2364 void Sigma2ffbar2HchgH12::sigmaKin() {
2365 
2366  // Common kinematics factora.
2367  sigma0 = 0.5 * (M_PI / sH2) * pow2(alpEM * thetaWRat * coupWHchgH12)
2368  * (uH * tH - s3 * s4) / ( pow2(sH - m2W) + pow2(mGammaW) );
2369 
2370 }
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  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2380  double sigma = (idUp > 0) ? sigma0 * openFracPos : sigma0 * openFracNeg;
2381  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
2382  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  int idUpChg = (abs(id1)%2 == 0) ? id1 : id2;
2394  int idHchg = (idUpChg > 0) ? 37 : -37;
2395  setId( id1, id2, idHchg, higgs12);
2396 
2397  // Colour flow topologies. Swap when antiquarks.
2398  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2399  else setColAcol( 0, 0, 0, 0, 0, 0);
2400  if (id1 < 0) swapColAcol();
2401 
2402 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
2413 
2414  // For Higgs decay hand over to standard routine.
2415  if (idMother == 25 || idMother == 35 || idMother == 36)
2416  return weightHiggsDecay( process, iResBeg, iResEnd);
2417 
2418  // For top decay hand over to standard routine.
2419  if (idMother == 6)
2420  return weightTopDecay( process, iResBeg, iResEnd);
2421 
2422  // Else done.
2423  return 1.;
2424 
2425 }
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  double mZ = particleDataPtr->m0(23);
2440  double GammaZ = particleDataPtr->mWidth(23);
2441  m2Z = mZ * mZ;
2442  mGammaZ = mZ * GammaZ;
2443  thetaWRat = 1. / (4. * couplingsPtr->sin2thetaW()
2444  * couplingsPtr->cos2thetaW());
2445 
2446  // Charged Higgs coupling to gamma and Z0.
2447  eH = -1.;
2448  lH = -1. + 2. * couplingsPtr->sin2thetaW();
2449 
2450  // Secondary open width fraction.
2451  openFrac = particleDataPtr->resOpenFrac(37, -37);
2452 
2453 }
2454 
2455 //--------------------------------------------------------------------------
2456 
2457 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2458 
2459 void Sigma2ffbar2HposHneg::sigmaKin() {
2460 
2461  // Common kinematics factora.
2462  double preFac = M_PI * pow2(alpEM) * ((uH * tH - s3 * s4) / sH2);
2463  double propZ = 1. / ( pow2(sH - m2Z) + pow2(mGammaZ) );
2464 
2465  // Separate parts for gamma*, interference and Z0.
2466  gamSig = preFac * 2. * pow2(eH) / sH2;
2467  intSig = preFac * 2. * eH * lH * thetaWRat * propZ * (sH - m2Z) / sH;
2468  resSig = preFac * pow2(lH * thetaWRat) * propZ;
2469 
2470 }
2471 
2472 //--------------------------------------------------------------------------
2473 
2474 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2475 
2476 double Sigma2ffbar2HposHneg::sigmaHat() {
2477 
2478  // Couplings for incoming flavour.
2479  int idAbs = abs(id1);
2480  double eIn = couplingsPtr->ef(idAbs);
2481  double lIn = couplingsPtr->lf(idAbs);
2482  double rIn = couplingsPtr->rf(idAbs);
2483 
2484  // Combine to total cross section. Colour factor.
2485  double sigma = (pow2(eIn) * gamSig + eIn * (lIn + rIn) * intSig
2486  + (pow2(lIn) + pow2(rIn)) * resSig) * openFrac;
2487  if (idAbs < 9) sigma /= 3.;
2488  return sigma;
2489 
2490 }
2491 
2492 //--------------------------------------------------------------------------
2493 
2494 // Select identity, colour and anticolour.
2495 
2496 void Sigma2ffbar2HposHneg::setIdColAcol() {
2497 
2498  // Flavours trivial
2499  setId( id1, id2, 37, -37);
2500 
2501  // Colour flow topologies. Swap when antiquarks.
2502  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2503  else setColAcol( 0, 0, 0, 0, 0, 0);
2504  if (id1 < 0) swapColAcol();
2505 
2506 }
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  int idMother = process[process[iResBeg].mother1()].idAbs();
2517 
2518  // For Higgs decay hand over to standard routine.
2519  if (idMother == 25 || idMother == 35 || idMother == 36)
2520  return weightHiggsDecay( process, iResBeg, iResEnd);
2521 
2522  // For top decay hand over to standard routine.
2523  if (idMother == 6)
2524  return weightTopDecay( process, iResBeg, iResEnd);
2525 
2526  // Else done.
2527  return 1.;
2528 
2529 }
2530 
2531 //==========================================================================
2532 
2533 } // end namespace Pythia8
Definition: AgUStep.h:26