// SigmaNewGaugeBosons.cc is a part of the PYTHIA event generator.
// Copyright (C) 2013 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// Function definitions (not found in the header) for the 
// leptoquark simulation classes. 

#include "SigmaNewGaugeBosons.h"

namespace Pythia8 {


//==========================================================================

// Sigma1ffbarZprimeWprime class.
// Collects common methods for f fbar -> Z'/W' -> WW/WZ -> 4 fermions.
// Copied from SigmaEW for gauge-boson-pair production.

//--------------------------------------------------------------------------

// Calculate and store internal products.

void Sigma1ffbarZprimeWprime::setupProd( Event& process, int i1, int i2, 
  int i3, int i4, int i5, int i6) {

  // Store incoming and outgoing momenta,
  pRot[1] = process[i1].p();
  pRot[2] = process[i2].p();
  pRot[3] = process[i3].p();
  pRot[4] = process[i4].p();
  pRot[5] = process[i5].p();
  pRot[6] = process[i6].p();

  // Do random rotation to avoid accidental zeroes in HA expressions. 
  bool smallPT = false;
  do {
    smallPT = false;
    double thetaNow = acos(2. * rndmPtr->flat() - 1.);
    double phiNow   = 2. * M_PI * rndmPtr->flat();
    for (int i = 1; i <= 6; ++i) { 
      pRot[i].rot( thetaNow, phiNow);
      if (pRot[i].pT2() < 1e-4 * pRot[i].pAbs2()) smallPT = true;
    }
  } while (smallPT); 

  // Calculate internal products.
  for (int i = 1; i < 6; ++i) {
    for (int j = i + 1; j <= 6; ++j) { 
      hA[i][j] = 
          sqrt( (pRot[i].e() - pRot[i].pz()) * (pRot[j].e() + pRot[j].pz()) 
        / pRot[i].pT2() ) * complex( pRot[i].px(), pRot[i].py() ) 
        - sqrt( (pRot[i].e() + pRot[i].pz()) * (pRot[j].e() - pRot[j].pz()) 
        / pRot[j].pT2() ) * complex( pRot[j].px(), pRot[j].py() ); 
      hC[i][j] = conj( hA[i][j] );
      if (i <= 2) {
        hA[i][j] *= complex( 0., 1.);
        hC[i][j] *= complex( 0., 1.);
      }
      hA[j][i] = - hA[i][j]; 
      hC[j][i] = - hC[i][j]; 
    }
  }

}

//--------------------------------------------------------------------------

// Evaluate the F function of Gunion and Kunszt.

complex Sigma1ffbarZprimeWprime::fGK(int j1, int j2, int j3, int j4, 
  int j5, int j6) {
 
  return 4. * hA[j1][j3] * hC[j2][j6] 
         * ( hA[j1][j5] * hC[j1][j4] + hA[j3][j5] * hC[j3][j4] ); 

}

//--------------------------------------------------------------------------

// Evaluate the Xi function of Gunion and Kunszt.

double Sigma1ffbarZprimeWprime::xiGK( double tHnow, double uHnow, 
  double s3now, double s4now) {

  return - 4. * s3now * s4now + tHnow * (3. * tHnow + 4. * uHnow) 
    + tHnow * tHnow * ( tHnow * uHnow / (s3now * s4now)
    - 2. * (1. / s3now + 1./s4now) * (tHnow + uHnow)  
    + 2. * (s3now / s4now + s4now / s3now) );

}

//--------------------------------------------------------------------------

// Evaluate the Xj function of Gunion and Kunszt.

double Sigma1ffbarZprimeWprime::xjGK( double tHnow, double uHnow,
  double s3now, double s4now) {

  return 8. * pow2(s3now + s4now) - 8. * (s3now + s4now) * (tHnow + uHnow)
    - 6. * tHnow * uHnow - 2. * tHnow * uHnow * ( tHnow * uHnow 
    / (s3now * s4now) - 2. * (1. / s3now + 1. / s4now) * (tHnow + uHnow) 
    + 2. * (s3now / s4now + s4now / s3now) );

}

//==========================================================================

// Sigma1ffbar2gmZZprime class.
// Cross section for f fbar -> gamma*/Z0/Z'0 (f is quark or lepton). 

//--------------------------------------------------------------------------

// Initialize process. 
  
void Sigma1ffbar2gmZZprime::initProc() {

  // Allow to pick only parts of full gamma*/Z0/Z'0 expression.
  gmZmode     = settingsPtr->mode("Zprime:gmZmode");

  // Store Z'0 mass and width for propagator. 
  mRes        = particleDataPtr->m0(32);
  GammaRes    = particleDataPtr->mWidth(32);
  m2Res       = mRes*mRes;
  GamMRat     = GammaRes / mRes;
  sin2tW      = couplingsPtr->sin2thetaW();
  cos2tW      = 1. - sin2tW;
  thetaWRat   = 1. / (16. * sin2tW * cos2tW);

  // Properties of Z0 resonance also needed.
  mZ          = particleDataPtr->m0(23);
  GammaZ      = particleDataPtr->mWidth(23);
  m2Z         = mZ*mZ;
  GamMRatZ    = GammaZ / mZ;

  // Ensure that arrays initially are empty.
  for (int i = 0; i < 20; ++i) afZp[i] = 0.;
  for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
  
  // Store first-generation axial and vector couplings.
  afZp[1]     = settingsPtr->parm("Zprime:ad");
  afZp[2]     = settingsPtr->parm("Zprime:au");
  afZp[11]    = settingsPtr->parm("Zprime:ae");
  afZp[12]    = settingsPtr->parm("Zprime:anue");
  vfZp[1]     = settingsPtr->parm("Zprime:vd");
  vfZp[2]     = settingsPtr->parm("Zprime:vu");
  vfZp[11]    = settingsPtr->parm("Zprime:ve");
  vfZp[12]    = settingsPtr->parm("Zprime:vnue");

  // Second and third generation could be carbon copy of this...
  if (settingsPtr->flag("Zprime:universality")) {
    for (int i = 3; i <= 6; ++i) {
      afZp[i]    = afZp[i-2];
      vfZp[i]    = vfZp[i-2];
      afZp[i+10] = afZp[i+8];
      vfZp[i+10] = vfZp[i+8];
    }

  // ... or could have different couplings.   
  } else {
    afZp[3]   = settingsPtr->parm("Zprime:as");
    afZp[4]   = settingsPtr->parm("Zprime:ac");
    afZp[5]   = settingsPtr->parm("Zprime:ab");
    afZp[6]   = settingsPtr->parm("Zprime:at");
    afZp[13]  = settingsPtr->parm("Zprime:amu");
    afZp[14]  = settingsPtr->parm("Zprime:anumu");
    afZp[15]  = settingsPtr->parm("Zprime:atau");
    afZp[16]  = settingsPtr->parm("Zprime:anutau");
    vfZp[3]   = settingsPtr->parm("Zprime:vs");
    vfZp[4]   = settingsPtr->parm("Zprime:vc");
    vfZp[5]   = settingsPtr->parm("Zprime:vb");
    vfZp[6]   = settingsPtr->parm("Zprime:vt");
    vfZp[13]  = settingsPtr->parm("Zprime:vmu");
    vfZp[14]  = settingsPtr->parm("Zprime:vnumu");
    vfZp[15]  = settingsPtr->parm("Zprime:vtau");
    vfZp[16]  = settingsPtr->parm("Zprime:vnutau");
  }

  // Coupling for Z' -> W+ W- and decay angular admixture.
  coupZpWW    = settingsPtr->parm("Zprime:coup2WW");
  anglesZpWW  = settingsPtr->parm("Zprime:anglesWW");

  // Set pointer to particle properties and decay table.
  particlePtr = particleDataPtr->particleDataEntryPtr(32);

} 

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat), part independent of incoming flavour. 

void Sigma1ffbar2gmZZprime::sigmaKin() { 

  // Common coupling factors.
  double colQ = 3. * (1. + alpS / M_PI);

  // Reset quantities to sum. Declare variables in loop.
  gamSum   = 0.;
  gamZSum  = 0.;
  ZSum     = 0.;
  gamZpSum = 0.;
  ZZpSum   = 0.;
  ZpSum    = 0.;
  int    idAbs, onMode;
  double mf, mr, ps, kinFacA, kinFacV, ef, af, vf, apf, vpf,
         ef2, efvf, vaf2, efvpf, vafvapf, vapf2, colf;

  // Loop over all open Z'0 decay channels. 
  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
    onMode = particlePtr->channel(i).onMode();
    if (onMode != 1 && onMode != 2) continue;
    idAbs = abs( particlePtr->channel(i).product(0) );

    // Contributions from three fermion generations.
    if ( (idAbs > 0 && idAbs < 7) || ( idAbs > 10 && idAbs < 17) ) {
      mf = particleDataPtr->m0(idAbs);

      // Check that above threshold. 
      if (mH > 2. * mf + MASSMARGIN) {
        mr        = pow2(mf / mH);
        ps        = sqrtpos(1. - 4. * mr); 

        // Couplings of gamma^*/Z^0/Z'^0  to final flavour 
        ef        = couplingsPtr->ef(idAbs);
        af        = couplingsPtr->af(idAbs);
        vf        = couplingsPtr->vf(idAbs);
        apf       = afZp[idAbs];
        vpf       = vfZp[idAbs];

        // Combine couplings with kinematical factors.
        kinFacA   = pow3(ps);
        kinFacV   = ps * (1. + 2. * mr);
        ef2       = ef * ef * kinFacV;
        efvf      = ef * vf * kinFacV;
        vaf2      = vf * vf * kinFacV + af * af * kinFacA; 
        efvpf     = ef * vpf * kinFacV;
        vafvapf   = vf * vpf * kinFacV + af * apf * kinFacA;
        vapf2     = vpf * vpf * kinFacV + apf * apf * kinFacA; 

        // Colour factor. Additionally secondary width for top.
        colf      = (idAbs < 9) ? colQ : 1.;
        if (idAbs == 6) colf *= particleDataPtr->resOpenFrac(6, -6);

        // Store sum of combinations.
        gamSum   += colf * ef2;
        gamZSum  += colf * efvf;
        ZSum     += colf * vaf2;
        gamZpSum += colf * efvpf;
        ZZpSum   += colf * vafvapf;
        ZpSum    += colf * vapf2;
      }

    // Optional contribution from W+ W-.
    } else if (idAbs == 24) {
      mf = particleDataPtr->m0(idAbs);
      if (mH > 2. * mf + MASSMARGIN) {
        mr        = pow2(mf / mH);
        ps        = sqrtpos(1. - 4. * mr); 
        ZpSum    += pow2(coupZpWW * cos2tW) * pow3(ps)
          	  * (1. + 20. * mr + 12. * mr*mr)
                  * particleDataPtr->resOpenFrac(24, -24);
      }
    }
  }

  // Calculate prefactors for gamma/Z0/Z'0 cross section terms.
  double propZ  = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
  double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
  gamNorm   = 4. * M_PI * pow2(alpEM) / (3. * sH); 
  gamZNorm  = gamNorm * 2. * thetaWRat * (sH - m2Z) * propZ;
  ZNorm     = gamNorm * pow2(thetaWRat) * sH * propZ;
  gamZpNorm = gamNorm * 2. * thetaWRat * (sH - m2Res) * propZp;
  ZZpNorm   = gamNorm * 2. * pow2(thetaWRat) * ((sH - m2Z) * (sH - m2Res)
              + sH * GamMRatZ * sH * GamMRat) * propZ * propZp;
  ZpNorm    = gamNorm * pow2(thetaWRat) * sH * propZp;

  // Optionally only keep some of gamma*, Z0 and Z' terms.
  if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.; 
    ZZpNorm = 0.; ZpNorm = 0.;}
  if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.; 
    ZZpNorm = 0.; ZpNorm = 0.;}
  if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
    gamZpNorm = 0.; ZZpNorm = 0.;}
  if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
  if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
  if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}

}

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat), including incoming flavour dependence. 

double Sigma1ffbar2gmZZprime::sigmaHat() { 

  // Couplings to an incoming flavour.
  int idAbs      = abs(id1); 
  double ei      = couplingsPtr->ef(idAbs);
  double ai      = couplingsPtr->af(idAbs);
  double vi      = couplingsPtr->vf(idAbs);
  double api     = afZp[idAbs];
  double vpi     = vfZp[idAbs];
  double ei2     = ei * ei;
  double eivi    = ei * vi;
  double vai2    = vi * vi + ai * ai; 
  double eivpi   = ei * vpi;
  double vaivapi = vi * vpi + ai * api;; 
  double vapi2   = vpi * vpi + api * api;

  // Combine gamma, interference and Z0 parts.
  double sigma = ei2 * gamNorm * gamSum + eivi * gamZNorm * gamZSum
               + vai2 * ZNorm * ZSum + eivpi * gamZpNorm * gamZpSum
               + vaivapi * ZZpNorm * ZZpSum + vapi2 * ZpNorm * ZpSum;

  // Colour factor. Answer.
  if (idAbs < 9) sigma /= 3.;
  return sigma;    

}

//--------------------------------------------------------------------------

// Select identity, colour and anticolour.

void Sigma1ffbar2gmZZprime::setIdColAcol() {

  // Flavours trivial.
  setId( id1, id2, 32);

  // Colour flow topologies. Swap when antiquarks.
  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
  else              setColAcol( 0, 0, 0, 0, 0, 0);
  if (id1 < 0) swapColAcol();

}

//--------------------------------------------------------------------------

// Evaluate weight for gamma*/Z0/Z'0 decay angle.
  
double Sigma1ffbar2gmZZprime::weightDecay( Event& process, int iResBeg, 
  int iResEnd) {

  // Default values, in- and out-flavours in process.
  double wt    = 1.;
  double wtMax = 1.;
  int idInAbs  = process[3].idAbs();
  int idOutAbs = process[6].idAbs();

  // Angular weight for outgoing fermion pair.
  if (iResBeg == 5 && iResEnd == 5 &&
    (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {

    // Couplings for in- and out-flavours.
    double ei  = couplingsPtr->ef(idInAbs);
    double vi  = couplingsPtr->vf(idInAbs);
    double ai  = couplingsPtr->af(idInAbs);
    double vpi = vfZp[idInAbs];
    double api = afZp[idInAbs];
    double ef  = couplingsPtr->ef(idOutAbs);
    double vf  = couplingsPtr->vf(idOutAbs);
    double af  = couplingsPtr->af(idOutAbs);
    double vpf = vfZp[idOutAbs];
    double apf = afZp[idOutAbs];

    // Phase space factors. (One power of beta left out in formulae.)
    double mr1 = pow2(process[6].m()) / sH;
    double mr2 = pow2(process[7].m()) / sH;
    double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 
    double mrAvg = 0.5 * (mr1 + mr2) - 0.25 * pow2(mr1 - mr2);

    // Coefficients of angular expression.
    double coefTran = ei*ei * gamNorm * ef*ef + ei * vi * gamZNorm * ef * vf
      + (vi*vi + ai*ai) * ZNorm * (vf*vf + ps*ps * af*af)
      + ei * vpi * gamZpNorm * ef * vpf
      + (vi * vpi + ai * api) * ZZpNorm * (vf * vpf + ps*ps * af * apf) 
      + (vpi*vpi + api*api) * ZpNorm * (vpf*vpf + ps*ps * apf*apf);
    double coefLong = 4. * mrAvg * ( ei*ei * gamNorm * ef*ef 
      + ei * vi * gamZNorm * ef * vf + (vi*vi + ai*ai) * ZNorm * vf*vf
      + ei * vpi * gamZpNorm * ef * vpf 
      + (vi * vpi + ai * api) * ZZpNorm * vf * vpf
      + (vpi*vpi + api*api) * ZpNorm * vpf*vpf );
    double coefAsym = ps * ( ei * ai * gamZNorm * ef * af 
      + 4. * vi * ai * ZNorm * vf * af + ei * api * gamZpNorm * ef * apf 
      + (vi * api + vpi * ai) * ZZpNorm * (vf * apf + vpf * af)
      + 4. * vpi * api * ZpNorm * vpf * apf );

    // Flip asymmetry for in-fermion + out-antifermion.
    if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;

    // Reconstruct decay angle and weight for it.
    double cosThe = (process[3].p() - process[4].p()) 
      * (process[7].p() - process[6].p()) / (sH * ps);
    wt    = coefTran * (1. + pow2(cosThe)) 
       + coefLong * (1. - pow2(cosThe)) + 2. * coefAsym * cosThe;
    wtMax = 2. * (coefTran + abs(coefAsym));
  }

  // Angular weight for Z' -> W+ W-.
  else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
    double mr1 = pow2(process[6].m()) / sH;
    double mr2 = pow2(process[7].m()) / sH;
    double ps  = sqrtpos(pow2(1. - mr1 -mr2) - 4. * mr1 * mr2);
    double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2 
      + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
    double cFlat = -cCos2 + 0.5 * (mr1 + mr2) 
      * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));

    // Reconstruct decay angle and weight for it.
    double cosThe = (process[3].p() - process[4].p()) 
      * (process[7].p() - process[6].p()) / (sH * ps);
    wt    = cFlat + cCos2 * cosThe*cosThe; 
    wtMax = cFlat + max(0., cCos2);
  }

  // Angular weight for f + fbar -> Z' -> W+ + W- -> 4 fermions.
  else if (iResBeg == 6 && iResEnd == 7 && idOutAbs == 24) {
 
    // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
    // with f' fbar' from W- and f" fbar" from W+.
    int i1 = (process[3].id() < 0) ? 3 : 4;
    int i2 = 7 - i1; 
    int i3 = (process[8].id() > 0) ? 8 : 9;
    int i4 = 17 - i3;
    int i5 = (process[10].id() > 0) ? 10 : 11;
    int i6 = 21 - i5;
    if (process[6].id() > 0) {swap(i3, i5); swap(i4, i6);}
  
    // Decay distribution like in f fbar -> Z^* -> W+ W-. 
    if (rndmPtr->flat() > anglesZpWW) {

      // Set up four-products and internal products.
      setupProd( process, i1, i2, i3, i4, i5, i6); 

      // tHat and uHat of fbar f -> W- W+, and their squared masses.
      int iNeg     = (process[6].id() < 0) ? 6 : 7;
      int iPos     = 13 - iNeg;
      double tHres = (process[i1].p() - process[iNeg].p()).m2Calc();
      double uHres = (process[i1].p() - process[iPos].p()).m2Calc(); 
      double s3now = process[iNeg].m2();
      double s4now = process[iPos].m2();

      // Kinematics combinations (norm(x) = |x|^2).
      double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
      double fGK253 = norm(fGK( 2, 1, 5, 6, 3, 4) - fGK( 2, 1, 3, 4, 5, 6) );
      double xiT    = xiGK( tHres, uHres, s3now, s4now);
      double xiU    = xiGK( uHres, tHres, s3now, s4now);
      double xjTU   = xjGK( tHres, uHres, s3now, s4now);

      //  Couplings of incoming (anti)fermion. Combine with kinematics.
      int idAbs     = process[i1].idAbs();
      double li     = 0.5 * (vfZp[idAbs] + afZp[idAbs]); 
      double ri     = 0.5 * (vfZp[idAbs] - afZp[idAbs]);  
      wt            = li*li * fGK135 + ri*ri * fGK253;
      wtMax         = 4. * s3now * s4now * (li*li + ri*ri) 
                    * (xiT + xiU - xjTU);
   
    // Decay distribution like in f fbar -> h^0 -> W+ W-. 
    } else {
      double p35  = 2. * process[i3].p() * process[i5].p(); 
      double p46  = 2. * process[i4].p() * process[i6].p(); 
      wt          = 16. * p35 * p46;
      wtMax       = sH2;
    }
  }

  // Angular weight in top decay by standard routine.
  else if (process[process[iResBeg].mother1()].idAbs() == 6)
    return weightTopDecay( process, iResBeg, iResEnd);


  // Done.
  return (wt / wtMax);

}

//==========================================================================

// Sigma1ffbar2Wprime class.
// Cross section for f fbar' -> W'+- (f is quark or lepton). 

//--------------------------------------------------------------------------

// Initialize process. 
  
void Sigma1ffbar2Wprime::initProc() {

  // Store W+- mass and width for propagator. 
  mRes     = particleDataPtr->m0(34);
  GammaRes = particleDataPtr->mWidth(34);
  m2Res    = mRes*mRes;
  GamMRat  = GammaRes / mRes;
  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());

  // Axial and vector couplings of fermions.
  aqWp      = settingsPtr->parm("Wprime:aq");
  vqWp      = settingsPtr->parm("Wprime:vq");
  alWp      = settingsPtr->parm("Wprime:al");
  vlWp      = settingsPtr->parm("Wprime:vl");

  // Coupling for W' -> W Z and decay angular admixture.
  coupWpWZ    = settingsPtr->parm("Wprime:coup2WZ");
  anglesWpWZ  = settingsPtr->parm("Wprime:anglesWZ");

  // Set pointer to particle properties and decay table.
  particlePtr = particleDataPtr->particleDataEntryPtr(34);

} 

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat), part independent of incoming flavour. 

void Sigma1ffbar2Wprime::sigmaKin() { 

  // Set up Breit-Wigner. Cross section for W+ and W- separately.
  double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
  double preFac = alpEM * thetaWRat * mH;
  sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(34, mH);    
  sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-34, mH);    

}

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat), including incoming flavour dependence. 

double Sigma1ffbar2Wprime::sigmaHat() {

  // Secondary width for W+ or W-. CKM and colour factors.
  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
  double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
  if (abs(id1) < 7) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;

  // Couplings.
  if (abs(id1) < 7) sigma *= 0.5 * (aqWp * aqWp + vqWp * vqWp);
  else              sigma *= 0.5 * (alWp * alWp + vlWp * vlWp);        

  // Answer.
  return sigma;    

}

//--------------------------------------------------------------------------

// Select identity, colour and anticolour.

void Sigma1ffbar2Wprime::setIdColAcol() {

  // Sign of outgoing W.
  int sign          = 1 - 2 * (abs(id1)%2);
  if (id1 < 0) sign = -sign;
  setId( id1, id2, 34 * sign);

  // Colour flow topologies. Swap when antiquarks.
  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
  else              setColAcol( 0, 0, 0, 0, 0, 0);
  if (id1 < 0) swapColAcol();

}

//--------------------------------------------------------------------------

// Evaluate weight for W decay angle.
  
double Sigma1ffbar2Wprime::weightDecay( Event& process, int iResBeg, 
  int iResEnd) {

  // Default values, in- and out-flavours in process.
  double wt    = 1.;
  double wtMax = 1.;
  int idInAbs  = process[3].idAbs();
  int idOutAbs = process[6].idAbs();

  // Angular weight for outgoing fermion pair.
  if (iResBeg == 5 && iResEnd == 5 &&
    (idOutAbs < 7 || ( idOutAbs > 10 && idOutAbs < 17)) ) {

    // Couplings for in- and out-flavours.
    double ai  = (idInAbs < 9) ? aqWp : alWp;
    double vi  = (idInAbs < 9) ? vqWp : vlWp;
    double af  = (idOutAbs < 9) ? aqWp : alWp;
    double vf  = (idOutAbs < 9) ? vqWp : vlWp;

    // Asymmetry expression.
    double coefAsym = 8. * vi * ai * vf * af 
      / ((vi*vi + ai*ai) * (vf*vf + af*af));

    // Flip asymmetry for in-fermion + out-antifermion.
    if (process[3].id() * process[6].id() < 0) coefAsym = -coefAsym;

    // Phase space factors.
    double mr1 = pow2(process[6].m()) / sH;
    double mr2 = pow2(process[7].m()) / sH;
    double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2); 

    // Reconstruct decay angle and weight for it.
    double cosThe = (process[3].p() - process[4].p()) 
      * (process[7].p() - process[6].p()) / (sH * ps);
    wt    = 1. + coefAsym * cosThe + cosThe * cosThe;
    wtMax = 2. + abs(coefAsym);
  }

  // Angular weight for W' -> W Z.
  else if (iResBeg == 5 && iResEnd == 5 && idOutAbs == 24) {
    double mr1 = pow2(process[6].m()) / sH;
    double mr2 = pow2(process[7].m()) / sH;
    double ps  = sqrtpos(pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
    double cCos2 = - (1./16.) * ps*ps * (1. - 2. * mr1 - 2. * mr2 
      + mr1*mr1 + mr2*mr2 + 10. * mr1 * mr2);
    double cFlat = -cCos2 + 0.5 * (mr1 + mr2) 
      * (1. - 2. * mr1 - 2. * mr2 + pow2(mr1 - mr2));

    // Reconstruct decay angle and weight for it.
    double cosThe = (process[3].p() - process[4].p()) 
      * (process[7].p() - process[6].p()) / (sH * ps);
    wt    = cFlat + cCos2 * cosThe*cosThe; 
    wtMax = cFlat + max(0., cCos2);
  }

  // Angular weight for f + fbar -> W' -> W + Z -> 4 fermions.
  else if (iResBeg == 6 && iResEnd == 7 
    && (idOutAbs == 24 || idOutAbs == 23)) {
 
    // Order so that fbar(1) f(2) -> f'(3) fbar'(4) f"(5) fbar"(6).
    // with f' fbar' from W and f" fbar" from Z.
    int i1 = (process[3].id() < 0) ? 3 : 4;
    int i2 = 7 - i1; 
    int i3 = (process[8].id() > 0) ? 8 : 9;
    int i4 = 17 - i3;
    int i5 = (process[10].id() > 0) ? 10 : 11;
    int i6 = 21 - i5;
    if (process[6].id() == 23) {swap(i3, i5); swap(i4, i6);}
  
    // Decay distribution like in f fbar -> Z^* -> W+ W-. 
    if (rndmPtr->flat() > anglesWpWZ) {

      // Set up four-products and internal products.
      setupProd( process, i1, i2, i3, i4, i5, i6); 

      // tHat and uHat of fbar f -> W Z, and their squared masses.
      int iW       = (process[6].id() == 23) ? 7 : 6;
      int iZ       = 13 - iW;
      double tHres = (process[i1].p() - process[iW].p()).m2Calc();
      double uHres = (process[i1].p() - process[iZ].p()).m2Calc(); 
      double s3now = process[iW].m2();
      double s4now = process[iZ].m2();

      // Kinematics combinations (norm(x) = |x|^2).
      double fGK135 = norm(fGK( 1, 2, 3, 4, 5, 6) - fGK( 1, 2, 5, 6, 3, 4) );
      double fGK136 = norm(fGK( 1, 2, 3, 4, 6, 5) - fGK( 1, 2, 6, 5, 3, 4) );
      double xiT    = xiGK( tHres, uHres, s3now, s4now);
      double xiU    = xiGK( uHres, tHres, s3now, s4now);
      double xjTU   = xjGK( tHres, uHres, s3now, s4now);

      //  Couplings of outgoing fermion from Z. Combine with kinematics.
      int idAbs     = process[i5].idAbs();
      double lfZ    = couplingsPtr->lf(idAbs); 
      double rfZ    = couplingsPtr->rf(idAbs); 
      wt            = lfZ*lfZ * fGK135 + rfZ*rfZ * fGK136;
      wtMax         = 4. * s3now * s4now * (lfZ*lfZ + rfZ*rfZ) 
                    * (xiT + xiU - xjTU);
   
    // Decay distribution like in f fbar -> H^+- -> W+- Z0. 
    } else {
      double p35  = 2. * process[i3].p() * process[i5].p(); 
      double p46  = 2. * process[i4].p() * process[i6].p(); 
      wt          = 16. * p35 * p46;
      wtMax       = sH2;
    }
  }

  // Angular weight in top decay by standard routine.
  else if (process[process[iResBeg].mother1()].idAbs() == 6)
    return weightTopDecay( process, iResBeg, iResEnd);
 
  // Done.
  return (wt / wtMax);

}


//==========================================================================

// Sigma1ffbar2Rhorizontal class.
// Cross section for f fbar' -> R^0 (f is a quark or lepton). 

//--------------------------------------------------------------------------

// Initialize process. 
  
void Sigma1ffbar2Rhorizontal::initProc() {

  // Store R^0 mass and width for propagator. 
  mRes     = particleDataPtr->m0(41);
  GammaRes = particleDataPtr->mWidth(41);
  m2Res    = mRes*mRes;
  GamMRat  = GammaRes / mRes;
  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());

  // Set pointer to particle properties and decay table.
  particlePtr = particleDataPtr->particleDataEntryPtr(41);

} 

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat), part independent of incoming flavour. 

void Sigma1ffbar2Rhorizontal::sigmaKin() { 

  // Set up Breit-Wigner. Cross section for W+ and W- separately.
  double sigBW  = 12. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) ); 
  double preFac = alpEM * thetaWRat * mH;
  sigma0Pos     = preFac * sigBW * particlePtr->resWidthOpen(41, mH);    
  sigma0Neg     = preFac * sigBW * particlePtr->resWidthOpen(-41, mH);    

}

//--------------------------------------------------------------------------

// Evaluate sigmaHat(sHat), including incoming flavour dependence. 

double Sigma1ffbar2Rhorizontal::sigmaHat() {

  // Check for allowed flavour combinations, one generation apart.
  if (id1 * id2 > 0 || abs(id1 + id2) != 2) return 0.;  

  // Find whether R0 or R0bar. Colour factors.
  double sigma = (id1 + id2 > 0) ? sigma0Pos : sigma0Neg;
  if (abs(id1) < 7) sigma /= 3.;

  // Answer.
  return sigma;    

}

//--------------------------------------------------------------------------

// Select identity, colour and anticolour.

void Sigma1ffbar2Rhorizontal::setIdColAcol() {

  // Outgoing R0 or R0bar.
  id3 = (id1 +id2 > 0) ? 41 : -41;
  setId( id1, id2, id3);

  // Colour flow topologies. Swap when antiquarks.
  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
  else              setColAcol( 0, 0, 0, 0, 0, 0);
  if (id1 < 0) swapColAcol();

}

//==========================================================================

} // end namespace Pythia8