Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-03-31 08:26:53

0001 // PowhegHooks.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2025 Richard Corke, Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // Author: Richard Corke, modified by Christian T Preuss.
0007 // This class is used to perform a vetoed shower, where emissions
0008 // already covered in a POWHEG NLO generator should be omitted.
0009 // To first approximation the handover should happen at the SCALE
0010 // of the LHA, but since the POWHEG-BOX uses a different pT definition
0011 // than PYTHIA, both for ISR and FSR, a more sophisticated treatment
0012 // is needed. See the online manual on POWHEG matching for details.
0013 
0014 #ifndef Pythia8_PowhegHooks_H
0015 #define Pythia8_PowhegHooks_H
0016 
0017 // Includes
0018 #include "Pythia8/Pythia.h"
0019 #include "Pythia8/Plugins.h"
0020 
0021 namespace Pythia8 {
0022 
0023 //==========================================================================
0024 
0025 // Use userhooks to veto PYTHIA emissions above the POWHEG scale.
0026 
0027 class PowhegHooks : public UserHooks {
0028 
0029 public:
0030 
0031   // Constructors and destructor.
0032   PowhegHooks() {}
0033   PowhegHooks(Pythia*, Settings*, Logger*) {}
0034   ~PowhegHooks() {}
0035 
0036   //--------------------------------------------------------------------------
0037 
0038   // Initialize settings, detailing merging strategy to use.
0039   bool initAfterBeams() {
0040     nFinal      = settingsPtr->mode("POWHEG:nFinal");
0041     vetoMode    = settingsPtr->mode("POWHEG:veto");
0042     vetoCount   = settingsPtr->mode("POWHEG:vetoCount");
0043     pThardMode  = settingsPtr->mode("POWHEG:pThard");
0044     pTemtMode   = settingsPtr->mode("POWHEG:pTemt");
0045     emittedMode = settingsPtr->mode("POWHEG:emitted");
0046     pTdefMode   = settingsPtr->mode("POWHEG:pTdef");
0047     MPIvetoMode = settingsPtr->mode("POWHEG:MPIveto");
0048     QEDvetoMode = settingsPtr->mode("POWHEG:QEDveto");
0049     showerModel = settingsPtr->mode("PartonShowers:model");
0050     return true;
0051   }
0052 
0053   //--------------------------------------------------------------------------
0054 
0055   // Routines to calculate the pT (according to pTdefMode) in a branching:
0056   //   ISR: i (radiator after)  -> j (emitted after) k (radiator before)
0057   //   FSR: i (radiator before) -> j (emitted after) k (radiator after)
0058   // For the Pythia pT definitions, a recoiler (after) must be specified.
0059 
0060   // Top-level wrapper for shower pTs.
0061   inline double pT(const Event& e, int RadAfterBranch,
0062     int EmtAfterBranch, int RecAfterBranch, bool FSR) {
0063     // VINCIA pT definition.
0064     if (showerModel == 2)
0065       return pTvincia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch);
0066     // Simple-shower pT definition.
0067     return pTpythia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch, FSR);
0068   }
0069 
0070   //--------------------------------------------------------------------------
0071 
0072   // Compute the Pythia pT separation.
0073   // Based on pTLund function in History.cc and identical to pTevol.
0074   inline double pTpythia(const Event& e, int RadAfterBranch,
0075     int EmtAfterBranch, int RecAfterBranch, bool FSR) {
0076 
0077     // Convenient shorthands for later
0078     Vec4 radVec = e[RadAfterBranch].p();
0079     Vec4 emtVec = e[EmtAfterBranch].p();
0080     Vec4 recVec = e[RecAfterBranch].p();
0081     int  radID  = e[RadAfterBranch].id();
0082 
0083     // Calculate virtuality of splitting
0084     double sign = (FSR) ? 1. : -1.;
0085     Vec4 Q(radVec + sign * emtVec);
0086     double Qsq = sign * Q.m2Calc();
0087 
0088     // Mass term of radiator
0089     double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
0090                    pow2(particleDataPtr->m0(radID)) : 0.;
0091 
0092     // z values for FSR and ISR
0093     double z, pTnow;
0094     if (FSR) {
0095       // Construct 2 -> 3 variables
0096       Vec4 sum = radVec + recVec + emtVec;
0097       double m2Dip = sum.m2Calc();
0098       double x1 = 2. * (sum * radVec) / m2Dip;
0099       double x3 = 2. * (sum * emtVec) / m2Dip;
0100       z     = x1 / (x1 + x3);
0101       pTnow = z * (1. - z);
0102 
0103     } else {
0104       // Construct dipoles before/after splitting
0105       Vec4 qBR(radVec - emtVec + recVec);
0106       Vec4 qAR(radVec + recVec);
0107       z     = qBR.m2Calc() / qAR.m2Calc();
0108       pTnow = (1. - z);
0109     }
0110 
0111     // Virtuality with correct sign
0112     pTnow *= (Qsq - sign * m2Rad);
0113 
0114     // Can get negative pT for massive splittings.
0115     if (pTnow < 0.) {
0116       loggerPtr->WARNING_MSG("negative pT");
0117       return -1.;
0118     }
0119 
0120     // Return pT
0121     return sqrt(pTnow);
0122   }
0123 
0124   //--------------------------------------------------------------------------
0125 
0126   // Compute the Vincia pT as in Eq. (2.63)-(2.66) in arXiv:2003.00702.
0127   // Branching is assumed to be {13} {23} -> 1 3 2.
0128   inline double pTvincia(const Event& event, int i1, int i3, int i2) {
0129 
0130     // Shorthands.
0131     Vec4 p1 = event[i1].p();
0132     Vec4 p3 = event[i3].p();
0133     Vec4 p2 = event[i2].p();
0134 
0135     // Fetch mothers of 1 and 2.
0136     int iMoth1 = event[i1].mother1();
0137     int iMoth2 = event[i2].mother1();
0138     if (iMoth1 == 0 || iMoth2 == 0) {
0139       loggerPtr->ABORT_MSG("could not find mothers of particles");
0140       exit(1);
0141     }
0142 
0143     // Invariants defined as in Eq. (5) in arXiv:2008.09468.
0144     double mMoth1Sq = event[iMoth1].m2();
0145     double mMoth2Sq = event[iMoth2].m2();
0146     double sgn1 = event[i1].isFinal() ? 1. : -1.;
0147     double sgn2 = event[i2].isFinal() ? 1. : -1.;
0148     double qSq13 = sgn1*(m2(sgn1*p1+p3) - mMoth1Sq);
0149     double qSq23 = sgn2*(m2(sgn2*p2+p3) - mMoth2Sq);
0150 
0151     // Normalisation as in Eq. (6) in arXiv:2008.09468.
0152     double sMax = -1.;
0153     if (event[i1].isFinal() && event[i2].isFinal()) {
0154       // FF.
0155       sMax = m2(p1+p2+p3) - mMoth1Sq - mMoth2Sq;
0156     } else if ((event[i1].isResonance() && event[i2].isFinal())
0157       || (!event[i1].isFinal() && event[i2].isFinal())) {
0158       // RF or IF.
0159       sMax = 2.*p1*p3 + 2.*p1*p2;
0160     } else if ((event[i1].isFinal() && event[i2].isResonance())
0161       || (event[i1].isFinal() && !event[i2].isFinal())) {
0162       // FR or FI.
0163       sMax = 2.*p2*p3 + 2.*p1*p2;
0164     } else if (!event[i1].isFinal() || !event[i2].isFinal()) {
0165       // II.
0166       sMax = 2.*p1*p2;
0167     } else {
0168       loggerPtr->ABORT_MSG("could not determine branching type");
0169       exit(1);
0170     }
0171 
0172     // Calculate pT2 as in Eq. (5) in arXiv:2008.09468.
0173     double pT2now = qSq13*qSq23/sMax;
0174 
0175     // Sanity check.
0176     if (pT2now < 0.) {
0177       loggerPtr->WARNING_MSG("negative pT");
0178       return -1.;
0179     }
0180 
0181     // Return pT.
0182     return sqrt(pT2now);
0183   }
0184 
0185   //--------------------------------------------------------------------------
0186 
0187   // Compute the POWHEG pT separation between i and j.
0188   inline double pTpowheg(const Event &e, int i, int j, bool FSR) {
0189 
0190     // pT value for FSR and ISR
0191     double pTnow = 0.;
0192     if (FSR) {
0193       // POWHEG d_ij (in CM frame). Note that the incoming beams have not
0194       // been updated in the parton systems pointer yet (i.e. prior to any
0195       // potential recoil).
0196       int iInA = partonSystemsPtr->getInA(0);
0197       int iInB = partonSystemsPtr->getInB(0);
0198       double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
0199                        ( e[iInA].e()  + e[iInB].e()  );
0200       Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
0201       iVecBst.bst(0., 0., betaZ);
0202       jVecBst.bst(0., 0., betaZ);
0203       pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
0204                     iVecBst.e() * jVecBst.e() /
0205                     pow2(iVecBst.e() + jVecBst.e()) );
0206 
0207     } else {
0208       // POWHEG pT_ISR is just kinematic pT
0209       pTnow = e[j].pT();
0210     }
0211 
0212     // Check result.
0213     if (pTnow < 0.) {
0214       loggerPtr->WARNING_MSG("negative pT");
0215       return -1.;
0216     }
0217 
0218     return pTnow;
0219   }
0220 
0221   //--------------------------------------------------------------------------
0222 
0223   // Calculate pT for a splitting based on pTdefMode.
0224   // If j is -1, all final-state partons are tried.
0225   // If i, k, r and xSR are -1, then all incoming and outgoing
0226   // partons are tried.
0227   // xSR set to 0 means ISR, while xSR set to 1 means FSR.
0228   inline double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
0229 
0230     // Loop over ISR and FSR if necessary
0231     double pTemt = -1., pTnow;
0232     int xSR1 = (xSRin == -1) ? 0 : xSRin;
0233     int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
0234     for (int xSR = xSR1; xSR < xSR2; xSR++) {
0235       // FSR flag
0236       bool FSR = (xSR == 0) ? false : true;
0237 
0238       // If all necessary arguments have been given, then directly calculate.
0239       // POWHEG ISR and FSR, need i and j.
0240       if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
0241         pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
0242 
0243       // Pythia ISR, need i, j and r.
0244       } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
0245         pTemt = pT(e, i, j, r, FSR);
0246 
0247       // Pythia FSR, need k, j and r.
0248       } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
0249         pTemt = pT(e, k, j, r, FSR);
0250 
0251       // Otherwise need to try all possible combinations.
0252       } else {
0253         // Start by finding incoming legs to the hard system after
0254         // branching (radiator after branching, i for ISR).
0255         // Use partonSystemsPtr to find incoming just prior to the
0256         // branching and track mothers.
0257         int iInA = partonSystemsPtr->getInA(0);
0258         int iInB = partonSystemsPtr->getInB(0);
0259         while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
0260         while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
0261 
0262         // If we do not have j, then try all final-state partons.
0263         int jNow = (j > 0) ? j : 0;
0264         int jMax = (j > 0) ? j + 1 : e.size();
0265         for (; jNow < jMax; jNow++) {
0266 
0267           // Final-state only.
0268           if (!e[jNow].isFinal()) continue;
0269           // Exclude photons (and W/Z!)
0270           if (QEDvetoMode==0 && e[jNow].colType() == 0) continue;
0271 
0272           // POWHEG.
0273           if (pTdefMode == 0 || pTdefMode == 1) {
0274 
0275             // ISR - only done once as just kinematical pT.
0276             if (!FSR) {
0277               pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
0278               if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
0279 
0280             // FSR - try all outgoing partons from system before branching
0281             // as i. Note that for the hard system, there is no
0282             // "before branching" information.
0283             } else {
0284 
0285               int outSize = partonSystemsPtr->sizeOut(0);
0286               for (int iMem = 0; iMem < outSize; iMem++) {
0287                 int iNow = partonSystemsPtr->getOut(0, iMem);
0288 
0289                 // i != jNow and no carbon copies
0290                 if (iNow == jNow ) continue;
0291                 // Exclude photons (and W/Z!)
0292                 if (QEDvetoMode==0 && e[iNow].colType() == 0) continue;
0293                 if (jNow == e[iNow].daughter1()
0294                   && jNow == e[iNow].daughter2()) continue;
0295 
0296                 pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
0297                   ? false : FSR);
0298                 if (pTnow > 0.) pTemt = (pTemt < 0)
0299                   ? pTnow : min(pTemt, pTnow);
0300               }
0301              // for (iMem)
0302             }
0303             // if (!FSR)
0304           // Pythia.
0305           } else if (pTdefMode == 2) {
0306 
0307             // ISR - other incoming as recoiler.
0308             if (!FSR) {
0309               pTnow = pT(e, iInA, jNow, iInB, FSR);
0310               if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
0311               pTnow = pT(e, iInB, jNow, iInA, FSR);
0312               if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
0313 
0314             // FSR - try all final-state coloured partons as radiator
0315             //       after emission (k).
0316             } else {
0317               for (int kNow = 0; kNow < e.size(); kNow++) {
0318                 if (kNow == jNow || !e[kNow].isFinal()) continue;
0319                 if (QEDvetoMode==0 && e[kNow].colType() == 0) continue;
0320 
0321                 // For this kNow, need to have a recoiler.
0322                 // Try two incoming.
0323                 pTnow = pT(e, kNow, jNow, iInA, FSR);
0324                 if (pTnow > 0.) pTemt = (pTemt < 0)
0325                   ? pTnow : min(pTemt, pTnow);
0326                 pTnow = pT(e, kNow, jNow, iInB, FSR);
0327                 if (pTnow > 0.) pTemt = (pTemt < 0)
0328                   ? pTnow : min(pTemt, pTnow);
0329 
0330                 // Try all other outgoing.
0331                 for (int rNow = 0; rNow < e.size(); rNow++) {
0332                   if (rNow == kNow || rNow == jNow ||
0333                       !e[rNow].isFinal()) continue;
0334                   if(QEDvetoMode==0 && e[rNow].colType() == 0) continue;
0335                   pTnow = pT(e, kNow, jNow, rNow, FSR);
0336                   if (pTnow > 0.) pTemt = (pTemt < 0)
0337                     ? pTnow : min(pTemt, pTnow);
0338                 }
0339               // for (rNow)
0340               }
0341             // for (kNow)
0342             }
0343           // if (!FSR)
0344           }
0345         // if (pTdefMode)
0346         }
0347       // for (j)
0348       }
0349     }
0350     // for (xSR)
0351 
0352     return pTemt;
0353   }
0354 
0355   //--------------------------------------------------------------------------
0356 
0357   // Extraction of pThard based on the incoming event.
0358   // Assume that all the final-state particles are in a continuous block
0359   // at the end of the event and the final entry is the POWHEG emission.
0360   // If there is no POWHEG emission, then pThard is set to SCALUP.
0361 
0362   inline bool canVetoMPIStep()    { return true; }
0363   inline int  numberVetoMPIStep() { return 1; }
0364   inline bool doVetoMPIStep(int nMPI, const Event &e) {
0365     // Extra check on nMPI
0366     if (nMPI > 1) return false;
0367     int iEmt = -1;
0368     double pT1(0.), pTsum(0.);
0369 
0370     // When nFinal is set, be strict about comparing the number of final-state
0371     // particles with expectation from Born and single-real emission states.
0372     // (Note: the default from 8.309 onwards is nFinal = -1).
0373     if (nFinal > 0) {
0374       // Find if there is a POWHEG emission. Go backwards through the
0375       // event record until there is a non-final particle. Also sum pT and
0376       // find pT_1 for possible MPI vetoing
0377       int count = 0;
0378       for (int i = e.size() - 1; i > 0; i--) {
0379         if (e[i].isFinal()) {
0380           count++;
0381           pT1    = e[i].pT();
0382           pTsum += e[i].pT();
0383         } else break;
0384       }
0385       // Extra check that we have the correct final state
0386       if (count != nFinal && count != nFinal + 1) {
0387         loggerPtr->ABORT_MSG("wrong number of final state particles in event");
0388         exit(1);
0389       }
0390       // Flag if POWHEG radiation present and index
0391       isEmt = (count == nFinal) ? false : true;
0392       iEmt  = (isEmt) ? e.size() - 1 : -1;
0393 
0394     // If nFinal == -1, then go through the event and extract only the
0395     // information on the emission and its pT, but do not enforce strict
0396     // comparisons of final state multiplicity.
0397     } else {
0398 
0399       // Flag whether POWHEG radiation is present, and save index of emission.
0400       isEmt = false;
0401       for (int i = e.size() - 1; i > 0; i--) {
0402         if (e[i].isFinal()) {
0403           if ( e[i].isParton() && iEmt < 0
0404             && e[e[i].mother1()].isParton() ) {
0405             isEmt = true;
0406             iEmt = i;
0407           }
0408           pT1    = e[i].pT();
0409           pTsum += e[i].pT();
0410         } else break;
0411       }
0412     }
0413 
0414     // If there is no radiation or if pThardMode is 0 then set pThard = SCALUP.
0415     if (!isEmt || pThardMode == 0) {
0416       pThard = infoPtr->scalup();
0417 
0418     // If pThardMode is 1 then the pT of the POWHEG emission is checked against
0419     // all other incoming and outgoing partons, with the minimal value taken.
0420     } else if (pThardMode == 1) {
0421       if (nFinal < 0) {
0422         loggerPtr->WARNING_MSG(
0423           "pThard == 1 not available for nFinal == -1, reset pThard = 0.");
0424         pThardMode = 0;
0425         pThard = infoPtr->scalup();
0426       } else {
0427         pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
0428       }
0429 
0430     // If pThardMode is 2, then the pT of all final-state partons is checked
0431     // against all other incoming and outgoing partons, with the minimal value
0432     // taken.
0433     } else if (pThardMode == 2) {
0434       if (nFinal < 0) {
0435         loggerPtr->WARNING_MSG(
0436           "pThard == 2 not available for nFinal == -1, reset pThard = 0.");
0437         pThardMode = 0;
0438         pThard = infoPtr->scalup();
0439       } else {
0440         pThard = pTcalc(e, -1, -1, -1, -1, -1);
0441       }
0442     }
0443 
0444     // Find MPI veto pT if necessary.
0445     if (MPIvetoMode == 1) {
0446       pTMPI = (isEmt) ? pTsum / 2. : pT1;
0447     }
0448 
0449     // Initialise other variables.
0450     accepted   = false;
0451     nAcceptSeq = nISRveto = nFSRveto = 0;
0452 
0453     // Do not veto the event
0454     return false;
0455   }
0456 
0457   //--------------------------------------------------------------------------
0458 
0459   // ISR veto.
0460 
0461   inline bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
0462   inline bool doVetoISREmission(int, const Event &e, int iSys) {
0463     // Must be radiation from the hard system
0464     if (iSys != 0) return false;
0465 
0466     // If we already have accepted 'vetoCount' emissions in a row, do nothing
0467     if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
0468 
0469     // Pythia radiator after, emitted and recoiler after.
0470     int iRadAft = -1, iEmt = -1, iRecAft = -1;
0471     for (int i = e.size() - 1; i > 0; i--) {
0472       if (showerModel == 1) {
0473         // Pythia default.
0474         if      (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
0475         else if (iEmt    == -1 && e[i].status() ==  43) iEmt    = i;
0476         else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
0477       } else if (showerModel == 2) {
0478         // Vincia.
0479         if      (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
0480         else if (iEmt    == -1 && e[i].status() ==  43) iEmt    = i;
0481         else if (iRecAft == -1
0482           && (e[i].status() == -41 || e[i].status() == 44)) iRecAft = i;
0483       }
0484       if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
0485     }
0486     if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
0487       loggerPtr->ABORT_MSG("could not find ISR emission");
0488       exit(1);
0489     }
0490 
0491     // pTemtMode == 0: pT of emitted w.r.t. radiator
0492     // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
0493     // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
0494     int xSR      = (pTemtMode == 0) ? 0       : -1;
0495     int i        = (pTemtMode == 0) ? iRadAft : -1;
0496     int j        = (pTemtMode != 2) ? iEmt    : -1;
0497     int k        = -1;
0498     int r        = (pTemtMode == 0) ? iRecAft : -1;
0499     double pTemt = pTcalc(e, i, j, k, r, xSR);
0500 
0501     // If a Born configuration, and a photon, and QEDvetoMode=2,
0502     //  then don't veto photons, W, or Z harder than pThard.
0503     bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
0504       ? false: true;
0505 
0506     // Veto if pTemt > pThard.
0507     if (pTemt > pThard) {
0508       if(!vetoParton) {
0509         // Don't veto ANY emissions afterwards.
0510         nAcceptSeq = vetoCount-1;
0511       } else {
0512         nAcceptSeq = 0;
0513         nISRveto++;
0514         return true;
0515       }
0516     }
0517 
0518     // Else mark that an emission has been accepted and continue.
0519     nAcceptSeq++;
0520     accepted = true;
0521     return false;
0522   }
0523 
0524   //--------------------------------------------------------------------------
0525 
0526   // FSR veto.
0527 
0528   inline bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
0529   inline bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
0530     // Must be radiation from the hard system.
0531     if (iSys != 0) return false;
0532 
0533     // If we already have accepted 'vetoCount' emissions in a row, do nothing.
0534     if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
0535 
0536     // Pythia radiator (before and after), emitted and recoiler (after).
0537     int iRecAft = e.size() - 1;
0538     int iEmt    = e.size() - 2;
0539     int iRadAft = e.size() - 3;
0540     int iRadBef = e[iEmt].mother1();
0541     bool stop = false;
0542     if (showerModel == 2) {
0543       // Vincia.
0544       if ( (e[iRecAft].status() != 51 && e[iRecAft].status() != 52) ||
0545         e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop = true;
0546     } else {
0547       // Pythia default.
0548       if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
0549         e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop = true;
0550     }
0551     if (stop) {
0552       loggerPtr->ABORT_MSG("could not find FSR emission");
0553       exit(1);
0554     }
0555 
0556     // Behaviour based on pTemtMode:
0557     //  0 - pT of emitted w.r.t. radiator before
0558     //  1 - min(pT of emitted w.r.t. all incoming/outgoing)
0559     //  2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
0560     int xSR = (pTemtMode == 0) ? 1       : -1;
0561     int i   = (pTemtMode == 0) ? iRadBef : -1;
0562     int k   = (pTemtMode == 0) ? iRadAft : -1;
0563     int r   = (pTemtMode == 0) ? iRecAft : -1;
0564 
0565     // When pTemtMode is 0 or 1, iEmt has been selected.
0566     double pTemt = 0.;
0567     if (pTemtMode == 0 || pTemtMode == 1) {
0568       // Which parton is emitted, based on emittedMode:
0569       //  0 - Pythia definition of emitted
0570       //  1 - Pythia definition of radiated after emission
0571       //  2 - Random selection of emitted or radiated after emission
0572       //  3 - Try both emitted and radiated after emission
0573       int j = iRadAft;
0574       if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
0575 
0576       for (int jLoop = 0; jLoop < 2; jLoop++) {
0577         if      (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
0578         else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
0579 
0580         // For emittedMode == 3, have tried iRadAft, now try iEmt
0581         if (emittedMode != 3) break;
0582         if (k != -1) swap(j, k); else j = iEmt;
0583       }
0584 
0585     // If pTemtMode is 2, then try all final-state partons as emitted
0586     } else if (pTemtMode == 2) {
0587       pTemt = pTcalc(e, i, -1, k, r, xSR);
0588 
0589     }
0590 
0591     // If a Born configuration, and a photon, and QEDvetoMode=2,
0592     //  then don't veto photons, W's or Z's harder than pThard
0593     bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
0594       ? false: true;
0595 
0596     // Veto if pTemt > pThard
0597     if (pTemt > pThard) {
0598       if(!vetoParton) {
0599         // Don't veto ANY emissions afterwards
0600         nAcceptSeq = vetoCount-1;
0601       } else {
0602         nAcceptSeq = 0;
0603         nFSRveto++;
0604         return true;
0605       }
0606     }
0607 
0608     // Else mark that an emission has been accepted and continue
0609     nAcceptSeq++;
0610     accepted = true;
0611     return false;
0612   }
0613 
0614   //--------------------------------------------------------------------------
0615 
0616   // MPI veto.
0617 
0618   inline bool canVetoMPIEmission() {return (MPIvetoMode == 0) ? false : true;}
0619   inline bool doVetoMPIEmission(int, const Event &e) {
0620     if (MPIvetoMode == 1) {
0621       if (e[e.size() - 1].pT() > pTMPI) return true;
0622     }
0623     return false;
0624   }
0625 
0626   //--------------------------------------------------------------------------
0627 
0628   // Functions to return information.
0629 
0630   inline int getNISRveto() { return nISRveto; }
0631   inline int getNFSRveto() { return nFSRveto; }
0632 
0633   //--------------------------------------------------------------------------
0634 
0635  protected:
0636   double pThard, pTMPI;
0637 
0638  private:
0639   int    showerModel, nFinal, vetoMode, MPIvetoMode, QEDvetoMode, vetoCount;
0640   int    pThardMode, pTemtMode, emittedMode, pTdefMode;
0641   bool   accepted, isEmt;
0642   // The number of accepted emissions (in a row)
0643   // Flag for PowHeg Born or Radiation
0644   int nAcceptSeq;
0645   // Statistics on vetos
0646   unsigned long int nISRveto, nFSRveto;
0647 
0648 };
0649 
0650 //--------------------------------------------------------------------------
0651 
0652 // Declare the plugin.
0653 
0654 PYTHIA8_PLUGIN_CLASS(UserHooks, PowhegHooks, false, false, false)
0655 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
0656 
0657 //==========================================================================
0658 
0659 } // end namespace Pythia8
0660 
0661 #endif // end Pythia8_PowhegHooks_H