LCOV - code coverage report
Current view: top level - PYTHIA8/pythia8210dev/src - SigmaLeftRightSym.cc (source / functions) Hit Total Coverage
Test: coverage.info Lines: 0 356 0.0 %
Date: 2016-06-14 17:26:59 Functions: 0 27 0.0 %

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

Generated by: LCOV version 1.11