PYTHIA8/pythia8140/src/SigmaHiggs.cxx

   1 // SigmaHiggs.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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
  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. * coupSMPtr->sin2thetaW()
 378                  * coupSMPtr->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 * coupSMPtr->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( coupSMPtr->lf(idAbs) );
 463   double riS   = pow2( coupSMPtr->rf(idAbs) );
 464   idAbs        = process[i3].idAbs();
 465   double lfS   = pow2( coupSMPtr->lf(idAbs) );
 466   double rfS   = pow2( coupSMPtr->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. * coupSMPtr->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 *= coupSMPtr->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 / (coupSMPtr->sin2thetaW()
 668            * coupSMPtr->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( coupSMPtr->lf(id1Abs) );
 706   double rf1S  = pow2( coupSMPtr->rf(id1Abs) );
 707   double lf2S  = pow2( coupSMPtr->lf(id2Abs) );
 708   double rf2S  = pow2( coupSMPtr->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 / coupSMPtr->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) * coupSMPtr->V2CKMsum(id1Abs)
 846     * coupSMPtr->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 = coupSMPtr->V2CKMpick(id1);
 868   id5 = coupSMPtr->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 / coupSMPtr->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 / coupSMPtr->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. * coupSMPtr->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. * coupSMPtr->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. * coupSMPtr->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. * coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW());
2252
2253   // Secondary open width fraction.
2254   openFrac   = particleDataPtr->resOpenFrac(36, higgs12);
2255
2256 }
2257
2258 //--------------------------------------------------------------------------
2259
2260 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2261
2262 void Sigma2ffbar2A3H12::sigmaKin() {
2263
2264   // Common kinematics factora.
2265   sigma0 = (M_PI / sH2) * pow2(alpEM * thetaWRat * coupZA3H12)
2266     * (uH * tH - s3 * s4) / ( pow2(sH - m2Z) + pow2(mGammaZ) );
2267
2268 }
2269
2270 //--------------------------------------------------------------------------
2271
2272 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2273
2274 double Sigma2ffbar2A3H12::sigmaHat() {
2275
2276   // Couplings for incoming flavour.
2277   int idAbs    = abs(id1);
2278   double lIn   = coupSMPtr->lf(idAbs);
2279   double rIn   = coupSMPtr->rf(idAbs);
2280
2281   // Combine to total cross section. Colour factor.
2282   double sigma = (pow2(lIn) + pow2(rIn)) * sigma0 * openFrac;
2283   if (idAbs < 9) sigma /= 3.;
2284   return sigma;
2285
2286 }
2287
2288 //--------------------------------------------------------------------------
2289
2290 // Select identity, colour and anticolour.
2291
2292 void Sigma2ffbar2A3H12::setIdColAcol() {
2293
2294   // Flavours trivial
2295   setId( id1, id2, 36, higgs12);
2296
2297   // Colour flow topologies. Swap when antiquarks.
2298   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2299   else              setColAcol( 0, 0, 0, 0, 0, 0);
2300   if (id1 < 0) swapColAcol();
2301
2302 }
2303
2304 //--------------------------------------------------------------------------
2305
2306 // Evaluate weight for decay angles.
2307
2308 double Sigma2ffbar2A3H12::weightDecay( Event& process, int iResBeg,
2309   int iResEnd) {
2310
2311   // Identity of mother of decaying reseonance(s).
2312   int idMother = process[process[iResBeg].mother1()].idAbs();
2313
2314   // For Higgs decay hand over to standard routine.
2315   if (idMother == 25 || idMother == 35 || idMother == 36)
2316     return weightHiggsDecay( process, iResBeg, iResEnd);
2317
2318   // For top decay hand over to standard routine.
2319   if (idMother == 6)
2320     return weightTopDecay( process, iResBeg, iResEnd);
2321
2322   // Else done.
2323   return 1.;
2324
2325 }
2326
2327 //==========================================================================
2328
2329 // Sigma2ffbar2HchgH12 class.
2330 // Cross section for f fbar -> H+- h0(H_1) or H+- H0(H_2).
2331
2332 //--------------------------------------------------------------------------
2333
2334 // Initialize process.
2335
2336 void Sigma2ffbar2HchgH12::initProc() {
2337
2338   // Set up whether h0(H_1) or H0(H_2).
2339   higgs12    = (higgsType == 1) ? 25 : 35;
2340   codeSave   = (higgsType == 1) ? 1083 : 1084;
2341   nameSave   = (higgsType == 1) ? "f fbar' -> H+- h0(H1)"
2342                                 : "f fbar' -> H+- H0(H2)";
2343   coupWHchgH12 = (higgsType == 1) ? settingsPtr->parm("HiggsHchg:coup2H1W")
2344                                   : settingsPtr->parm("HiggsHchg:coup2H2W");
2345
2346   // Standard parameters.
2347   double mW  = particleDataPtr->m0(24);
2348   double GammaW = particleDataPtr->mWidth(24);
2349   m2W        = mW * mW;
2350   mGammaW    = mW * GammaW;
2351   thetaWRat  = 1. / (2. * coupSMPtr->sin2thetaW());
2352
2353   // Secondary open width fraction.
2354   openFracPos   = particleDataPtr->resOpenFrac( 37, higgs12);
2355   openFracNeg   = particleDataPtr->resOpenFrac(-37, higgs12);
2356
2357 }
2358
2359 //--------------------------------------------------------------------------
2360
2361 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2362
2363 void Sigma2ffbar2HchgH12::sigmaKin() {
2364
2365   // Common kinematics factora.
2366   sigma0 = 0.5 * (M_PI / sH2) * pow2(alpEM * thetaWRat * coupWHchgH12)
2367     * (uH * tH - s3 * s4) / ( pow2(sH - m2W) + pow2(mGammaW) );
2368
2369 }
2370
2371 //--------------------------------------------------------------------------
2372
2373 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2374
2375 double Sigma2ffbar2HchgH12::sigmaHat() {
2376
2377   // Combine to total cross section. CKM and colour factor.
2378   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
2379   double sigma = (idUp > 0) ? sigma0 * openFracPos : sigma0 * openFracNeg;
2380   if (abs(id1) < 9) sigma *= coupSMPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
2381   return sigma;
2382
2383 }
2384
2385 //--------------------------------------------------------------------------
2386
2387 // Select identity, colour and anticolour.
2388
2389 void Sigma2ffbar2HchgH12::setIdColAcol() {
2390
2391   // Charge of Higgs. Fill flavours.
2392   int idUpChg = (abs(id1)%2 == 0) ? id1 : id2;
2393   int idHchg  = (idUpChg > 0) ? 37 : -37;
2394   setId( id1, id2, idHchg, higgs12);
2395
2396   // Colour flow topologies. Swap when antiquarks.
2397   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2398   else              setColAcol( 0, 0, 0, 0, 0, 0);
2399   if (id1 < 0) swapColAcol();
2400
2401 }
2402
2403 //--------------------------------------------------------------------------
2404
2405 // Evaluate weight for decay angles.
2406
2407 double Sigma2ffbar2HchgH12::weightDecay( Event& process, int iResBeg,
2408   int iResEnd) {
2409
2410   // Identity of mother of decaying reseonance(s).
2411   int idMother = process[process[iResBeg].mother1()].idAbs();
2412
2413   // For Higgs decay hand over to standard routine.
2414   if (idMother == 25 || idMother == 35 || idMother == 36)
2415     return weightHiggsDecay( process, iResBeg, iResEnd);
2416
2417   // For top decay hand over to standard routine.
2418   if (idMother == 6)
2419     return weightTopDecay( process, iResBeg, iResEnd);
2420
2421   // Else done.
2422   return 1.;
2423
2424 }
2425
2426 //==========================================================================
2427
2428 // Sigma2ffbar2HposHneg class.
2429 // Cross section for q g -> H+- q'.
2430
2431 //--------------------------------------------------------------------------
2432
2433 // Initialize process.
2434
2435 void Sigma2ffbar2HposHneg::initProc() {
2436
2437   // Standard parameters.
2438   double mZ = particleDataPtr->m0(23);
2439   double GammaZ = particleDataPtr->mWidth(23);
2440   m2Z       = mZ * mZ;
2441   mGammaZ   = mZ * GammaZ;
2442   thetaWRat = 1. / (4. * coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW());
2443
2444   // Charged Higgs coupling to gamma and Z0.
2445   eH        = -1.;
2446   lH        = -1. + 2. * coupSMPtr->sin2thetaW();
2447
2448   // Secondary open width fraction.
2449   openFrac  = particleDataPtr->resOpenFrac(37, -37);
2450
2451 }
2452
2453 //--------------------------------------------------------------------------
2454
2455 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
2456
2457 void Sigma2ffbar2HposHneg::sigmaKin() {
2458
2459   // Common kinematics factora.
2460   double preFac = M_PI * pow2(alpEM) * ((uH * tH - s3 * s4) / sH2);
2461   double propZ  = 1. / ( pow2(sH - m2Z) + pow2(mGammaZ) );
2462
2463   // Separate parts for gamma*, interference and Z0.
2464   gamSig    = preFac * 2. * pow2(eH) / sH2;
2465   intSig    = preFac * 2. * eH * lH * thetaWRat * propZ * (sH - m2Z) / sH;
2466   resSig    = preFac * pow2(lH * thetaWRat) * propZ;
2467
2468 }
2469
2470 //--------------------------------------------------------------------------
2471
2472 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
2473
2474 double Sigma2ffbar2HposHneg::sigmaHat() {
2475
2476   // Couplings for incoming flavour.
2477   int idAbs    = int(id1);
2478   double eIn   = coupSMPtr->ef(idAbs);
2479   double lIn   = coupSMPtr->lf(idAbs);
2480   double rIn   = coupSMPtr->rf(idAbs);
2481
2482   // Combine to total cross section. Colour factor.
2483   double sigma = (pow2(eIn) * gamSig + eIn * (lIn + rIn) * intSig
2484     + (pow2(lIn) + pow2(rIn)) * resSig) * openFrac;
2485   if (idAbs < 9) sigma /= 3.;
2486   return sigma;
2487
2488 }
2489
2490 //--------------------------------------------------------------------------
2491
2492 // Select identity, colour and anticolour.
2493
2494 void Sigma2ffbar2HposHneg::setIdColAcol() {
2495
2496   // Flavours trivial
2497   setId( id1, id2, 37, -37);
2498
2499   // Colour flow topologies. Swap when antiquarks.
2500   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
2501   else              setColAcol( 0, 0, 0, 0, 0, 0);
2502   if (id1 < 0) swapColAcol();
2503
2504 }
2505
2506 //--------------------------------------------------------------------------
2507
2508 // Evaluate weight for decay angles.
2509
2510 double Sigma2ffbar2HposHneg::weightDecay( Event& process, int iResBeg,
2511   int iResEnd) {
2512
2513   // Identity of mother of decaying reseonance(s).
2514   int idMother = process[process[iResBeg].mother1()].idAbs();
2515
2516   // For Higgs decay hand over to standard routine.
2517   if (idMother == 25 || idMother == 35 || idMother == 36)
2518     return weightHiggsDecay( process, iResBeg, iResEnd);
2519
2520   // For top decay hand over to standard routine.
2521   if (idMother == 6)
2522     return weightTopDecay( process, iResBeg, iResEnd);
2523
2524   // Else done.
2525   return 1.;
2526
2527 }
2528
2529 //==========================================================================
2530
2531 } // end namespace Pythia8