EmissionVetoHook1.cc

Go to the documentation of this file.

#include "Pythia8/Pythia.h"

using namespace Pythia8;

#include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
void EmissionVetoHook1::fatalEmissionVeto(std::string message) {
  throw edm::Exception(edm::errors::Configuration, "Pythia8Interface") << "EmissionVeto: " << message << std::endl;
}

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

// Routines to calculate the pT (according to pTdefMode) in a splitting:
//   ISR: i (radiator after)  -> j (emitted after) k (radiator before)
//   FSR: i (radiator before) -> j (emitted after) k (radiator after)
// For the Pythia pT definition, a recoiler (after) must be specified.

// Compute the Pythia pT separation. Based on pTLund function in History.cc
double EmissionVetoHook1::pTpythia(
    const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR) {
  // Convenient shorthands for later
  Pythia8::Vec4 radVec = e[RadAfterBranch].p();
  Pythia8::Vec4 emtVec = e[EmtAfterBranch].p();
  Pythia8::Vec4 recVec = e[RecAfterBranch].p();
  int radID = e[RadAfterBranch].id();

  // Calculate virtuality of splitting
  double sign = (FSR) ? 1. : -1.;
  Pythia8::Vec4 Q(radVec + sign * emtVec);
  double Qsq = sign * Q.m2Calc();

  // Mass term of radiator
  double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? Pythia8::pow2(particleDataPtr->m0(radID)) : 0.;

  // z values for FSR and ISR
  double z, pTnow;
  if (FSR) {
    // Construct 2 -> 3 variables
    Pythia8::Vec4 sum = radVec + recVec + emtVec;
    double m2Dip = sum.m2Calc();
    double x1 = 2. * (sum * radVec) / m2Dip;
    double x3 = 2. * (sum * emtVec) / m2Dip;
    z = x1 / (x1 + x3);
    pTnow = z * (1. - z);

  } else {
    // Construct dipoles before/after splitting
    Pythia8::Vec4 qBR(radVec - emtVec + recVec);
    Pythia8::Vec4 qAR(radVec + recVec);
    z = qBR.m2Calc() / qAR.m2Calc();
    pTnow = (1. - z);
  }

  // Virtuality with correct sign
  pTnow *= (Qsq - sign * m2Rad);

  // Can get negative pT for massive splittings
  if (pTnow < 0.) {
    std::cout << "Warning: pTpythia was negative" << std::endl;
    return -1.;
  }

#ifdef DBGOUTPUT
  std::cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " << EmtAfterBranch << ", rec = " << RecAfterBranch
            << ", pTnow = " << sqrt(pTnow) << std::endl;
#endif

  // Return pT
  return sqrt(pTnow);
}

// Compute the POWHEG pT separation between i and j
double EmissionVetoHook1::pTpowheg(const Pythia8::Event &e, int i, int j, bool FSR) {
  // pT value for FSR and ISR
  double pTnow = 0.;
  if (FSR) {
    // POWHEG d_ij (in CM frame). Note that the incoming beams have not
    // been updated in the parton systems pointer yet (i.e. prior to any
    // potential recoil).
    int iInA = partonSystemsPtr->getInA(0);
    int iInB = partonSystemsPtr->getInB(0);
    double betaZ = -(e[iInA].pz() + e[iInB].pz()) / (e[iInA].e() + e[iInB].e());
    Pythia8::Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
    iVecBst.bst(0., 0., betaZ);
    jVecBst.bst(0., 0., betaZ);
    pTnow = sqrt((iVecBst + jVecBst).m2Calc() * iVecBst.e() * jVecBst.e() / Pythia8::pow2(iVecBst.e() + jVecBst.e()));

  } else {
    // POWHEG pT_ISR is just kinematic pT
    pTnow = e[j].pT();
  }

  // Check result
  if (pTnow < 0.) {
    std::cout << "Warning: pTpowheg was negative" << std::endl;
    return -1.;
  }

#ifdef DBGOUTPUT
  std::cout << "pTpowheg: i = " << i << ", j = " << j << ", pTnow = " << pTnow << std::endl;
#endif

  return pTnow;
}

// Calculate pT for a splitting based on pTdefMode.
// If j is -1, all final-state partons are tried.
// If i, k, r and xSR are -1, then all incoming and outgoing
// partons are tried.
// xSR set to 0 means ISR, while xSR set to 1 means FSR
double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, int r, int xSRin) {
  // Loop over ISR and FSR if necessary
  double pTemt = -1., pTnow;
  int xSR1 = (xSRin == -1) ? 0 : xSRin;
  int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
  for (int xSR = xSR1; xSR < xSR2; xSR++) {
    // FSR flag
    bool FSR = (xSR == 0) ? false : true;

    // If all necessary arguments have been given, then directly calculate.
    // POWHEG ISR and FSR, need i and j.
    if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
      pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);

      // Pythia ISR, need i, j and r.
    } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
      pTemt = pTpythia(e, i, j, r, FSR);

      // Pythia FSR, need k, j and r.
    } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
      pTemt = pTpythia(e, k, j, r, FSR);

      // Otherwise need to try all possible combintations.
    } else {
      // Start by finding incoming legs to the hard system after
      // branching (radiator after branching, i for ISR).
      // Use partonSystemsPtr to find incoming just prior to the
      // branching and track mothers.
      int iInA = partonSystemsPtr->getInA(0);
      int iInB = partonSystemsPtr->getInB(0);
      while (e[iInA].mother1() != 1) {
        iInA = e[iInA].mother1();
      }
      while (e[iInB].mother1() != 2) {
        iInB = e[iInB].mother1();
      }

      // If we do not have j, then try all final-state partons
      int jNow = (j > 0) ? j : 0;
      int jMax = (j > 0) ? j + 1 : e.size();
      for (; jNow < jMax; jNow++) {
        // Final-state only
        if (!e[jNow].isFinal())
          continue;
        // Exclude photons (and W/Z!)
        if (QEDvetoMode == 0 && e[jNow].colType() == 0)
          continue;

        // POWHEG
        if (pTdefMode == 0 || pTdefMode == 1) {
          // ISR - only done once as just kinematical pT
          if (!FSR) {
            pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
            if (pTnow > 0.)
              pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);

            // FSR - try all outgoing partons from system before branching
            // as i. Note that for the hard system, there is no
            // "before branching" information.
          } else {
            int outSize = partonSystemsPtr->sizeOut(0);
            for (int iMem = 0; iMem < outSize; iMem++) {
              int iNow = partonSystemsPtr->getOut(0, iMem);

              // if i != jNow and no carbon copies
              if (iNow == jNow)
                continue;
              // Exlude photons (and W/Z!)
              if (QEDvetoMode == 0 && e[iNow].colType() == 0)
                continue;
              if (jNow == e[iNow].daughter1() && jNow == e[iNow].daughter2())
                continue;

              pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) ? false : FSR);
              if (pTnow > 0.)
                pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
            }  // for (iMem)

          }  // if (!FSR)

          // Pythia
        } else if (pTdefMode == 2) {
          // ISR - other incoming as recoiler
          if (!FSR) {
            pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
            if (pTnow > 0.)
              pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
            pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
            if (pTnow > 0.)
              pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);

            // FSR - try all final-state coloured partons as radiator
            //       after emission (k).
          } else {
            for (int kNow = 0; kNow < e.size(); kNow++) {
              if (kNow == jNow || !e[kNow].isFinal())
                continue;
              if (QEDvetoMode == 0 && e[kNow].colType() == 0)
                continue;

              // For this kNow, need to have a recoiler.
              // Try two incoming.
              pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
              if (pTnow > 0.)
                pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
              pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
              if (pTnow > 0.)
                pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);

              // Try all other outgoing.
              for (int rNow = 0; rNow < e.size(); rNow++) {
                if (rNow == kNow || rNow == jNow || !e[rNow].isFinal())
                  continue;
                if (QEDvetoMode == 0 && e[rNow].colType() == 0)
                  continue;
                pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
                if (pTnow > 0.)
                  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
              }  // for (rNow)

            }  // for (kNow)
          }    // if (!FSR)
        }      // if (pTdefMode)
      }        // for (j)
    }
  }  // for (xSR)

#ifdef DBGOUTPUT
  std::cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k << ", r = " << r << ", xSR = " << xSRin
            << ", pTemt = " << pTemt << std::endl;
#endif

  return pTemt;
}

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

// Extraction of pThard based on the incoming event.
// Assume that all the final-state particles are in a continuous block
// at the end of the event and the final entry is the POWHEG emission.
// If there is no POWHEG emission, then pThard is set to QRen.
bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) {
  if (nFinalMode == 3 && pThardMode != 0)
    fatalEmissionVeto(std::string(
        "When nFinalMode is set to 3, ptHardMode should be set to 0, since the emission variables in doVetoMPIStep are "
        "not set correctly case when there are three possible particle Born particle counts."));

  // Extra check on nMPI
  if (nMPI > 1)
    return false;

  // Find if there is a POWHEG emission. Go backwards through the
  // event record until there is a non-final particle. Also sum pT and
  // find pT_1 for possible MPI vetoing
  // Flag if POWHEG radiation is present and index at the same time
  int count = 0, inonfinal = 0;
  double pT1 = 0., pTsum = 0.;
  isEmt = false;
  int iEmt = -1;

  for (int i = e.size() - 1; i > 0; i--) {
    inonfinal = i;
    if (e[i].isFinal()) {
      count++;
      pT1 = e[i].pT();
      pTsum += e[i].pT();
      // the following was added for bbbar4l and will be triggered by specifying nfinalmode == 2
      // the solution provided by Tomas may not be process independent but should work for hvq and bb4l
      // if the particle is a light quark or a gluon and
      // comes for a light quark or a gluon (not a resonance, not a top)
      // then it is the POWHEG emission (hvq) or the POWHEG emission in production (bb4l)
      if (nFinalMode == 2) {
        if ((abs(e[i].id()) < 6 || e[i].id() == 21) &&
            (abs(e[e[i].mother1()].id()) < 6 || e[e[i].mother1()].id() == 21)) {
          isEmt = true;
          iEmt = i;
        }
      }
    } else
      break;
  }

  nFinal = nFinalExt;
  if (nFinal < 0 || nFinalMode == 1) {  // nFinal is not specified from external, try to find out
    int first = -1, myid;
    int last = -1;
    for (int ip = 2; ip < e.size(); ip++) {
      myid = e[ip].id();
      if (abs(myid) < 6 || abs(myid) == 21)
        continue;
      first = ip;
      break;
    }
    if (first < 0)
      fatalEmissionVeto(std::string("signal particles not found"));
    for (int ip = first; ip < e.size(); ip++) {
      myid = e[ip].id();
      if (abs(myid) < 6 || abs(myid) == 21)
        continue;
      last = ip;
    }
    nFinal = last - inonfinal;
  }

  // don't perform a cross check in case of nfinalmode == 2
  if (nFinalMode != 2) {
    // Extra check that we have the correct final state
    // In POWHEG WGamma, both w+0/1 jets and w+gamma+0/1 jets are generated at the same time, which leads to three different possible numbers of particles
    // Normally, there would be only two possible numbers of particles for X and X+1 jet
    if (nFinalMode == 3) {
      if (count != nFinal && count != nFinal + 1 && count != nFinal - 1)
        fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
    } else {
      if (count != nFinal && count != nFinal + 1)
        fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
    }
    // Flag if POWHEG radiation present and index
    if (count == nFinal + 1)
      isEmt = true;
    if (isEmt)
      iEmt = e.size() - 1;
  }

  // If there is no radiation or if pThardMode is 0 then set pThard to QRen.
  if (!isEmt || pThardMode == 0) {
    pThard = infoPtr->QRen();

    // If pThardMode is 1 then the pT of the POWHEG emission is checked against
    // all other incoming and outgoing partons, with the minimal value taken
  } else if (pThardMode == 1) {
    pThard = pTcalc(e, -1, iEmt, -1, -1, -1);

    // If pThardMode is 2, then the pT of all final-state partons is checked
    // against all other incoming and outgoing partons, with the minimal value
    // taken
  } else if (pThardMode == 2) {
    pThard = pTcalc(e, -1, -1, -1, -1, -1);
  }

  // Find MPI veto pT if necessary
  if (MPIvetoOn) {
    pTMPI = (isEmt) ? pTsum / 2. : pT1;
  }

  if (Verbosity)
    std::cout << "doVetoMPIStep: QFac = " << infoPtr->QFac() << ", QRen = " << infoPtr->QRen()
              << ", pThard = " << pThard << std::endl
              << std::endl;

  // Initialise other variables
  accepted = false;
  nAcceptSeq = 0;  // should not  reset nISRveto = nFSRveto = nFSRvetoBB4l here

  // Do not veto the event
  return false;
}

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

// ISR veto

bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys) {
  // Must be radiation from the hard system
  if (iSys != 0)
    return false;

  // If we already have accepted 'vetoCount' emissions in a row, do nothing
  if (vetoOn && nAcceptSeq >= vetoCount)
    return false;

  // Pythia radiator after, emitted and recoiler after.
  int iRadAft = -1, iEmt = -1, iRecAft = -1;
  for (int i = e.size() - 1; i > 0; i--) {
    if (iRadAft == -1 && e[i].status() == -41)
      iRadAft = i;
    else if (iEmt == -1 && e[i].status() == 43)
      iEmt = i;
    else if (iRecAft == -1 && e[i].status() == -42)
      iRecAft = i;
    if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
      break;
  }
  if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
    e.list();
    fatalEmissionVeto(std::string("Couldn't find Pythia ISR emission"));
  }

  // pTemtMode == 0: pT of emitted w.r.t. radiator
  // pTemtMode == 1: std::min(pT of emitted w.r.t. all incoming/outgoing)
  // pTemtMode == 2: std::min(pT of all outgoing w.r.t. all incoming/outgoing)
  int xSR = (pTemtMode == 0) ? 0 : -1;
  int i = (pTemtMode == 0) ? iRadAft : -1;
  int j = (pTemtMode != 2) ? iEmt : -1;
  int k = -1;
  int r = (pTemtMode == 0) ? iRecAft : -1;
  double pTemt = pTcalc(e, i, j, k, r, xSR);

#ifdef DBGOUTPUT
  std::cout << "doVetoISREmission: pTemt = " << pTemt << std::endl << std::endl;
#endif

  // If a Born configuration, and a colorless FS, and QEDvetoMode=2,
  // then don't veto photons, W, or Z harder than pThard
  bool vetoParton = (!isEmt && e[iEmt].colType() == 0 && QEDvetoMode == 2) ? false : true;

  // Veto if pTemt > pThard
  if (pTemt > pThard) {
    if (!vetoParton) {
      // Don't veto ANY emissions afterwards
      nAcceptSeq = vetoCount - 1;
    } else {
      nAcceptSeq = 0;
      nISRveto++;
      return true;
    }
  }

  // Else mark that an emission has been accepted and continue
  nAcceptSeq++;
  accepted = true;
  return false;
}

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

// FSR veto

bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys, bool inResonance) {
  // Must be radiation from the hard system
  if (iSys != 0)
    return false;

  // only use for outside resonance vetos in combination with bb4l:FSREmission:veto
  if (inResonance && settingsPtr->flag("POWHEG:bb4l:FSREmission:veto") == 1)
    return false;

  // If we already have accepted 'vetoCount' emissions in a row, do nothing
  if (vetoOn && nAcceptSeq >= vetoCount)
    return false;

  // Pythia radiator (before and after), emitted and recoiler (after)
  int iRecAft = e.size() - 1;
  int iEmt = e.size() - 2;
  int iRadAft = e.size() - 3;
  int iRadBef = e[iEmt].mother1();
  if ((e[iRecAft].status() != 52 && e[iRecAft].status() != -53) || e[iEmt].status() != 51 ||
      e[iRadAft].status() != 51) {
    e.list();
    fatalEmissionVeto(std::string("Couldn't find Pythia FSR emission"));
  }

  // Behaviour based on pTemtMode:
  //  0 - pT of emitted w.r.t. radiator before
  //  1 - std::min(pT of emitted w.r.t. all incoming/outgoing)
  //  2 - std::min(pT of all outgoing w.r.t. all incoming/outgoing)
  int xSR = (pTemtMode == 0) ? 1 : -1;
  int i = (pTemtMode == 0) ? iRadBef : -1;
  int k = (pTemtMode == 0) ? iRadAft : -1;
  int r = (pTemtMode == 0) ? iRecAft : -1;

  // When pTemtMode is 0 or 1, iEmt has been selected
  double pTemt = -1.;
  if (pTemtMode == 0 || pTemtMode == 1) {
    // Which parton is emitted, based on emittedMode:
    //  0 - Pythia definition of emitted
    //  1 - Pythia definition of radiated after emission
    //  2 - Random selection of emitted or radiated after emission
    //  3 - Try both emitted and radiated after emission
    int j = iRadAft;
    if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5))
      j++;

    for (int jLoop = 0; jLoop < 2; jLoop++) {
      if (jLoop == 0)
        pTemt = pTcalc(e, i, j, k, r, xSR);
      else if (jLoop == 1)
        pTemt = std::min(pTemt, pTcalc(e, i, j, k, r, xSR));

      // For emittedMode == 3, have tried iRadAft, now try iEmt
      if (emittedMode != 3)
        break;
      if (k != -1)
        std::swap(j, k);
      else
        j = iEmt;
    }

    // If pTemtMode is 2, then try all final-state partons as emitted
  } else if (pTemtMode == 2) {
    pTemt = pTcalc(e, i, -1, k, r, xSR);
  }

#ifdef DBGOUTPUT
  std::cout << "doVetoFSREmission: pTemt = " << pTemt << std::endl << std::endl;
#endif

  // If a Born configuration, and a colorless FS, and QEDvetoMode=2,
  // then don't veto photons, W, or Z harder than pThard
  bool vetoParton = (!isEmt && e[iEmt].colType() == 0 && QEDvetoMode == 2) ? false : true;

  // Veto if pTemt > pThard
  if (pTemt > pThard) {
    if (!vetoParton) {
      // Don't veto ANY emissions afterwards
      nAcceptSeq = vetoCount - 1;
    } else {
      nAcceptSeq = 0;
      nFSRveto++;
      return true;
    }
  }

  // Else mark that an emission has been accepted and continue
  nAcceptSeq++;
  accepted = true;
  return false;
}

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

// MPI veto

bool EmissionVetoHook1::doVetoMPIEmission(int, const Pythia8::Event &e) {
  if (MPIvetoOn) {
    if (e[e.size() - 1].pT() > pTMPI) {
#ifdef DBGOUTPUT
      std::cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() << ", pTMPI = " << pTMPI << std::endl
                << std::endl;
#endif
      return true;
    }
  }
  return false;
}