Pythia8Interface/plugins/EmissionVetoHook1.cc

0001 #include "Pythia8/Pythia.h"
0002
0003 using namespace Pythia8;
0004
0005 #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h"
0006
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 void EmissionVetoHook1::fatalEmissionVeto(std::string message) {
0009   throw edm::Exception(edm::errors::Configuration, "Pythia8Interface") << "EmissionVeto: " << message << std::endl;
0010 }
0011
0012 //--------------------------------------------------------------------------
0013
0014 // Routines to calculate the pT (according to pTdefMode) in a splitting:
0015 //   ISR: i (radiator after)  -> j (emitted after) k (radiator before)
0016 //   FSR: i (radiator before) -> j (emitted after) k (radiator after)
0017 // For the Pythia pT definition, a recoiler (after) must be specified.
0018
0019 // Compute the Pythia pT separation. Based on pTLund function in History.cc
0020 double EmissionVetoHook1::pTpythia(
0021     const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR) {
0022   // Convenient shorthands for later
0023   Pythia8::Vec4 radVec = e[RadAfterBranch].p();
0024   Pythia8::Vec4 emtVec = e[EmtAfterBranch].p();
0025   Pythia8::Vec4 recVec = e[RecAfterBranch].p();
0026   int radID = e[RadAfterBranch].id();
0027
0028   // Calculate virtuality of splitting
0029   double sign = (FSR) ? 1. : -1.;
0030   Pythia8::Vec4 Q(radVec + sign * emtVec);
0031   double Qsq = sign * Q.m2Calc();
0032
0033   // Mass term of radiator
0034   double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? Pythia8::pow2(particleDataPtr->m0(radID)) : 0.;
0035
0036   // z values for FSR and ISR
0037   double z, pTnow;
0038   if (FSR) {
0039     // Construct 2 -> 3 variables
0040     Pythia8::Vec4 sum = radVec + recVec + emtVec;
0041     double m2Dip = sum.m2Calc();
0042     double x1 = 2. * (sum * radVec) / m2Dip;
0043     double x3 = 2. * (sum * emtVec) / m2Dip;
0044     z = x1 / (x1 + x3);
0045     pTnow = z * (1. - z);
0046
0047   } else {
0048     // Construct dipoles before/after splitting
0049     Pythia8::Vec4 qBR(radVec - emtVec + recVec);
0050     Pythia8::Vec4 qAR(radVec + recVec);
0051     z = qBR.m2Calc() / qAR.m2Calc();
0052     pTnow = (1. - z);
0053   }
0054
0055   // Virtuality with correct sign
0056   pTnow *= (Qsq - sign * m2Rad);
0057
0058   // Can get negative pT for massive splittings
0059   if (pTnow < 0.) {
0060     std::cout << "Warning: pTpythia was negative" << std::endl;
0061     return -1.;
0062   }
0063
0064 #ifdef DBGOUTPUT
0065   std::cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " << EmtAfterBranch << ", rec = " << RecAfterBranch
0066             << ", pTnow = " << sqrt(pTnow) << std::endl;
0067 #endif
0068
0069   // Return pT
0070   return sqrt(pTnow);
0071 }
0072
0073 // Compute the POWHEG pT separation between i and j
0074 double EmissionVetoHook1::pTpowheg(const Pythia8::Event &e, int i, int j, bool FSR) {
0075   // pT value for FSR and ISR
0076   double pTnow = 0.;
0077   if (FSR) {
0078     // POWHEG d_ij (in CM frame). Note that the incoming beams have not
0079     // been updated in the parton systems pointer yet (i.e. prior to any
0080     // potential recoil).
0081     int iInA = partonSystemsPtr->getInA(0);
0082     int iInB = partonSystemsPtr->getInB(0);
0083     double betaZ = -(e[iInA].pz() + e[iInB].pz()) / (e[iInA].e() + e[iInB].e());
0084     Pythia8::Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
0085     iVecBst.bst(0., 0., betaZ);
0086     jVecBst.bst(0., 0., betaZ);
0087     pTnow = sqrt((iVecBst + jVecBst).m2Calc() * iVecBst.e() * jVecBst.e() / Pythia8::pow2(iVecBst.e() + jVecBst.e()));
0088
0089   } else {
0090     // POWHEG pT_ISR is just kinematic pT
0091     pTnow = e[j].pT();
0092   }
0093
0094   // Check result
0095   if (pTnow < 0.) {
0096     std::cout << "Warning: pTpowheg was negative" << std::endl;
0097     return -1.;
0098   }
0099
0100 #ifdef DBGOUTPUT
0101   std::cout << "pTpowheg: i = " << i << ", j = " << j << ", pTnow = " << pTnow << std::endl;
0102 #endif
0103
0104   return pTnow;
0105 }
0106
0107 // Calculate pT for a splitting based on pTdefMode.
0108 // If j is -1, all final-state partons are tried.
0109 // If i, k, r and xSR are -1, then all incoming and outgoing
0110 // partons are tried.
0111 // xSR set to 0 means ISR, while xSR set to 1 means FSR
0112 double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, int r, int xSRin) {
0113   // Loop over ISR and FSR if necessary
0114   double pTemt = -1., pTnow;
0115   int xSR1 = (xSRin == -1) ? 0 : xSRin;
0116   int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
0117   for (int xSR = xSR1; xSR < xSR2; xSR++) {
0118     // FSR flag
0119     bool FSR = (xSR == 0) ? false : true;
0120
0121     // If all necessary arguments have been given, then directly calculate.
0122     // POWHEG ISR and FSR, need i and j.
0123     if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
0124       pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
0125
0126       // Pythia ISR, need i, j and r.
0127     } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
0128       pTemt = pTpythia(e, i, j, r, FSR);
0129
0130       // Pythia FSR, need k, j and r.
0131     } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
0132       pTemt = pTpythia(e, k, j, r, FSR);
0133
0134       // Otherwise need to try all possible combintations.
0135     } else {
0136       // Start by finding incoming legs to the hard system after
0137       // branching (radiator after branching, i for ISR).
0138       // Use partonSystemsPtr to find incoming just prior to the
0139       // branching and track mothers.
0140       int iInA = partonSystemsPtr->getInA(0);
0141       int iInB = partonSystemsPtr->getInB(0);
0142       while (e[iInA].mother1() != 1) {
0143         iInA = e[iInA].mother1();
0144       }
0145       while (e[iInB].mother1() != 2) {
0146         iInB = e[iInB].mother1();
0147       }
0148
0149       // If we do not have j, then try all final-state partons
0150       int jNow = (j > 0) ? j : 0;
0151       int jMax = (j > 0) ? j + 1 : e.size();
0152       for (; jNow < jMax; jNow++) {
0153         // Final-state only
0154         if (!e[jNow].isFinal())
0155           continue;
0156         // Exclude photons (and W/Z!)
0157         if (QEDvetoMode == 0 && e[jNow].colType() == 0)
0158           continue;
0159
0160         // POWHEG
0161         if (pTdefMode == 0 || pTdefMode == 1) {
0162           // ISR - only done once as just kinematical pT
0163           if (!FSR) {
0164             pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
0165             if (pTnow > 0.)
0166               pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0167
0168             // FSR - try all outgoing partons from system before branching
0169             // as i. Note that for the hard system, there is no
0170             // "before branching" information.
0171           } else {
0172             int outSize = partonSystemsPtr->sizeOut(0);
0173             for (int iMem = 0; iMem < outSize; iMem++) {
0174               int iNow = partonSystemsPtr->getOut(0, iMem);
0175
0176               // if i != jNow and no carbon copies
0177               if (iNow == jNow)
0178                 continue;
0179               // Exlude photons (and W/Z!)
0180               if (QEDvetoMode == 0 && e[iNow].colType() == 0)
0181                 continue;
0182               if (jNow == e[iNow].daughter1() && jNow == e[iNow].daughter2())
0183                 continue;
0184
0185               pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) ? false : FSR);
0186               if (pTnow > 0.)
0187                 pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0188             }  // for (iMem)
0189
0190           }  // if (!FSR)
0191
0192           // Pythia
0193         } else if (pTdefMode == 2) {
0194           // ISR - other incoming as recoiler
0195           if (!FSR) {
0196             pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
0197             if (pTnow > 0.)
0198               pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0199             pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
0200             if (pTnow > 0.)
0201               pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0202
0203             // FSR - try all final-state coloured partons as radiator
0204             //       after emission (k).
0205           } else {
0206             for (int kNow = 0; kNow < e.size(); kNow++) {
0207               if (kNow == jNow || !e[kNow].isFinal())
0208                 continue;
0209               if (QEDvetoMode == 0 && e[kNow].colType() == 0)
0210                 continue;
0211
0212               // For this kNow, need to have a recoiler.
0213               // Try two incoming.
0214               pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
0215               if (pTnow > 0.)
0216                 pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0217               pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
0218               if (pTnow > 0.)
0219                 pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0220
0221               // Try all other outgoing.
0222               for (int rNow = 0; rNow < e.size(); rNow++) {
0223                 if (rNow == kNow || rNow == jNow || !e[rNow].isFinal())
0224                   continue;
0225                 if (QEDvetoMode == 0 && e[rNow].colType() == 0)
0226                   continue;
0227                 pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
0228                 if (pTnow > 0.)
0229                   pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
0230               }  // for (rNow)
0231
0232             }  // for (kNow)
0233           }  // if (!FSR)
0234         }  // if (pTdefMode)
0235       }  // for (j)
0236     }
0237   }  // for (xSR)
0238
0239 #ifdef DBGOUTPUT
0240   std::cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k << ", r = " << r << ", xSR = " << xSRin
0241             << ", pTemt = " << pTemt << std::endl;
0242 #endif
0243
0244   return pTemt;
0245 }
0246
0247 //--------------------------------------------------------------------------
0248
0249 // Extraction of pThard based on the incoming event.
0250 // Assume that all the final-state particles are in a continuous block
0251 // at the end of the event and the final entry is the POWHEG emission.
0252 // If there is no POWHEG emission, then pThard is set to QRen.
0253 bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) {
0254   if (nFinalMode == 3 && pThardMode != 0)
0255     fatalEmissionVeto(std::string(
0256         "When nFinalMode is set to 3, ptHardMode should be set to 0, since the emission variables in doVetoMPIStep are "
0257         "not set correctly case when there are three possible particle Born particle counts."));
0258
0259   // Extra check on nMPI
0260   if (nMPI > 1)
0261     return false;
0262
0263   // Find if there is a POWHEG emission. Go backwards through the
0264   // event record until there is a non-final particle. Also sum pT and
0265   // find pT_1 for possible MPI vetoing
0266   // Flag if POWHEG radiation is present and index at the same time
0267   int count = 0, inonfinal = 0;
0268   double pT1 = 0., pTsum = 0.;
0269   isEmt = false;
0270   int iEmt = -1;
0271
0272   for (int i = e.size() - 1; i > 0; i--) {
0273     inonfinal = i;
0274     if (e[i].isFinal()) {
0275       count++;
0276       pT1 = e[i].pT();
0277       pTsum += e[i].pT();
0278       // the following was added for bbbar4l and will be triggered by specifying nfinalmode == 2
0279       // the solution provided by Tomas may not be process independent but should work for hvq and bb4l
0280       // if the particle is a light quark or a gluon and
0281       // comes for a light quark or a gluon (not a resonance, not a top)
0282       // then it is the POWHEG emission (hvq) or the POWHEG emission in production (bb4l)
0283       if (nFinalMode == 2) {
0284         if ((abs(e[i].id()) < 6 || e[i].id() == 21) &&
0285             (abs(e[e[i].mother1()].id()) < 6 || e[e[i].mother1()].id() == 21)) {
0286           isEmt = true;
0287           iEmt = i;
0288         }
0289       }
0290     } else
0291       break;
0292   }
0293
0294   nFinal = nFinalExt;
0295   if (nFinal < 0 || nFinalMode == 1) {  // nFinal is not specified from external, try to find out
0296     int first = -1, myid;
0297     int last = -1;
0298     for (int ip = 2; ip < e.size(); ip++) {
0299       myid = e[ip].id();
0300       if (abs(myid) < 6 || abs(myid) == 21)
0301         continue;
0302       first = ip;
0303       break;
0304     }
0305     if (first < 0)
0306       fatalEmissionVeto(std::string("signal particles not found"));
0307     for (int ip = first; ip < e.size(); ip++) {
0308       myid = e[ip].id();
0309       if (abs(myid) < 6 || abs(myid) == 21)
0310         continue;
0311       last = ip;
0312     }
0313     nFinal = last - inonfinal;
0314   }
0315
0316   // don't perform a cross check in case of nfinalmode == 2
0317   if (nFinalMode != 2) {
0318     // Extra check that we have the correct final state
0319     // 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
0320     // Normally, there would be only two possible numbers of particles for X and X+1 jet
0321     if (nFinalMode == 3) {
0322       if (count != nFinal && count != nFinal + 1 && count != nFinal - 1)
0323         fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
0324     } else {
0325       if (count != nFinal && count != nFinal + 1)
0326         fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
0327     }
0328     // Flag if POWHEG radiation present and index
0329     if (count == nFinal + 1)
0330       isEmt = true;
0331     if (isEmt)
0332       iEmt = e.size() - 1;
0333   }
0334
0335   // If there is no radiation or if pThardMode is 0 then set pThard to QRen.
0336   if (!isEmt || pThardMode == 0) {
0337     pThard = infoPtr->QRen();
0338
0339     // If pThardMode is 1 then the pT of the POWHEG emission is checked against
0340     // all other incoming and outgoing partons, with the minimal value taken
0341   } else if (pThardMode == 1) {
0342     pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
0343
0344     // If pThardMode is 2, then the pT of all final-state partons is checked
0345     // against all other incoming and outgoing partons, with the minimal value
0346     // taken
0347   } else if (pThardMode == 2) {
0348     pThard = pTcalc(e, -1, -1, -1, -1, -1);
0349   }
0350
0351   // Find MPI veto pT if necessary
0352   if (MPIvetoOn) {
0353     pTMPI = (isEmt) ? pTsum / 2. : pT1;
0354   }
0355
0356   if (Verbosity)
0357     std::cout << "doVetoMPIStep: QFac = " << infoPtr->QFac() << ", QRen = " << infoPtr->QRen()
0358               << ", pThard = " << pThard << std::endl
0359               << std::endl;
0360
0361   // Initialise other variables
0362   accepted = false;
0363   nAcceptSeq = 0;  // should not  reset nISRveto = nFSRveto = nFSRvetoBB4l here
0364
0365   // Do not veto the event
0366   return false;
0367 }
0368
0369 //--------------------------------------------------------------------------
0370
0371 // ISR veto
0372
0373 bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys) {
0374   // Must be radiation from the hard system
0375   if (iSys != 0)
0376     return false;
0377
0378   // If we already have accepted 'vetoCount' emissions in a row, do nothing
0379   if (vetoOn && nAcceptSeq >= vetoCount)
0380     return false;
0381
0382   // Pythia radiator after, emitted and recoiler after.
0383   int iRadAft = -1, iEmt = -1, iRecAft = -1;
0384   for (int i = e.size() - 1; i > 0; i--) {
0385     if (iRadAft == -1 && e[i].status() == -41)
0386       iRadAft = i;
0387     else if (iEmt == -1 && e[i].status() == 43)
0388       iEmt = i;
0389     else if (iRecAft == -1 && e[i].status() == -42)
0390       iRecAft = i;
0391     if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
0392       break;
0393   }
0394   if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
0395     e.list();
0396     fatalEmissionVeto(std::string("Couldn't find Pythia ISR emission"));
0397   }
0398
0399   // pTemtMode == 0: pT of emitted w.r.t. radiator
0400   // pTemtMode == 1: std::min(pT of emitted w.r.t. all incoming/outgoing)
0401   // pTemtMode == 2: std::min(pT of all outgoing w.r.t. all incoming/outgoing)
0402   int xSR = (pTemtMode == 0) ? 0 : -1;
0403   int i = (pTemtMode == 0) ? iRadAft : -1;
0404   int j = (pTemtMode != 2) ? iEmt : -1;
0405   int k = -1;
0406   int r = (pTemtMode == 0) ? iRecAft : -1;
0407   double pTemt = pTcalc(e, i, j, k, r, xSR);
0408
0409 #ifdef DBGOUTPUT
0410   std::cout << "doVetoISREmission: pTemt = " << pTemt << std::endl << std::endl;
0411 #endif
0412
0413   // If a Born configuration, and a colorless FS, and QEDvetoMode=2,
0414   // then don't veto photons, W, or Z harder than pThard
0415   bool vetoParton = (!isEmt && e[iEmt].colType() == 0 && QEDvetoMode == 2) ? false : true;
0416
0417   // Veto if pTemt > pThard
0418   if (pTemt > pThard) {
0419     if (!vetoParton) {
0420       // Don't veto ANY emissions afterwards
0421       nAcceptSeq = vetoCount - 1;
0422     } else {
0423       nAcceptSeq = 0;
0424       nISRveto++;
0425       return true;
0426     }
0427   }
0428
0429   // Else mark that an emission has been accepted and continue
0430   nAcceptSeq++;
0431   accepted = true;
0432   return false;
0433 }
0434
0435 //--------------------------------------------------------------------------
0436
0437 // FSR veto
0438
0439 bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys, bool inResonance) {
0440   // Must be radiation from the hard system
0441   if (iSys != 0)
0442     return false;
0443
0444   // only use for outside resonance vetos in combination with bb4l:FSREmission:veto
0445   if (inResonance && settingsPtr->flag("POWHEG:bb4l:FSREmission:veto") == 1)
0446     return false;
0447
0448   // If we already have accepted 'vetoCount' emissions in a row, do nothing
0449   if (vetoOn && nAcceptSeq >= vetoCount)
0450     return false;
0451
0452   // Pythia radiator (before and after), emitted and recoiler (after)
0453   int iRecAft = e.size() - 1;
0454   int iEmt = e.size() - 2;
0455   int iRadAft = e.size() - 3;
0456   int iRadBef = e[iEmt].mother1();
0457   if ((e[iRecAft].status() != 52 && e[iRecAft].status() != -53) || e[iEmt].status() != 51 ||
0458       e[iRadAft].status() != 51) {
0459     e.list();
0460     fatalEmissionVeto(std::string("Couldn't find Pythia FSR emission"));
0461   }
0462
0463   // Behaviour based on pTemtMode:
0464   //  0 - pT of emitted w.r.t. radiator before
0465   //  1 - std::min(pT of emitted w.r.t. all incoming/outgoing)
0466   //  2 - std::min(pT of all outgoing w.r.t. all incoming/outgoing)
0467   int xSR = (pTemtMode == 0) ? 1 : -1;
0468   int i = (pTemtMode == 0) ? iRadBef : -1;
0469   int k = (pTemtMode == 0) ? iRadAft : -1;
0470   int r = (pTemtMode == 0) ? iRecAft : -1;
0471
0472   // When pTemtMode is 0 or 1, iEmt has been selected
0473   double pTemt = -1.;
0474   if (pTemtMode == 0 || pTemtMode == 1) {
0475     // Which parton is emitted, based on emittedMode:
0476     //  0 - Pythia definition of emitted
0477     //  1 - Pythia definition of radiated after emission
0478     //  2 - Random selection of emitted or radiated after emission
0479     //  3 - Try both emitted and radiated after emission
0480     int j = iRadAft;
0481     if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5))
0482       j++;
0483
0484     for (int jLoop = 0; jLoop < 2; jLoop++) {
0485       if (jLoop == 0)
0486         pTemt = pTcalc(e, i, j, k, r, xSR);
0487       else if (jLoop == 1)
0488         pTemt = std::min(pTemt, pTcalc(e, i, j, k, r, xSR));
0489
0490       // For emittedMode == 3, have tried iRadAft, now try iEmt
0491       if (emittedMode != 3)
0492         break;
0493       if (k != -1)
0494         std::swap(j, k);
0495       else
0496         j = iEmt;
0497     }
0498
0499     // If pTemtMode is 2, then try all final-state partons as emitted
0500   } else if (pTemtMode == 2) {
0501     pTemt = pTcalc(e, i, -1, k, r, xSR);
0502   }
0503
0504 #ifdef DBGOUTPUT
0505   std::cout << "doVetoFSREmission: pTemt = " << pTemt << std::endl << std::endl;
0506 #endif
0507
0508   // If a Born configuration, and a colorless FS, and QEDvetoMode=2,
0509   // then don't veto photons, W, or Z harder than pThard
0510   bool vetoParton = (!isEmt && e[iEmt].colType() == 0 && QEDvetoMode == 2) ? false : true;
0511
0512   // Veto if pTemt > pThard
0513   if (pTemt > pThard) {
0514     if (!vetoParton) {
0515       // Don't veto ANY emissions afterwards
0516       nAcceptSeq = vetoCount - 1;
0517     } else {
0518       nAcceptSeq = 0;
0519       nFSRveto++;
0520       return true;
0521     }
0522   }
0523
0524   // Else mark that an emission has been accepted and continue
0525   nAcceptSeq++;
0526   accepted = true;
0527   return false;
0528 }
0529
0530 //--------------------------------------------------------------------------
0531
0532 // MPI veto
0533
0534 bool EmissionVetoHook1::doVetoMPIEmission(int, const Pythia8::Event &e) {
0535   if (MPIvetoOn) {
0536     if (e[e.size() - 1].pT() > pTMPI) {
0537 #ifdef DBGOUTPUT
0538       std::cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() << ", pTMPI = " << pTMPI << std::endl
0539                 << std::endl;
0540 #endif
0541       return true;
0542     }
0543   }
0544   return false;
0545 }