PYTHIA8/pythia8145/src/SigmaNewGaugeBosons.cxx

   1 // SigmaNewGaugeBosons.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 // leptoquark simulation classes.
   8
   9 #include "SigmaNewGaugeBosons.h"
  10
  11 namespace Pythia8 {
  12
  13
  14 //==========================================================================
  15
  16 // Sigma1ffbarZprimeWprime class.
  17 // Collects common methods for f fbar -> Z'/W' -> WW/WZ -> 4 fermions.
  18 // Copied from SigmaEW for gauge-boson-pair production.
  19
  20 //--------------------------------------------------------------------------
  21
  22 // Calculate and store internal products.
  23
  24 void Sigma1ffbarZprimeWprime::setupProd( Event& process, int i1, int i2,
  25   int i3, int i4, int i5, int i6) {
  26
  27   // Store incoming and outgoing momenta,
  28   pRot[1] = process[i1].p();
  29   pRot[2] = process[i2].p();
  30   pRot[3] = process[i3].p();
  31   pRot[4] = process[i4].p();
  32   pRot[5] = process[i5].p();
  33   pRot[6] = process[i6].p();
  34
  35   // Do random rotation to avoid accidental zeroes in HA expressions.
  36   bool smallPT = false;
  37   do {
  38     smallPT = false;
  39     double thetaNow = acos(2. * rndmPtr->flat() - 1.);
  40     double phiNow   = 2. * M_PI * rndmPtr->flat();
  41     for (int i = 1; i <= 6; ++i) {
  42       pRot[i].rot( thetaNow, phiNow);
  43       if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
  44     }
  45   } while (smallPT);
  46
  47   // Calculate internal products.
  48   for (int i = 1; i < 6; ++i) {
  49     for (int j = i + 1; j <= 6; ++j) {
  50       hA[i][j] =
  51           sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz())
  52         / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() )
  53         - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz())
  54         / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() );
  55       hC[i][j] = conj( hA[i][j] );
  56       if (i <= 2) {
  57         hA[i][j] *= complex( 0., 1.);
  58         hC[i][j] *= complex( 0., 1.);
  59       }
  60       hA[j][i] = - hA[i][j];
  61       hC[j][i] = - hC[i][j];
  62     }
  63   }
  64
  65 }
  66
  67 //--------------------------------------------------------------------------
  68
  69 // Evaluate the F function of Gunion and Kunszt.
  70
  71 complex Sigma1ffbarZprimeWprime::fGK(int j1, int j2, int j3, int j4,
  72   int j5, int j6) {
  73
  74   return 4. * hA[j1][j3] * hC[j2][j6]
  75          * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] );
  76
  77 }
  78
  79 //--------------------------------------------------------------------------
  80
  81 // Evaluate the Xi function of Gunion and Kunszt.
  82
  83 double Sigma1ffbarZprimeWprime::xiGK( double tHnow, double uHnow,
  84   double s3now, double s4now) {
  85
  86   return - 4. * s3now * s4now + tHnow * (3. * tHnow + 4. * uHnow)
  87     + tHnow * tHnow * ( tHnow * uHnow / (s3now * s4now)
  88     - 2. * (1. / s3now + 1./s4now) * (tHnow + uHnow)
  89     + 2. * (s3now / s4now + s4now / s3now) );
  90
  91 }
  92
  93 //--------------------------------------------------------------------------
  94
  95 // Evaluate the Xj function of Gunion and Kunszt.
  96
  97 double Sigma1ffbarZprimeWprime::xjGK( double tHnow, double uHnow,
  98   double s3now, double s4now) {
  99
 100   return 8. * pow2(s3now + s4now) - 8. * (s3now + s4now) * (tHnow + uHnow)
 101     - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow
 102     / (s3now * s4now) - 2. * (1. / s3now + 1. / s4now) * (tHnow + uHnow)
 103     + 2. * (s3now / s4now + s4now / s3now) );
 104
 105 }
 106
 107 //==========================================================================
 108
 109 // Sigma1ffbar2gmZZprime class.
 110 // Cross section for f fbar -> gamma*/Z0/Z'0 (f is quark or lepton).
 111
 112 //--------------------------------------------------------------------------
 113
 114 // Initialize process.
 115
 116 void Sigma1ffbar2gmZZprime::initProc() {
 117
 118   // Allow to pick only parts of full gamma*/Z0/Z'0 expression.
 119   gmZmode     = settingsPtr->mode("Zprime:gmZmode");
 120
 121   // Store Z'0 mass and width for propagator.
 122   mRes        = particleDataPtr->m0(32);
 123   GammaRes    = particleDataPtr->mWidth(32);
 124   m2Res       = mRes*mRes;
 125   GamMRat     = GammaRes / mRes;
 126   sin2tW      = couplingsPtr->sin2thetaW();
 127   cos2tW      = 1. - sin2tW;
 128   thetaWRat   = 1. / (16. * sin2tW * cos2tW);
 129
 130   // Properties of Z0 resonance also needed.
 131   mZ          = particleDataPtr->m0(23);
 132   GammaZ      = particleDataPtr->mWidth(23);
 133   m2Z         = mZ*mZ;
 134   GamMRatZ    = GammaZ / mZ;
 135
 136   // Ensure that arrays initially are empty.
 137   for (int i = 0; i < 20; ++i) afZp[i] = 0.;
 138   for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
 139
 140   // Store first-generation axial and vector couplings.
 141   afZp[1]     = settingsPtr->parm("Zprime:ad");
 142   afZp[2]     = settingsPtr->parm("Zprime:au");
 143   afZp[11]    = settingsPtr->parm("Zprime:ae");
 144   afZp[12]    = settingsPtr->parm("Zprime:anue");
 145   vfZp[1]     = settingsPtr->parm("Zprime:vd");
 146   vfZp[2]     = settingsPtr->parm("Zprime:vu");
 147   vfZp[11]    = settingsPtr->parm("Zprime:ve");
 148   vfZp[12]    = settingsPtr->parm("Zprime:vnue");
 149
 150   // Second and third generation could be carbon copy of this...
 151   if (settingsPtr->flag("Zprime:universality")) {
 152     for (int i = 3; i <= 6; ++i) {
 153       afZp[i]    = afZp[i-2];
 154       vfZp[i]    = vfZp[i-2];
 155       afZp[i+10] = afZp[i+8];
 156       vfZp[i+10] = vfZp[i+8];
 157     }
 158
 159   // ... or could have different couplings.
 160   } else {
 161     afZp[3]   = settingsPtr->parm("Zprime:as");
 162     afZp[4]   = settingsPtr->parm("Zprime:ac");
 163     afZp[5]   = settingsPtr->parm("Zprime:ab");
 164     afZp[6]   = settingsPtr->parm("Zprime:at");
 165     afZp[13]  = settingsPtr->parm("Zprime:amu");
 166     afZp[14]  = settingsPtr->parm("Zprime:anumu");
 167     afZp[15]  = settingsPtr->parm("Zprime:atau");
 168     afZp[16]  = settingsPtr->parm("Zprime:anutau");
 169     vfZp[3]   = settingsPtr->parm("Zprime:vs");
 170     vfZp[4]   = settingsPtr->parm("Zprime:vc");
 171     vfZp[5]   = settingsPtr->parm("Zprime:vb");
 172     vfZp[6]   = settingsPtr->parm("Zprime:vt");
 173     vfZp[13]  = settingsPtr->parm("Zprime:vmu");
 174     vfZp[14]  = settingsPtr->parm("Zprime:vnumu");
 175     vfZp[15]  = settingsPtr->parm("Zprime:vtau");
 176     vfZp[16]  = settingsPtr->parm("Zprime:vnutau");
 177   }
 178
 179   // Coupling for Z' -> W+ W- and decay angular admixture.
 180   coupZpWW    = settingsPtr->parm("Zprime:coup2WW");
 181   anglesZpWW  = settingsPtr->parm("Zprime:anglesWW");
 182
 183   // Set pointer to particle properties and decay table.
 184   particlePtr = particleDataPtr->particleDataEntryPtr(32);
 185
 186 }
 187
 188 //--------------------------------------------------------------------------
 189
 190 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 191
 192 void Sigma1ffbar2gmZZprime::sigmaKin() {
 193
 194   // Common coupling factors.
 195   double colQ = 3. * (1. + alpS / M_PI);
 196
 197   // Reset quantities to sum. Declare variables in loop.
 198   gamSum   = 0.;
 199   gamZSum  = 0.;
 200   ZSum     = 0.;
 201   gamZpSum = 0.;
 202   ZZpSum   = 0.;
 203   ZpSum    = 0.;
 204   int    idAbs, onMode;
 205   double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf,
 206          ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf;
 207
 208   // Loop over all open Z'0 decay channels.
 209   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
 210     onMode = particlePtr->channel(i).onMode();
 211     if (onMode != 1 && onMode != 2) continue;
 212     idAbs = abs( particlePtr->channel(i).product(0) );
 213
 214     // Contributions from three fermion generations.
 215     if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) {
 216       mf = particleDataPtr->m0(idAbs);
 217
 218       // Check that above threshold.
 219       if (mH > 2. * mf + MASSMARGIN) {
 220         mr        = pow2(mf / mH);
 221         ps        = sqrtpos(1. - 4. * mr);
 222
 223         // Couplings of gamma^*/Z^0/Z'^0  to final flavour
 224         ef        = couplingsPtr->ef(idAbs);
 225         af        = couplingsPtr->af(idAbs);
 226         vf        = couplingsPtr->vf(idAbs);
 227         apf       = afZp[idAbs];
 228         vpf       = vfZp[idAbs];
 229
 230         // Combine couplings with kinematical factors.
 231         kinFacA   = pow3(ps);
 232         kinFacV   = ps * (1. + 2. * mr);
 233         ef2       = ef * ef * kinFacV;
 234         efvf      = ef * vf * kinFacV;
 235         vaf2      = vf * vf * kinFacV + af * af * kinFacA;
 236         efvpf     = ef * vpf * kinFacV;
 237         vafvapf   = vf * vpf * kinFacV + af * apf * kinFacA;
 238         vapf2     = vpf * vpf * kinFacV + apf * apf * kinFacA;
 239
 240         // Colour factor. Additionally secondary width for top.
 241         colf      = (idAbs < 9) ? colQ : 1.;
 242         if (idAbs == 6) colf *= particleDataPtr->resOpenFrac(6, -6);
 243
 244         // Store sum of combinations.
 245         gamSum   += colf * ef2;
 246         gamZSum  += colf * efvf;
 247         ZSum     += colf * vaf2;
 248         gamZpSum += colf * efvpf;
 249         ZZpSum   += colf * vafvapf;
 250         ZpSum    += colf * vapf2;
 251       }
 252
 253     // Optional contribution from W+ W-.
 254     } else if (idAbs == 24) {
 255       mf = particleDataPtr->m0(idAbs);
 256       if (mH > 2. * mf + MASSMARGIN) {
 257         mr        = pow2(mf / mH);
 258         ps        = sqrtpos(1. - 4. * mr);
 259         ZpSum    += pow2(coupZpWW * cos2tW) * pow3(ps)
 260                   * (1. + 20. * mr + 12. * mr*mr)
 261                   * particleDataPtr->resOpenFrac(24, -24);
 262       }
 263     }
 264   }
 265
 266   // Calculate prefactors for gamma/Z0/Z'0 cross section terms.
 267   double propZ  = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
 268   double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 269   gamNorm   = 4. * M_PI * pow2(alpEM) / (3. * sH);
 270   gamZNorm  = gamNorm * 2. * thetaWRat * (sH - m2Z) * propZ;
 271   ZNorm     = gamNorm * pow2(thetaWRat) * sH * propZ;
 272   gamZpNorm = gamNorm * 2. * thetaWRat * (sH - m2Res) * propZp;
 273   ZZpNorm   = gamNorm * 2. * pow2(thetaWRat) * ((sH - m2Z) * (sH - m2Res)
 274               + sH * GamMRatZ * sH * GamMRat) * propZ * propZp;
 275   ZpNorm    = gamNorm * pow2(thetaWRat) * sH * propZp;
 276
 277   // Optionally only keep some of gamma*, Z0 and Z' terms.
 278   if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.;
 279     ZZpNorm = 0.; ZpNorm = 0.;}
 280   if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;
 281     ZZpNorm = 0.; ZpNorm = 0.;}
 282   if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
 283     gamZpNorm = 0.; ZZpNorm = 0.;}
 284   if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
 285   if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
 286   if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
 287
 288 }
 289
 290 //--------------------------------------------------------------------------
 291
 292 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 293
 294 double Sigma1ffbar2gmZZprime::sigmaHat() {
 295
 296   // Couplings to an incoming flavour.
 297   int idAbs      = abs(id1);
 298   double ei      = couplingsPtr->ef(idAbs);
 299   double ai      = couplingsPtr->af(idAbs);
 300   double vi      = couplingsPtr->vf(idAbs);
 301   double api     = afZp[idAbs];
 302   double vpi     = vfZp[idAbs];
 303   double ei2     = ei * ei;
 304   double eivi    = ei * vi;
 305   double vai2    = vi * vi + ai * ai;
 306   double eivpi   = ei * vpi;
 307   double vaivapi = vi * vpi + ai * api;;
 308   double vapi2   = vpi * vpi + api * api;
 309
 310   // Combine gamma, interference and Z0 parts.
 311   double sigma = ei2 * gamNorm * gamSum + eivi * gamZNorm * gamZSum
 312                + vai2 * ZNorm * ZSum + eivpi * gamZpNorm * gamZpSum
 313                + vaivapi * ZZpNorm * ZZpSum + vapi2 * ZpNorm * ZpSum;
 314
 315   // Colour factor. Answer.
 316   if (idAbs < 9) sigma /= 3.;
 317   return sigma;
 318
 319 }
 320
 321 //--------------------------------------------------------------------------
 322
 323 // Select identity, colour and anticolour.
 324
 325 void Sigma1ffbar2gmZZprime::setIdColAcol() {
 326
 327   // Flavours trivial.
 328   setId( id1, id2, 32);
 329
 330   // Colour flow topologies. Swap when antiquarks.
 331   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
 332   else              setColAcol( 0, 0, 0, 0, 0, 0);
 333   if (id1 < 0) swapColAcol();
 334
 335 }
 336
 337 //--------------------------------------------------------------------------
 338
 339 // Evaluate weight for gamma*/Z0/Z'0 decay angle.
 340
 341 double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg,
 342   int iResEnd) {
 343
 344   // Default values, in- and out-flavours in process.
 345   double wt    = 1.;
 346   double wtMax = 1.;
 347   int idInAbs  = process[3].idAbs();
 348   int idOutAbs = process[6].idAbs();
 349
 350   // Angular weight for outgoing fermion pair.
 351   if (iResBeg == 5 && iResEnd == 5 &&
 352     (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
 353
 354     // Couplings for in- and out-flavours.
 355     double ei  = couplingsPtr->ef(idInAbs);
 356     double vi  = couplingsPtr->vf(idInAbs);
 357     double ai  = couplingsPtr->af(idInAbs);
 358     double vpi = vfZp[idInAbs];
 359     double api = afZp[idInAbs];
 360     double ef  = couplingsPtr->ef(idOutAbs);
 361     double vf  = couplingsPtr->vf(idOutAbs);
 362     double af  = couplingsPtr->af(idOutAbs);
 363     double vpf = vfZp[idOutAbs];
 364     double apf = afZp[idOutAbs];
 365
 366     // Phase space factors. (One power of beta left out in formulae.)
 367     double mr1 = pow2(process[6].m()) / sH;
 368     double mr2 = pow2(process[7].m()) / sH;
 369     double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 370     double mrAvg = 0.5 * (mr1 + mr2) - 0.25 * pow2(mr1 - mr2);
 371
 372     // Coefficients of angular expression.
 373     double coefTran = ei*ei * gamNorm * ef*ef + ei * vi * gamZNorm * ef * vf
 374       + (vi*vi + ai*ai) * ZNorm * (vf*vf + ps*ps * af*af)
 375       + ei * vpi * gamZpNorm * ef * vpf
 376       + (vi * vpi + ai * api) * ZZpNorm * (vf * vpf + ps*ps * af * apf)
 377       + (vpi*vpi + api*api) * ZpNorm * (vpf*vpf + ps*ps * apf*apf);
 378     double coefLong = 4. * mrAvg * ( ei*ei * gamNorm * ef*ef
 379       + ei * vi * gamZNorm * ef * vf + (vi*vi + ai*ai) * ZNorm * vf*vf
 380       + ei * vpi * gamZpNorm * ef * vpf
 381       + (vi * vpi + ai * api) * ZZpNorm * vf * vpf
 382       + (vpi*vpi + api*api) * ZpNorm * vpf*vpf );
 383     double coefAsym = ps * ( ei * ai * gamZNorm * ef * af
 384       + 4. * vi * ai * ZNorm * vf * af + ei * api * gamZpNorm * ef * apf
 385       + (vi * api + vpi * ai) * ZZpNorm * (vf * apf + vpf * af)
 386       + 4. * vpi * api * ZpNorm * vpf * apf );
 387
 388     // Flip asymmetry for in-fermion + out-antifermion.
 389     if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
 390
 391     // Reconstruct decay angle and weight for it.
 392     double cosThe = (process[3].p() - process[4].p())
 393       * (process[7].p() - process[6].p()) / (sH * ps);
 394     wt    = coefTran * (1. + pow2(cosThe))
 395        + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
 396     wtMax = 2. * (coefTran + abs(coefAsym));
 397   }
 398
 399   // Angular weight for Z' -> W+ W-.
 400   else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
 401     double mr1 = pow2(process[6].m()) / sH;
 402     double mr2 = pow2(process[7].m()) / sH;
 403     double ps  = sqrtpos(pow2(1. - mr1 -mr2) - 4. * mr1 * mr2);
 404     double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2
 405       + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
 406     double cFlat = -cCos2 + 0.5 * (mr1 + mr2)
 407       * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
 408
 409     // Reconstruct decay angle and weight for it.
 410     double cosThe = (process[3].p() - process[4].p())
 411       * (process[7].p() - process[6].p()) / (sH * ps);
 412     wt    = cFlat + cCos2 * cosThe*cosThe;
 413     wtMax = cFlat + max(0., cCos2);
 414   }
 415
 416   // Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions.
 417   else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) {
 418
 419     // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
 420     // with f' fbar' from W- and f" fbar" from W+.
 421     int i1 = (process[3].id() < 0) ? 3 : 4;
 422     int i2 = 7 - i1;
 423     int i3 = (process[8].id() > 0) ? 8 : 9;
 424     int i4 = 17 - i3;
 425     int i5 = (process[10].id() > 0) ? 10 : 11;
 426     int i6 = 21 - i5;
 427     if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);}
 428
 429     // Decay distribution like in f fbar -> Z^* -> W+ W-.
 430     if (rndmPtr->flat() > anglesZpWW) {
 431
 432       // Set up four-products and internal products.
 433       setupProd( process, i1, i2, i3, i4, i5, i6);
 434
 435       // tHat and uHat of fbar f -> W- W+, and their squared masses.
 436       int iNeg     = (process[6].id() < 0) ? 6 : 7;
 437       int iPos     = 13 - iNeg;
 438       double tHres = (process[i1].p() - process[iNeg].p()).m2Calc();
 439       double uHres = (process[i1].p() - process[iPos].p()).m2Calc();
 440       double s3now = process[iNeg].m2();
 441       double s4now = process[iPos].m2();
 442
 443       // Kinematics combinations (norm(x) = |x|^2).
 444       double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
 445       double fGK253 = norm(fGK( 2, 1, 5, 6, 3, 4) - fGK( 2, 1, 3, 4, 5, 6) );
 446       double xiT    = xiGK( tHres, uHres, s3now, s4now);
 447       double xiU    = xiGK( uHres, tHres, s3now, s4now);
 448       double xjTU   = xjGK( tHres, uHres, s3now, s4now);
 449
 450       //  Couplings of incoming (anti)fermion. Combine with kinematics.
 451       int idAbs     = process[i1].idAbs();
 452       double li     = 0.5 * (vfZp[idAbs] + afZp[idAbs]);
 453       double ri     = 0.5 * (vfZp[idAbs] - afZp[idAbs]);
 454       wt            = li*li * fGK135 + ri*ri * fGK253;
 455       wtMax         = 4. * s3now * s4now * (li*li + ri*ri)
 456                     * (xiT + xiU - xjTU);
 457
 458     // Decay distribution like in f fbar -> h^0 -> W+ W-.
 459     } else {
 460       double p35  = 2. * process[i3].p() * process[i5].p();
 461       double p46  = 2. * process[i4].p() * process[i6].p();
 462       wt          = 16. * p35 * p46;
 463       wtMax       = sH2;
 464     }
 465   }
 466
 467   // Angular weight in top decay by standard routine.
 468   else if (process[process[iResBeg].mother1()].idAbs() == 6)
 469     return weightTopDecay( process, iResBeg, iResEnd);
 470
 471
 472   // Done.
 473   return (wt / wtMax);
 474
 475 }
 476
 477 //==========================================================================
 478
 479 // Sigma1ffbar2Wprime class.
 480 // Cross section for f fbar' -> W'+- (f is quark or lepton).
 481
 482 //--------------------------------------------------------------------------
 483
 484 // Initialize process.
 485
 486 void Sigma1ffbar2Wprime::initProc() {
 487
 488   // Store W+- mass and width for propagator.
 489   mRes     = particleDataPtr->m0(34);
 490   GammaRes = particleDataPtr->mWidth(34);
 491   m2Res    = mRes*mRes;
 492   GamMRat  = GammaRes / mRes;
 493   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
 494
 495   // Axial and vector couplings of fermions.
 496   aqWp      = settingsPtr->parm("Wprime:aq");
 497   vqWp      = settingsPtr->parm("Wprime:vq");
 498   alWp      = settingsPtr->parm("Wprime:al");
 499   vlWp      = settingsPtr->parm("Wprime:vl");
 500
 501   // Coupling for W' -> W Z and decay angular admixture.
 502   coupWpWZ    = settingsPtr->parm("Wprime:coup2WZ");
 503   anglesWpWZ  = settingsPtr->parm("Wprime:anglesWZ");
 504
 505   // Set pointer to particle properties and decay table.
 506   particlePtr = particleDataPtr->particleDataEntryPtr(34);
 507
 508 }
 509
 510 //--------------------------------------------------------------------------
 511
 512 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 513
 514 void Sigma1ffbar2Wprime::sigmaKin() {
 515
 516   // Set up Breit-Wigner. Cross section for W+ and W- separately.
 517   double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 518   double preFac = alpEM * thetaWRat * mH;
 519   sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(34, mH);
 520   sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-34, mH);
 521
 522 }
 523
 524 //--------------------------------------------------------------------------
 525
 526 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 527
 528 double Sigma1ffbar2Wprime::sigmaHat() {
 529
 530   // Secondary width for W+ or W-. CKM and colour factors.
 531   int idUp = (abs(id1)%2 == 0) ? id1 : id2;
 532   double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
 533   if (abs(id1) < 7) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
 534
 535   // Couplings.
 536   if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp);
 537   else              sigma *= 0.5 * (alWp * alWp + vlWp * vlWp);
 538
 539   // Answer.
 540   return sigma;
 541
 542 }
 543
 544 //--------------------------------------------------------------------------
 545
 546 // Select identity, colour and anticolour.
 547
 548 void Sigma1ffbar2Wprime::setIdColAcol() {
 549
 550   // Sign of outgoing W.
 551   int sign          = 1 - 2 * (abs(id1)%2);
 552   if (id1 < 0) sign = -sign;
 553   setId( id1, id2, 34 * sign);
 554
 555   // Colour flow topologies. Swap when antiquarks.
 556   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
 557   else              setColAcol( 0, 0, 0, 0, 0, 0);
 558   if (id1 < 0) swapColAcol();
 559
 560 }
 561
 562 //--------------------------------------------------------------------------
 563
 564 // Evaluate weight for W decay angle.
 565
 566 double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg,
 567   int iResEnd) {
 568
 569   // Default values, in- and out-flavours in process.
 570   double wt    = 1.;
 571   double wtMax = 1.;
 572   int idInAbs  = process[3].idAbs();
 573   int idOutAbs = process[6].idAbs();
 574
 575   // Angular weight for outgoing fermion pair.
 576   if (iResBeg == 5 && iResEnd == 5 &&
 577     (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {
 578
 579     // Couplings for in- and out-flavours.
 580     double ai  = (idInAbs < 9) ? aqWp : alWp;
 581     double vi  = (idInAbs < 9) ? vqWp : vlWp;
 582     double af  = (idOutAbs < 9) ? aqWp : alWp;
 583     double vf  = (idOutAbs < 9) ? vqWp : vlWp;
 584
 585     // Asymmetry expression.
 586     double coefAsym = 8. * vi * ai * vf * af
 587       / ((vi*vi + ai*ai) * (vf*vf + af*af));
 588
 589     // Flip asymmetry for in-fermion + out-antifermion.
 590     if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;
 591
 592     // Phase space factors.
 593     double mr1 = pow2(process[6].m()) / sH;
 594     double mr2 = pow2(process[7].m()) / sH;
 595     double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 596
 597     // Reconstruct decay angle and weight for it.
 598     double cosThe = (process[3].p() - process[4].p())
 599       * (process[7].p() - process[6].p()) / (sH * ps);
 600     wt    = 1. + coefAsym * cosThe + cosThe * cosThe;
 601     wtMax = 2. + abs(coefAsym);
 602   }
 603
 604   // Angular weight for W' -> W Z.
 605   else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
 606     double mr1 = pow2(process[6].m()) / sH;
 607     double mr2 = pow2(process[7].m()) / sH;
 608     double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
 609     double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2
 610       + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
 611     double cFlat = -cCos2 + 0.5 * (mr1 + mr2)
 612       * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));
 613
 614     // Reconstruct decay angle and weight for it.
 615     double cosThe = (process[3].p() - process[4].p())
 616       * (process[7].p() - process[6].p()) / (sH * ps);
 617     wt    = cFlat + cCos2 * cosThe*cosThe;
 618     wtMax = cFlat + max(0., cCos2);
 619   }
 620
 621   // Angular weight for f + fbar -> W' -> W + Z -> 4 fermions.
 622   else if (iResBeg == 6 && iResEnd == 7
 623     && (idOutAbs == 24 || idOutAbs == 23)) {
 624
 625     // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
 626     // with f' fbar' from W and f" fbar" from Z.
 627     int i1 = (process[3].id() < 0) ? 3 : 4;
 628     int i2 = 7 - i1;
 629     int i3 = (process[8].id() > 0) ? 8 : 9;
 630     int i4 = 17 - i3;
 631     int i5 = (process[10].id() > 0) ? 10 : 11;
 632     int i6 = 21 - i5;
 633     if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);}
 634
 635     // Decay distribution like in f fbar -> Z^* -> W+ W-.
 636     if (rndmPtr->flat() > anglesWpWZ) {
 637
 638       // Set up four-products and internal products.
 639       setupProd( process, i1, i2, i3, i4, i5, i6);
 640
 641       // tHat and uHat of fbar f -> W Z, and their squared masses.
 642       int iW       = (process[6].id() == 23) ? 7 : 6;
 643       int iZ       = 13 - iW;
 644       double tHres = (process[i1].p() - process[iW].p()).m2Calc();
 645       double uHres = (process[i1].p() - process[iZ].p()).m2Calc();
 646       double s3now = process[iW].m2();
 647       double s4now = process[iZ].m2();
 648
 649       // Kinematics combinations (norm(x) = |x|^2).
 650       double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
 651       double fGK136 = norm(fGK( 1, 2, 3, 4, 6, 5) - fGK( 1, 2, 6, 5, 3, 4) );
 652       double xiT    = xiGK( tHres, uHres, s3now, s4now);
 653       double xiU    = xiGK( uHres, tHres, s3now, s4now);
 654       double xjTU   = xjGK( tHres, uHres, s3now, s4now);
 655
 656       //  Couplings of outgoing fermion from Z. Combine with kinematics.
 657       int idAbs     = process[i5].idAbs();
 658       double lfZ    = couplingsPtr->lf(idAbs);
 659       double rfZ    = couplingsPtr->rf(idAbs);
 660       wt            = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136;
 661       wtMax         = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ)
 662                     * (xiT + xiU - xjTU);
 663
 664     // Decay distribution like in f fbar -> H^+- -> W+- Z0.
 665     } else {
 666       double p35  = 2. * process[i3].p() * process[i5].p();
 667       double p46  = 2. * process[i4].p() * process[i6].p();
 668       wt          = 16. * p35 * p46;
 669       wtMax       = sH2;
 670     }
 671   }
 672
 673   // Angular weight in top decay by standard routine.
 674   else if (process[process[iResBeg].mother1()].idAbs() == 6)
 675     return weightTopDecay( process, iResBeg, iResEnd);
 676
 677   // Done.
 678   return (wt / wtMax);
 679
 680 }
 681
 682
 683 //==========================================================================
 684
 685 // Sigma1ffbar2Rhorizontal class.
 686 // Cross section for f fbar' -> R^0 (f is a quark or lepton).
 687
 688 //--------------------------------------------------------------------------
 689
 690 // Initialize process.
 691
 692 void Sigma1ffbar2Rhorizontal::initProc() {
 693
 694   // Store R^0 mass and width for propagator.
 695   mRes     = particleDataPtr->m0(41);
 696   GammaRes = particleDataPtr->mWidth(41);
 697   m2Res    = mRes*mRes;
 698   GamMRat  = GammaRes / mRes;
 699   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
 700
 701   // Set pointer to particle properties and decay table.
 702   particlePtr = particleDataPtr->particleDataEntryPtr(41);
 703
 704 }
 705
 706 //--------------------------------------------------------------------------
 707
 708 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
 709
 710 void Sigma1ffbar2Rhorizontal::sigmaKin() {
 711
 712   // Set up Breit-Wigner. Cross section for W+ and W- separately.
 713   double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
 714   double preFac = alpEM * thetaWRat * mH;
 715   sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(41, mH);
 716   sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-41, mH);
 717
 718 }
 719
 720 //--------------------------------------------------------------------------
 721
 722 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
 723
 724 double Sigma1ffbar2Rhorizontal::sigmaHat() {
 725
 726   // Check for allowed flavour combinations, one generation apart.
 727   if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.;
 728
 729   // Find whether R0 or R0bar. Colour factors.
 730   double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg;
 731   if (abs(id1) < 7) sigma /= 3.;
 732
 733   // Answer.
 734   return sigma;
 735
 736 }
 737
 738 //--------------------------------------------------------------------------
 739
 740 // Select identity, colour and anticolour.
 741
 742 void Sigma1ffbar2Rhorizontal::setIdColAcol() {
 743
 744   // Outgoing R0 or R0bar.
 745   id3 = (id1 +id2 > 0) ? 41 : -41;
 746   setId( id1, id2, id3);
 747
 748   // Colour flow topologies. Swap when antiquarks.
 749   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
 750   else              setColAcol( 0, 0, 0, 0, 0, 0);
 751   if (id1 < 0) swapColAcol();
 752
 753 }
 754
 755 //==========================================================================
 756
 757 } // end namespace Pythia8