PYTHIA8/pythia8145/src/SigmaLeftRightSym.cxx

   1 // SigmaLeftRightSym.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the
   7 // left-right-symmetry simulation classes.
   8
   9 #include "SigmaLeftRightSym.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // Sigma1ffbar2ZRight class.
  16 // Cross section for f fbar -> Z_R^0 (righthanded gauge boson).
  17
  18 //--------------------------------------------------------------------------
  19
  20 // Initialize process.
  21
  22 void Sigma1ffbar2ZRight::initProc() {
  23
  24   // Store Z_R mass and width for propagator.
  25   idZR     = 9900023;
  26   mRes     = particleDataPtr->m0(idZR);
  27   GammaRes = particleDataPtr->mWidth(idZR);
  28   m2Res    = mRes*mRes;
  29   GamMRat  = GammaRes / mRes;
  30   sin2tW   = couplingsPtr->sin2thetaW();
  31
  32   // Set pointer to particle properties and decay table.
  33   ZRPtr    = particleDataPtr->particleDataEntryPtr(idZR);
  34
  35 }
  36
  37 //--------------------------------------------------------------------------
  38
  39 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
  40
  41 void Sigma1ffbar2ZRight::sigmaKin() {
  42
  43   // Set up Breit-Wigner. Width out only includes open channels.
  44   double sigBW    = 12. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
  45   double widthOut = ZRPtr->resWidthOpen(idZR, mH);
  46
  47   // Prefactor for incoming widths. Combine. Done.
  48   double preFac   = alpEM * mH / ( 48. * sin2tW * (1. - sin2tW)
  49                   * (1. - 2. * sin2tW) );
  50   sigma0          = preFac * sigBW * widthOut;
  51
  52 }
  53
  54 //--------------------------------------------------------------------------
  55
  56 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
  57
  58 double Sigma1ffbar2ZRight::sigmaHat() {
  59
  60   // Vector and axial couplings of incoming fermion pair.
  61   int    idAbs = abs(id1);
  62   double af = 0.;
  63   double vf = 0.;
  64   if (idAbs < 9 && idAbs%2 == 1) {
  65     af      = -1. + 2. * sin2tW;
  66     vf      = -1. + 4. * sin2tW / 3.;
  67   } else if (idAbs < 9) {
  68     af      = 1. - 2. * sin2tW;
  69     vf      = 1. - 8. * sin2tW / 3.;
  70   } else if (idAbs < 19 && idAbs%2 == 1) {
  71     af      = -1. + 2. * sin2tW;
  72     vf      = -1. + 4. * sin2tW;
  73   }
  74
  75   // Colour factor. Answer.
  76   double sigma = (vf*vf + af*af) * sigma0;
  77   if (idAbs < 9) sigma /= 3.;
  78   return sigma;
  79
  80 }
  81
  82 //--------------------------------------------------------------------------
  83
  84 // Select identity, colour and anticolour.
  85
  86 void Sigma1ffbar2ZRight::setIdColAcol() {
  87
  88   // Flavours trivial.
  89   setId( id1, id2, idZR);
  90
  91   // Colour flow topologies. Swap when antiquarks.
  92   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
  93   else              setColAcol( 0, 0, 0, 0, 0, 0);
  94   if (id1 < 0) swapColAcol();
  95
  96 }
  97
  98 //--------------------------------------------------------------------------
  99
 100 // Evaluate weight for Z_R decay angle.
 101
 102 double Sigma1ffbar2ZRight::weightDecay( Event& process, int iResBeg,
 103   int iResEnd) {
 104
 105   // Identity of mother of decaying reseonance(s).
 106   int idMother = process[process[iResBeg].mother1()].idAbs();
 107
 108   // For top decay hand over to standard routine.
 109   if (idMother == 6)
 110     return weightTopDecay( process, iResBeg, iResEnd);
 111
 112   // Z_R should sit in entry 5.
 113   if (iResBeg != 5 || iResEnd != 5) return 1.;
 114
 115   // Couplings for in- and out-flavours.
 116   double ai, vi, af, vf;
 117   int idInAbs   = process[3].idAbs();
 118   if (idInAbs < 9 && idInAbs%2 == 1) {
 119     ai          = -1. + 2. * sin2tW;
 120     vi          = -1. + 4. * sin2tW / 3.;
 121   } else if (idInAbs < 9) {
 122     ai          = 1. - 2. * sin2tW;
 123     vi          = 1. - 8. * sin2tW / 3.;
 124   } else {
 125     ai          = -1. + 2. * sin2tW;
 126     vi          = -1. + 4. * sin2tW;
 127   }
 128   int idOutAbs = process[6].idAbs();
 129   if (idOutAbs < 9 && idOutAbs%2 == 1) {
 130     af          = -1. + 2. * sin2tW;
 131     vf          = -1. + 4. * sin2tW / 3.;
 132   } else if (idOutAbs < 9) {
 133     af          = 1. - 2. * sin2tW;
 134     vf          = 1. - 8. * sin2tW / 3.;
 135   } else {
 136     af          = -1. + 2. * sin2tW;
 137     vf          = -1. + 4. * sin2tW;
 138   }
 139
 140   // Phase space factors. Reconstruct decay angle.
 141   double mr1    = pow2(process[6].m()) / sH;
 142   double mr2    = pow2(process[7].m()) / sH;
 143   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 144   double cosThe = (process[3].p() - process[4].p())
 145     * (process[7].p() - process[6].p()) / (sH * betaf);
 146
 147   // Angular weight and its maximum.
 148   double wt1    = (vi*vi + ai*ai) * (vf*vf + af*af * betaf*betaf);
 149   double wt2    = (1. - betaf*betaf) * (vi*vi + ai*ai) * vf*vf;
 150   double wt3    = betaf * 4. * vi * ai * vf * af;
 151   if (process[3].id() * process[6].id() < 0) wt3 = -wt3;
 152   double wt     = wt1 * (1. + cosThe*cosThe) + wt2 * (1. - cosThe*cosThe)
 153                 + 2. * wt3 * cosThe;
 154   double wtMax  = 2. * (wt1 + abs(wt3));
 155
 156   // Done.
 157   return wt / wtMax;
 158
 159 }
 160
 161 //==========================================================================
 162
 163 // Sigma1ffbar2WRight class.
 164 // Cross section for f fbar' -> W_R^+- (righthanded gauge boson).
 165
 166 //--------------------------------------------------------------------------
 167
 168 // Initialize process.
 169
 170 void Sigma1ffbar2WRight::initProc() {
 171
 172   // Store W_R^+- mass and width for propagator.
 173   idWR     = 9900024;
 174   mRes     = particleDataPtr->m0(idWR);
 175   GammaRes = particleDataPtr->mWidth(idWR);
 176   m2Res    = mRes*mRes;
 177   GamMRat  = GammaRes / mRes;
 178   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
 179
 180   // Set pointer to particle properties and decay table.
 181   particlePtr = particleDataPtr->particleDataEntryPtr(idWR);
 182
 183 }
 184
 185 //--------------------------------------------------------------------------
 186
 187 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 188
 189 void Sigma1ffbar2WRight::sigmaKin() {
 190
 191   // Common coupling factors.
 192   double colQ   = 3. * (1. + alpS / M_PI);
 193
 194   // Reset quantities to sum. Declare variables inside loop.
 195   double widOutPos = 0.;
 196   double widOutNeg = 0.;
 197   int    id1Now, id2Now, id1Abs, id2Abs, id1Neg, id2Neg, onMode;
 198   double widNow, widSecPos, widSecNeg, mf1, mf2, mr1, mr2, kinFac;
 199
 200   // Loop over all W_R^+- decay channels.
 201   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
 202     widNow = 0.;
 203     id1Now      = particlePtr->channel(i).product(0);
 204     id2Now      = particlePtr->channel(i).product(1);
 205     id1Abs      = abs(id1Now);
 206     id2Abs      = abs(id2Now);
 207
 208     // Check that above threshold. Phase space.
 209     mf1 = particleDataPtr->m0(id1Abs);
 210     mf2 = particleDataPtr->m0(id2Abs);
 211     if (mH > mf1 + mf2 + MASSMARGIN) {
 212       mr1    = pow2(mf1 / mH);
 213       mr2    = pow2(mf2 / mH);
 214       kinFac = (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
 215              * sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
 216
 217       // Combine kinematics with colour factor and CKM couplings.
 218       widNow = kinFac;
 219       if (id1Abs < 9) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
 220
 221       // Secondary width from top and righthanded neutrino decay.
 222       id1Neg    = (id1Abs < 19) ? -id1Now : id1Abs;
 223       id2Neg    = (id2Abs < 19) ? -id2Now : id2Abs;
 224       widSecPos = particleDataPtr->resOpenFrac(id1Now, id2Now);
 225       widSecNeg = particleDataPtr->resOpenFrac(id1Neg, id2Neg);
 226
 227       // Add weight for channels on for all, W_R^+ and W_R^-, respectively.
 228       onMode = particlePtr->channel(i).onMode();
 229       if (onMode == 1 || onMode == 2) widOutPos += widNow * widSecPos;
 230       if (onMode == 1 || onMode == 3) widOutNeg += widNow * widSecNeg;
 231
 232     // End loop over fermions.
 233     }
 234   }
 235
 236   // Set up Breit-Wigner. Cross section for W_R^+ and W_R^- separately.
 237   double sigBW = 12. * M_PI * pow2(alpEM * thetaWRat) * sH
 238                / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 239   sigma0Pos = sigBW * widOutPos;
 240   sigma0Neg = sigBW * widOutNeg;
 241
 242 }
 243
 244 //--------------------------------------------------------------------------
 245
 246 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 247
 248 double Sigma1ffbar2WRight::sigmaHat() {
 249
 250   // Secondary width for W_R^+ or W_R^-. CKM and colour factors.
 251   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
 252   double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
 253   if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
 254
 255   // Answer.
 256   return sigma;
 257
 258 }
 259
 260 //--------------------------------------------------------------------------
 261
 262 // Select identity, colour and anticolour.
 263
 264 void Sigma1ffbar2WRight::setIdColAcol() {
 265
 266   // Sign of outgoing W_R.
 267   int sign          = (abs(id1)%2 == 0) ? 1 : -1;
 268   if (id1 < 0) sign = -sign;
 269   setId( id1, id2, idWR * sign);
 270
 271   // Colour flow topologies. Swap when antiquarks.
 272   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
 273   else              setColAcol( 0, 0, 0, 0, 0, 0);
 274   if (id1 < 0) swapColAcol();
 275
 276 }
 277
 278 //--------------------------------------------------------------------------
 279
 280 // Evaluate weight for W_R decay angle.
 281
 282 double Sigma1ffbar2WRight::weightDecay( Event& process, int iResBeg,
 283   int iResEnd) {
 284
 285   // Identity of mother of decaying reseonance(s).
 286   int idMother = process[process[iResBeg].mother1()].idAbs();
 287
 288   // For top decay hand over to standard routine.
 289   if (idMother == 6)
 290     return weightTopDecay( process, iResBeg, iResEnd);
 291
 292   // W_R should sit in entry 5.
 293   if (iResBeg != 5 || iResEnd != 5) return 1.;
 294
 295   // Phase space factors.
 296   double mr1    = pow2(process[6].m()) / sH;
 297   double mr2    = pow2(process[7].m()) / sH;
 298   double betaf  = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 299
 300   // Sign of asymmetry.
 301   double eps    = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
 302
 303   // Reconstruct decay angle and weight for it.
 304   double cosThe = (process[3].p() - process[4].p())
 305     * (process[7].p() - process[6].p()) / (sH * betaf);
 306   double wtMax  = 4.;
 307   double wt     = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
 308
 309   // Done.
 310   return (wt / wtMax);
 311
 312 }
 313
 314 //==========================================================================
 315
 316 // Sigma1ll2Hchgchg class.
 317 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
 318
 319 //--------------------------------------------------------------------------
 320
 321 // Initialize process.
 322
 323 void Sigma1ll2Hchgchg::initProc() {
 324
 325   // Set process properties: H_L^++-- or H_R^++--.
 326   if (leftRight == 1) {
 327     idHLR    = 9900041;
 328     codeSave = 3121;
 329     nameSave = "l l -> H_L^++--";
 330   } else {
 331     idHLR    = 9900042;
 332     codeSave = 3141;
 333     nameSave = "l l -> H_R^++--";
 334   }
 335
 336   // Read in Yukawa matrix for couplings to a lepton pair.
 337   yukawa[1][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHee");
 338   yukawa[2][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
 339   yukawa[2][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
 340   yukawa[3][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
 341   yukawa[3][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
 342   yukawa[3][3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
 343
 344   // Store H_L/R mass and width for propagator.
 345   mRes     = particleDataPtr->m0(idHLR);
 346   GammaRes = particleDataPtr->mWidth(idHLR);
 347   m2Res    = mRes*mRes;
 348   GamMRat  = GammaRes / mRes;
 349
 350   // Set pointer to particle properties and decay table.
 351   particlePtr = particleDataPtr->particleDataEntryPtr(idHLR);
 352
 353 }
 354
 355 //--------------------------------------------------------------------------
 356
 357 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 358
 359 double Sigma1ll2Hchgchg::sigmaHat() {
 360
 361   // Initial state must consist of two identical-sign leptons.
 362   if (id1 * id2 < 0) return 0.;
 363   int id1Abs = abs(id1);
 364   int id2Abs = abs(id2);
 365   if (id1Abs != 11 && id1Abs != 13 && id1Abs != 15) return 0.;
 366   if (id2Abs != 11 && id2Abs != 13 && id2Abs != 15) return 0.;
 367
 368   // Set up Breit-Wigner, inwidth and outwidth.
 369   double sigBW  = 8. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 370   double widIn  = pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2])
 371                 * mH / (8. * M_PI);
 372   int idSgn     = (id1 < 0) ? idHLR : -idHLR;
 373   double widOut = particlePtr->resWidthOpen( idSgn, mH);
 374
 375   // Answer.
 376   return widIn * sigBW * widOut;
 377
 378 }
 379
 380 //--------------------------------------------------------------------------
 381
 382 // Select identity, colour and anticolour.
 383
 384 void Sigma1ll2Hchgchg::setIdColAcol() {
 385
 386   // Sign of outgoing H_L/R.
 387   int idSgn     = (id1 < 0) ? idHLR : -idHLR;
 388   setId( id1, id2, idSgn);
 389
 390   // No colours whatsoever.
 391   setColAcol( 0, 0, 0, 0, 0, 0);
 392
 393 }
 394
 395 //--------------------------------------------------------------------------
 396
 397 // Evaluate weight for H_L/R sequential decay angles.
 398
 399 double Sigma1ll2Hchgchg::weightDecay( Event& process, int iResBeg,
 400   int iResEnd) {
 401
 402   // Identity of mother of decaying reseonance(s).
 403   int idMother = process[process[iResBeg].mother1()].idAbs();
 404
 405   // For top decay hand over to standard routine.
 406   if (idMother == 6)
 407     return weightTopDecay( process, iResBeg, iResEnd);
 408
 409   // Else isotropic decay.
 410   return 1.;
 411
 412 }
 413
 414 //==========================================================================
 415
 416 // Sigma2lgm2Hchgchgl class.
 417 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
 418
 419 //--------------------------------------------------------------------------
 420
 421 // Initialize process.
 422
 423 void Sigma2lgm2Hchgchgl::initProc() {
 424
 425   // Set process properties: H_L^++-- or H_R^++-- and e/mu/tau.
 426   idHLR        = (leftRight == 1) ? 9900041 : 9900042;
 427   codeSave     = (leftRight == 1) ? 3122 : 3142;
 428   if (idLep == 13) codeSave += 2;
 429   if (idLep == 15) codeSave += 4;
 430   if      (codeSave == 3122) nameSave = "l^+- gamma -> H_L^++-- e^-+";
 431   else if (codeSave == 3123) nameSave = "l^+- gamma -> H_L^++-- mu^-+";
 432   else if (codeSave == 3124) nameSave = "l^+- gamma -> H_L^++-- tau^-+";
 433   else if (codeSave == 3142) nameSave = "l^+- gamma -> H_R^++-- e^-+";
 434   else if (codeSave == 3143) nameSave = "l^+- gamma -> H_R^++-- mu^-+";
 435   else                       nameSave = "l^+- gamma -> H_R^++-- tau^-+";
 436
 437   // Read in relevantYukawa matrix for couplings to a lepton pair.
 438   if (idLep == 11) {
 439     yukawa[1]  = settingsPtr->parm("LeftRightSymmmetry:coupHee");
 440     yukawa[2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
 441     yukawa[3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
 442   } else if (idLep == 13) {
 443     yukawa[1]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
 444     yukawa[2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
 445     yukawa[3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
 446   } else {
 447     yukawa[1]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
 448     yukawa[2]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
 449     yukawa[3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
 450   }
 451
 452   // Secondary open width fractions.
 453   openFracPos  = particleDataPtr->resOpenFrac( idHLR);
 454   openFracNeg  = particleDataPtr->resOpenFrac(-idHLR);
 455
 456 }
 457
 458 //--------------------------------------------------------------------------
 459
 460 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 461
 462 double Sigma2lgm2Hchgchgl::sigmaHat() {
 463
 464   // Initial state must consist of a lepton and a photon.
 465   int idIn     = (id2 == 22) ? id1 : id2;
 466   int idInAbs  = abs(idIn);
 467   if (idInAbs != 11 && idInAbs != 13 && idInAbs != 15) return 0.;
 468
 469   // Incoming squared lepton mass.
 470   double s1    = pow2( particleDataPtr->m0(idInAbs) );
 471
 472   // Kinematical expressions.
 473   double smm1  = 8. * (sH + tH - s3) * (sH + tH - 2. * s3 - s1 - s4)
 474                / pow2(uH - s3);
 475   double smm2  = 2. * ( (2. * s3 - 3. * s1) * s4 + (s1 - 2. * s4) * tH
 476                - (tH - s4) * sH ) / pow2(tH - s4);
 477   double smm3  = 2. * ( (2. * s3 - 3. * s4 + tH) * s1
 478                - (2. * s1 - s4 + tH) * sH ) / pow2(sH - s1);
 479   double smm12 = 4. * ( (2. * s1 - s4 - 2. * s3 + tH) * sH
 480                + (tH - 3. * s3 - 3. * s4) * tH + (2. * s3 - 2. * s1
 481                + 3. * s4) * s3 ) / ( (uH - s3) * (tH - s4) );
 482   double smm13 = -4. * ( (tH + s1 - 2. * s4) * tH - (s3 + 3. * s1 - 2. * s4)
 483                * s3 + (s3 + 3. * s1 + tH) * sH - pow2(tH - s3 + sH) )
 484                / ( (uH - s3) * (sH - s1) );
 485   double smm23 = -4. * ( (s1 - s4 + s3) * tH - s3*s3 + s3 * (s1 + s4)
 486                - 3. * s1 * s4 - (s1 - s4 - s3 + tH) * sH)
 487                / ( (sH - s1) * (tH - s4) );
 488   double sigma = alpEM * pow2(sH / (sH - s1) ) * (smm1 + smm2 + smm3
 489                + smm12 + smm13 + smm23) / (4. * sH2);
 490
 491   // Lepton Yukawa and secondary widths.
 492   sigma       *= pow2(yukawa[(idInAbs-9)/2]);
 493   sigma       *= (idIn < 0) ? openFracPos : openFracNeg;
 494
 495   // Answer.
 496   return sigma;
 497
 498 }
 499
 500 //--------------------------------------------------------------------------
 501
 502 // Select identity, colour and anticolour.
 503
 504 void Sigma2lgm2Hchgchgl::setIdColAcol() {
 505
 506   // Sign of outgoing H_L/R.
 507   int idIn     = (id2 == 22) ? id1 : id2;
 508   int idSgn    = (idIn < 0) ? idHLR : -idHLR;
 509   int idOut    = (idIn < 0) ? idLep : -idLep;
 510   setId( id1, id2, idSgn, idOut);
 511
 512   // tHat is defined between incoming lepton and outgoing Higgs.
 513   if (id1 == 22) swapTU = true;
 514
 515   // No colours whatsoever.
 516   setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 517
 518 }
 519
 520 //--------------------------------------------------------------------------
 521
 522 // Evaluate weight for H_L/R sequential decay angles.
 523
 524 double Sigma2lgm2Hchgchgl::weightDecay( Event& process, int iResBeg,
 525   int iResEnd) {
 526
 527   // Identity of mother of decaying reseonance(s).
 528   int idMother = process[process[iResBeg].mother1()].idAbs();
 529
 530   // For top decay hand over to standard routine.
 531   if (idMother == 6)
 532     return weightTopDecay( process, iResBeg, iResEnd);
 533
 534   // Else isotropic decay.
 535   return 1.;
 536
 537 }
 538
 539 //==========================================================================
 540
 541 // Sigma3ff2HchgchgfftWW class.
 542 // Cross section for  f_1 f_2 -> H_(L/R)^++-- f_3 f_4 (W+- W+- fusion).
 543
 544 //--------------------------------------------------------------------------
 545
 546 // Initialize process.
 547
 548 void Sigma3ff2HchgchgfftWW::initProc() {
 549
 550   // Set process properties: H_L^++-- or H_R^++--.
 551   if (leftRight == 1) {
 552     idHLR      = 9900041;
 553     codeSave   = 3125;
 554     nameSave   = "f_1 f_2 -> H_L^++-- f_3 f_4 (W+- W+- fusion)";
 555   } else {
 556     idHLR      = 9900042;
 557     codeSave   = 3145;
 558     nameSave   = "f_1 f_2 -> H_R^++-- f_3 f_4 (W+- W+- fusion)";
 559   }
 560
 561   // Common fixed mass and coupling factor.
 562   double mW    = particleDataPtr->m0(24);
 563   double mWR   = particleDataPtr->m0(9900024);
 564   mWS          = (leftRight == 1) ? pow2(mW) : pow2(mWR);
 565   double gL    = settingsPtr->parm("LeftRightSymmmetry:gL");
 566   double gR    = settingsPtr->parm("LeftRightSymmmetry:gR");
 567   double vL    = settingsPtr->parm("LeftRightSymmmetry:vL");
 568   prefac       = (leftRight == 1) ? pow2(pow4(gL) * vL)
 569                                   : 2. * pow2(pow3(gR) * mWR);
 570   // Secondary open width fractions.
 571   openFracPos  = particleDataPtr->resOpenFrac( idHLR);
 572   openFracNeg  = particleDataPtr->resOpenFrac(-idHLR);
 573
 574 }
 575
 576 //--------------------------------------------------------------------------
 577
 578 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 579
 580 void Sigma3ff2HchgchgfftWW::sigmaKin() {
 581
 582   // Required four-vector products.
 583   double pp12  = 0.5 * sH;
 584   double pp14  = 0.5 * mH * p4cm.pNeg();
 585   double pp15  = 0.5 * mH * p5cm.pNeg();
 586   double pp24  = 0.5 * mH * p4cm.pPos();
 587   double pp25  = 0.5 * mH * p5cm.pPos();
 588   double pp45  = p4cm * p5cm;
 589
 590   // Cross section: kinematics part. Combine with couplings.
 591   double propT = 1. / ( (2. * pp14 + mWS) * (2. * pp25 + mWS) );
 592   double propU = 1. / ( (2. * pp24 + mWS) * (2. * pp15 + mWS) );
 593   sigma0TU     = prefac * pp12 * pp45 * pow2(propT + propU);
 594   sigma0T      = prefac * pp12 * pp45 * 2. * pow2(propT);
 595
 596 }
 597
 598 //--------------------------------------------------------------------------
 599
 600 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 601
 602 double Sigma3ff2HchgchgfftWW::sigmaHat() {
 603
 604   // Do not allow creation of righthanded neutrinos for H_R.
 605   int id1Abs   = abs(id1);
 606   int id2Abs   = abs(id2);
 607   if ( leftRight == 2 && (id1Abs > 10 || id2Abs > 10) ) return 0.;
 608
 609   // Many flavour combinations not possible because of charge.
 610   int chg1     = (( id1Abs%2 == 0 && id1 > 0)
 611                || (id1Abs%2 == 1 && id1 < 0) ) ? 1 : -1;
 612   int chg2     = (( id2Abs%2 == 0 && id2 > 0)
 613                || (id2Abs%2 == 1 && id2 < 0) ) ? 1 : -1;
 614   if (abs(chg1 + chg2) != 2) return 0.;
 615
 616   // Basic cross section. CKM factors for final states.
 617   double sigma = (id2 == id1 && id1Abs > 10) ? sigma0TU : sigma0T;
 618   sigma       *= couplingsPtr->V2CKMsum(id1Abs) * couplingsPtr->V2CKMsum(id2Abs);
 619
 620   // Secondary width for H0.
 621   sigma       *= (chg1 + chg2 == 2) ? openFracPos : openFracNeg;
 622
 623   // Spin-state extra factor 2 per incoming neutrino.
 624   if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
 625   if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
 626
 627   // Answer.
 628   return sigma;
 629
 630 }
 631
 632 //--------------------------------------------------------------------------
 633
 634 // Select identity, colour and anticolour.
 635
 636 void Sigma3ff2HchgchgfftWW::setIdColAcol() {
 637
 638   // Pick out-flavours by relative CKM weights.
 639   int id1Abs   = abs(id1);
 640   int id2Abs   = abs(id2);
 641   id4          = couplingsPtr->V2CKMpick(id1);
 642   id5          = couplingsPtr->V2CKMpick(id2);
 643
 644   // Find charge of Higgs .
 645   id3 = (( id1Abs%2 == 0 && id1 > 0) || (id1Abs%2 == 1 && id1 < 0) )
 646       ? idHLR : -idHLR;
 647   setId( id1, id2, id3, id4, id5);
 648
 649   // Colour flow topologies. Swap when antiquarks.
 650   if (id1Abs < 9 && id2Abs < 9 && id1*id2 > 0)
 651                        setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
 652   else if (id1Abs < 9 && id2Abs < 9)
 653                        setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
 654   else if (id1Abs < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
 655   else if (id2Abs < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0);
 656   else                 setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 657   if ( (id1Abs < 9 && id1 < 0) || (id1Abs > 10 && id2 < 0) )
 658     swapColAcol();
 659
 660 }
 661
 662 //--------------------------------------------------------------------------
 663
 664 // Evaluate weight for decay angles.
 665
 666 double Sigma3ff2HchgchgfftWW::weightDecay( Event& process, int iResBeg,
 667   int iResEnd) {
 668
 669   // Identity of mother of decaying reseonance(s).
 670   int idMother = process[process[iResBeg].mother1()].idAbs();
 671
 672   // For top decay hand over to standard routine.
 673   if (idMother == 6)
 674     return weightTopDecay( process, iResBeg, iResEnd);
 675
 676   // Else done.
 677   return 1.;
 678
 679 }
 680
 681 //==========================================================================
 682
 683 // Sigma2ffbar2HchgchgHchgchg class.
 684 // Cross section for f fbar -> H_(L/R)^++ H_(L/R)^-- (doubly charged Higgs).
 685
 686 //--------------------------------------------------------------------------
 687
 688 // Initialize process.
 689
 690 void Sigma2ffbar2HchgchgHchgchg::initProc() {
 691
 692   // Set process properties: H_L^++ H_L^-- or H_R^++ H_R^--.
 693   if (leftRight == 1) {
 694     idHLR      = 9900041;
 695     codeSave   = 3126;
 696     nameSave   = "f fbar -> H_L^++ H_L^--";
 697   } else {
 698     idHLR      = 9900042;
 699     codeSave   = 3146;
 700     nameSave   = "f fbar -> H_R^++ H_R^--";
 701   }
 702
 703   // Read in Yukawa matrix for couplings to a lepton pair.
 704   yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
 705   yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
 706   yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
 707   yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
 708   yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
 709   yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
 710
 711   // Electroweak parameters.
 712   mRes         = particleDataPtr->m0(23);
 713   GammaRes     = particleDataPtr->mWidth(23);
 714   m2Res        = mRes*mRes;
 715   GamMRat      = GammaRes / mRes;
 716   sin2tW       = couplingsPtr->sin2thetaW();
 717   preFac       = (1. - 2. * sin2tW) / ( 8. * sin2tW * (1. - sin2tW) );
 718
 719   // Open fraction from secondary widths.
 720   openFrac     = particleDataPtr->resOpenFrac( idHLR, -idHLR);
 721
 722 }
 723
 724 //--------------------------------------------------------------------------
 725
 726 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 727
 728 double Sigma2ffbar2HchgchgHchgchg::sigmaHat() {
 729
 730   // Electroweak couplings to gamma^*/Z^0.
 731   int    idAbs   = abs(id1);
 732   double ei      = couplingsPtr->ef(idAbs);
 733   double vi      = couplingsPtr->vf(idAbs);
 734   double ai      = couplingsPtr->af(idAbs);
 735
 736   // Part via gamma^*/Z^0 propagator.No Z^0 coupling to H_R.
 737   double resProp = 1. / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 738   double sigma   = 8. * pow2(alpEM) * ei*ei / sH2;
 739   if (leftRight == 1) sigma += 8. * pow2(alpEM)
 740     * (2. * ei * vi * preFac * (sH - m2Res) * resProp / sH
 741     + (vi * vi + ai * ai) * pow2(preFac) * resProp);
 742
 743   // Part via t-channel lepton + interference; sum over possibilities.
 744   if (idAbs == 11 || idAbs == 13 || idAbs == 15) {
 745     double yuk2Sum;
 746     if (idAbs == 11) yuk2Sum
 747       = pow2(yukawa[1][1]) + pow2(yukawa[2][1]) + pow2(yukawa[3][1]);
 748     else if (idAbs == 13) yuk2Sum
 749       = pow2(yukawa[2][1]) + pow2(yukawa[2][2]) + pow2(yukawa[3][2]);
 750     else yuk2Sum
 751       = pow2(yukawa[3][1]) + pow2(yukawa[3][2]) + pow2(yukawa[3][3]);
 752     yuk2Sum /= 4. * M_PI;
 753     sigma   += 8. * alpEM * ei * yuk2Sum / (sH * tH)
 754       + 4. * pow2(yuk2Sum) / tH2;
 755     if (leftRight == 1) sigma += 8. * alpEM * (vi + ai) * yuk2Sum
 756       * preFac * (sH - m2Res) * resProp / tH;
 757   }
 758
 759   // Common kinematical factor. Colour factor.
 760   sigma *= M_PI * (tH * uH - s3 * s4) / sH2;
 761   if (idAbs < 9) sigma /= 3.;
 762
 763   // Answer.
 764   return sigma;
 765
 766 }
 767
 768 //--------------------------------------------------------------------------
 769
 770 // Select identity, colour and anticolour.
 771
 772 void Sigma2ffbar2HchgchgHchgchg::setIdColAcol() {
 773
 774   // Outgoing flavours trivial.
 775   setId( id1, id2, idHLR, -idHLR);
 776
 777   // tHat is defined between incoming fermion and outgoing H--.
 778   if (id1 > 0) swapTU = true;
 779
 780   // No colours at all or one flow topology. Swap if first is antiquark.
 781   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
 782   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 783   if (id1 < 0) swapColAcol();
 784
 785 }
 786
 787 //--------------------------------------------------------------------------
 788
 789 // Evaluate weight for H_L/R sequential decay angles.
 790
 791 double Sigma2ffbar2HchgchgHchgchg::weightDecay( Event& process,
 792   int iResBeg, int iResEnd) {
 793
 794   // Identity of mother of decaying reseonance(s).
 795   int idMother = process[process[iResBeg].mother1()].idAbs();
 796
 797   // For top decay hand over to standard routine.
 798   if (idMother == 6)
 799     return weightTopDecay( process, iResBeg, iResEnd);
 800
 801   // Else isotropic decay.
 802   return 1.;
 803
 804 }
 805
 806 //==========================================================================
 807
 808 } // end namespace Pythia8