PYTHIA  8.316
PowhegHooks.h
1 // PowhegHooks.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2025 Richard Corke, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Author: Richard Corke, modified by Christian T Preuss.
7 // This class is used to perform a vetoed shower, where emissions
8 // already covered in a POWHEG NLO generator should be omitted.
9 // To first approximation the handover should happen at the SCALE
10 // of the LHA, but since the POWHEG-BOX uses a different pT definition
11 // than PYTHIA, both for ISR and FSR, a more sophisticated treatment
12 // is needed. See the online manual on POWHEG matching for details.
13 
14 #ifndef Pythia8_PowhegHooks_H
15 #define Pythia8_PowhegHooks_H
16 
17 // Includes
18 #include "Pythia8/Pythia.h"
19 #include "Pythia8/Plugins.h"
20 
21 namespace Pythia8 {
22 
23 //==========================================================================
24 
25 // Use userhooks to veto PYTHIA emissions above the POWHEG scale.
26 
27 class PowhegHooks : public UserHooks {
28 
29 public:
30 
31  // Constructors and destructor.
34  ~PowhegHooks() {}
35 
36  //--------------------------------------------------------------------------
37 
38  // Initialize settings, detailing merging strategy to use.
39  bool initAfterBeams() {
40  nFinal = settingsPtr->mode("POWHEG:nFinal");
41  vetoMode = settingsPtr->mode("POWHEG:veto");
42  vetoCount = settingsPtr->mode("POWHEG:vetoCount");
43  pThardMode = settingsPtr->mode("POWHEG:pThard");
44  pTemtMode = settingsPtr->mode("POWHEG:pTemt");
45  emittedMode = settingsPtr->mode("POWHEG:emitted");
46  pTdefMode = settingsPtr->mode("POWHEG:pTdef");
47  MPIvetoMode = settingsPtr->mode("POWHEG:MPIveto");
48  QEDvetoMode = settingsPtr->mode("POWHEG:QEDveto");
49  showerModel = settingsPtr->mode("PartonShowers:model");
50  return true;
51  }
52 
53  //--------------------------------------------------------------------------
54 
55  // Routines to calculate the pT (according to pTdefMode) in a branching:
56  // ISR: i (radiator after) -> j (emitted after) k (radiator before)
57  // FSR: i (radiator before) -> j (emitted after) k (radiator after)
58  // For the Pythia pT definitions, a recoiler (after) must be specified.
59 
60  // Top-level wrapper for shower pTs.
61  inline double pT(const Event& e, int RadAfterBranch,
62  int EmtAfterBranch, int RecAfterBranch, bool FSR) {
63  // VINCIA pT definition.
64  if (showerModel == 2)
65  return pTvincia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch);
66  // Simple-shower pT definition.
67  return pTpythia(e, RadAfterBranch, EmtAfterBranch, RecAfterBranch, FSR);
68  }
69 
70  //--------------------------------------------------------------------------
71 
72  // Compute the Pythia pT separation.
73  // Based on pTLund function in History.cc and identical to pTevol.
74  inline double pTpythia(const Event& e, int RadAfterBranch,
75  int EmtAfterBranch, int RecAfterBranch, bool FSR) {
76 
77  // Convenient shorthands for later
78  Vec4 radVec = e[RadAfterBranch].p();
79  Vec4 emtVec = e[EmtAfterBranch].p();
80  Vec4 recVec = e[RecAfterBranch].p();
81  int radID = e[RadAfterBranch].id();
82 
83  // Calculate virtuality of splitting
84  double sign = (FSR) ? 1. : -1.;
85  Vec4 Q(radVec + sign * emtVec);
86  double Qsq = sign * Q.m2Calc();
87 
88  // Mass term of radiator
89  double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
90  pow2(particleDataPtr->m0(radID)) : 0.;
91 
92  // z values for FSR and ISR
93  double z, pTnow;
94  if (FSR) {
95  // Construct 2 -> 3 variables
96  Vec4 sum = radVec + recVec + emtVec;
97  double m2Dip = sum.m2Calc();
98  double x1 = 2. * (sum * radVec) / m2Dip;
99  double x3 = 2. * (sum * emtVec) / m2Dip;
100  z = x1 / (x1 + x3);
101  pTnow = z * (1. - z);
102 
103  } else {
104  // Construct dipoles before/after splitting
105  Vec4 qBR(radVec - emtVec + recVec);
106  Vec4 qAR(radVec + recVec);
107  z = qBR.m2Calc() / qAR.m2Calc();
108  pTnow = (1. - z);
109  }
110 
111  // Virtuality with correct sign
112  pTnow *= (Qsq - sign * m2Rad);
113 
114  // Can get negative pT for massive splittings.
115  if (pTnow < 0.) {
116  loggerPtr->WARNING_MSG("negative pT");
117  return -1.;
118  }
119 
120  // Return pT
121  return sqrt(pTnow);
122  }
123 
124  //--------------------------------------------------------------------------
125 
126  // Compute the Vincia pT as in Eq. (2.63)-(2.66) in arXiv:2003.00702.
127  // Branching is assumed to be {13} {23} -> 1 3 2.
128  inline double pTvincia(const Event& event, int i1, int i3, int i2) {
129 
130  // Shorthands.
131  Vec4 p1 = event[i1].p();
132  Vec4 p3 = event[i3].p();
133  Vec4 p2 = event[i2].p();
134 
135  // Fetch mothers of 1 and 2.
136  int iMoth1 = event[i1].mother1();
137  int iMoth2 = event[i2].mother1();
138  if (iMoth1 == 0 || iMoth2 == 0) {
139  loggerPtr->ABORT_MSG("could not find mothers of particles");
140  exit(1);
141  }
142 
143  // Invariants defined as in Eq. (5) in arXiv:2008.09468.
144  double mMoth1Sq = event[iMoth1].m2();
145  double mMoth2Sq = event[iMoth2].m2();
146  double sgn1 = event[i1].isFinal() ? 1. : -1.;
147  double sgn2 = event[i2].isFinal() ? 1. : -1.;
148  double qSq13 = sgn1*(m2(sgn1*p1+p3) - mMoth1Sq);
149  double qSq23 = sgn2*(m2(sgn2*p2+p3) - mMoth2Sq);
150 
151  // Normalisation as in Eq. (6) in arXiv:2008.09468.
152  double sMax = -1.;
153  if (event[i1].isFinal() && event[i2].isFinal()) {
154  // FF.
155  sMax = m2(p1+p2+p3) - mMoth1Sq - mMoth2Sq;
156  } else if ((event[i1].isResonance() && event[i2].isFinal())
157  || (!event[i1].isFinal() && event[i2].isFinal())) {
158  // RF or IF.
159  sMax = 2.*p1*p3 + 2.*p1*p2;
160  } else if ((event[i1].isFinal() && event[i2].isResonance())
161  || (event[i1].isFinal() && !event[i2].isFinal())) {
162  // FR or FI.
163  sMax = 2.*p2*p3 + 2.*p1*p2;
164  } else if (!event[i1].isFinal() || !event[i2].isFinal()) {
165  // II.
166  sMax = 2.*p1*p2;
167  } else {
168  loggerPtr->ABORT_MSG("could not determine branching type");
169  exit(1);
170  }
171 
172  // Calculate pT2 as in Eq. (5) in arXiv:2008.09468.
173  double pT2now = qSq13*qSq23/sMax;
174 
175  // Sanity check.
176  if (pT2now < 0.) {
177  loggerPtr->WARNING_MSG("negative pT");
178  return -1.;
179  }
180 
181  // Return pT.
182  return sqrt(pT2now);
183  }
184 
185  //--------------------------------------------------------------------------
186 
187  // Compute the POWHEG pT separation between i and j.
188  inline double pTpowheg(const Event &e, int i, int j, bool FSR) {
189 
190  // pT value for FSR and ISR
191  double pTnow = 0.;
192  if (FSR) {
193  // POWHEG d_ij (in CM frame). Note that the incoming beams have not
194  // been updated in the parton systems pointer yet (i.e. prior to any
195  // potential recoil).
196  int iInA = partonSystemsPtr->getInA(0);
197  int iInB = partonSystemsPtr->getInB(0);
198  double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
199  ( e[iInA].e() + e[iInB].e() );
200  Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
201  iVecBst.bst(0., 0., betaZ);
202  jVecBst.bst(0., 0., betaZ);
203  pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
204  iVecBst.e() * jVecBst.e() /
205  pow2(iVecBst.e() + jVecBst.e()) );
206 
207  } else {
208  // POWHEG pT_ISR is just kinematic pT
209  pTnow = e[j].pT();
210  }
211 
212  // Check result.
213  if (pTnow < 0.) {
214  loggerPtr->WARNING_MSG("negative pT");
215  return -1.;
216  }
217 
218  return pTnow;
219  }
220 
221  //--------------------------------------------------------------------------
222 
223  // Calculate pT for a splitting based on pTdefMode.
224  // If j is -1, all final-state partons are tried.
225  // If i, k, r and xSR are -1, then all incoming and outgoing
226  // partons are tried.
227  // xSR set to 0 means ISR, while xSR set to 1 means FSR.
228  inline double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
229 
230  // Loop over ISR and FSR if necessary
231  double pTemt = -1., pTnow;
232  int xSR1 = (xSRin == -1) ? 0 : xSRin;
233  int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
234  for (int xSR = xSR1; xSR < xSR2; xSR++) {
235  // FSR flag
236  bool FSR = (xSR == 0) ? false : true;
237 
238  // If all necessary arguments have been given, then directly calculate.
239  // POWHEG ISR and FSR, need i and j.
240  if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
241  pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
242 
243  // Pythia ISR, need i, j and r.
244  } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
245  pTemt = pT(e, i, j, r, FSR);
246 
247  // Pythia FSR, need k, j and r.
248  } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
249  pTemt = pT(e, k, j, r, FSR);
250 
251  // Otherwise need to try all possible combinations.
252  } else {
253  // Start by finding incoming legs to the hard system after
254  // branching (radiator after branching, i for ISR).
255  // Use partonSystemsPtr to find incoming just prior to the
256  // branching and track mothers.
257  int iInA = partonSystemsPtr->getInA(0);
258  int iInB = partonSystemsPtr->getInB(0);
259  while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
260  while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
261 
262  // If we do not have j, then try all final-state partons.
263  int jNow = (j > 0) ? j : 0;
264  int jMax = (j > 0) ? j + 1 : e.size();
265  for (; jNow < jMax; jNow++) {
266 
267  // Final-state only.
268  if (!e[jNow].isFinal()) continue;
269  // Exclude photons (and W/Z!)
270  if (QEDvetoMode==0 && e[jNow].colType() == 0) continue;
271 
272  // POWHEG.
273  if (pTdefMode == 0 || pTdefMode == 1) {
274 
275  // ISR - only done once as just kinematical pT.
276  if (!FSR) {
277  pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
278  if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
279 
280  // FSR - try all outgoing partons from system before branching
281  // as i. Note that for the hard system, there is no
282  // "before branching" information.
283  } else {
284 
285  int outSize = partonSystemsPtr->sizeOut(0);
286  for (int iMem = 0; iMem < outSize; iMem++) {
287  int iNow = partonSystemsPtr->getOut(0, iMem);
288 
289  // i != jNow and no carbon copies
290  if (iNow == jNow ) continue;
291  // Exclude photons (and W/Z!)
292  if (QEDvetoMode==0 && e[iNow].colType() == 0) continue;
293  if (jNow == e[iNow].daughter1()
294  && jNow == e[iNow].daughter2()) continue;
295 
296  pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
297  ? false : FSR);
298  if (pTnow > 0.) pTemt = (pTemt < 0)
299  ? pTnow : min(pTemt, pTnow);
300  }
301  // for (iMem)
302  }
303  // if (!FSR)
304  // Pythia.
305  } else if (pTdefMode == 2) {
306 
307  // ISR - other incoming as recoiler.
308  if (!FSR) {
309  pTnow = pT(e, iInA, jNow, iInB, FSR);
310  if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
311  pTnow = pT(e, iInB, jNow, iInA, FSR);
312  if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
313 
314  // FSR - try all final-state coloured partons as radiator
315  // after emission (k).
316  } else {
317  for (int kNow = 0; kNow < e.size(); kNow++) {
318  if (kNow == jNow || !e[kNow].isFinal()) continue;
319  if (QEDvetoMode==0 && e[kNow].colType() == 0) continue;
320 
321  // For this kNow, need to have a recoiler.
322  // Try two incoming.
323  pTnow = pT(e, kNow, jNow, iInA, FSR);
324  if (pTnow > 0.) pTemt = (pTemt < 0)
325  ? pTnow : min(pTemt, pTnow);
326  pTnow = pT(e, kNow, jNow, iInB, FSR);
327  if (pTnow > 0.) pTemt = (pTemt < 0)
328  ? pTnow : min(pTemt, pTnow);
329 
330  // Try all other outgoing.
331  for (int rNow = 0; rNow < e.size(); rNow++) {
332  if (rNow == kNow || rNow == jNow ||
333  !e[rNow].isFinal()) continue;
334  if(QEDvetoMode==0 && e[rNow].colType() == 0) continue;
335  pTnow = pT(e, kNow, jNow, rNow, FSR);
336  if (pTnow > 0.) pTemt = (pTemt < 0)
337  ? pTnow : min(pTemt, pTnow);
338  }
339  // for (rNow)
340  }
341  // for (kNow)
342  }
343  // if (!FSR)
344  }
345  // if (pTdefMode)
346  }
347  // for (j)
348  }
349  }
350  // for (xSR)
351 
352  return pTemt;
353  }
354 
355  //--------------------------------------------------------------------------
356 
357  // Extraction of pThard based on the incoming event.
358  // Assume that all the final-state particles are in a continuous block
359  // at the end of the event and the final entry is the POWHEG emission.
360  // If there is no POWHEG emission, then pThard is set to SCALUP.
361 
362  inline bool canVetoMPIStep() { return true; }
363  inline int numberVetoMPIStep() { return 1; }
364  inline bool doVetoMPIStep(int nMPI, const Event &e) {
365  // Extra check on nMPI
366  if (nMPI > 1) return false;
367  int iEmt = -1;
368  double pT1(0.), pTsum(0.);
369 
370  // When nFinal is set, be strict about comparing the number of final-state
371  // particles with expectation from Born and single-real emission states.
372  // (Note: the default from 8.309 onwards is nFinal = -1).
373  if (nFinal > 0) {
374  // Find if there is a POWHEG emission. Go backwards through the
375  // event record until there is a non-final particle. Also sum pT and
376  // find pT_1 for possible MPI vetoing
377  int count = 0;
378  for (int i = e.size() - 1; i > 0; i--) {
379  if (e[i].isFinal()) {
380  count++;
381  pT1 = e[i].pT();
382  pTsum += e[i].pT();
383  } else break;
384  }
385  // Extra check that we have the correct final state
386  if (count != nFinal && count != nFinal + 1) {
387  loggerPtr->ABORT_MSG("wrong number of final state particles in event");
388  exit(1);
389  }
390  // Flag if POWHEG radiation present and index
391  isEmt = (count == nFinal) ? false : true;
392  iEmt = (isEmt) ? e.size() - 1 : -1;
393 
394  // If nFinal == -1, then go through the event and extract only the
395  // information on the emission and its pT, but do not enforce strict
396  // comparisons of final state multiplicity.
397  } else {
398 
399  // Flag whether POWHEG radiation is present, and save index of emission.
400  isEmt = false;
401  for (int i = e.size() - 1; i > 0; i--) {
402  if (e[i].isFinal()) {
403  if ( e[i].isParton() && iEmt < 0
404  && e[e[i].mother1()].isParton() ) {
405  isEmt = true;
406  iEmt = i;
407  }
408  pT1 = e[i].pT();
409  pTsum += e[i].pT();
410  } else break;
411  }
412  }
413 
414  // If there is no radiation or if pThardMode is 0 then set pThard = SCALUP.
415  if (!isEmt || pThardMode == 0) {
416  pThard = infoPtr->scalup();
417 
418  // If pThardMode is 1 then the pT of the POWHEG emission is checked against
419  // all other incoming and outgoing partons, with the minimal value taken.
420  } else if (pThardMode == 1) {
421  if (nFinal < 0) {
422  loggerPtr->WARNING_MSG(
423  "pThard == 1 not available for nFinal == -1, reset pThard = 0.");
424  pThardMode = 0;
425  pThard = infoPtr->scalup();
426  } else {
427  pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
428  }
429 
430  // If pThardMode is 2, then the pT of all final-state partons is checked
431  // against all other incoming and outgoing partons, with the minimal value
432  // taken.
433  } else if (pThardMode == 2) {
434  if (nFinal < 0) {
435  loggerPtr->WARNING_MSG(
436  "pThard == 2 not available for nFinal == -1, reset pThard = 0.");
437  pThardMode = 0;
438  pThard = infoPtr->scalup();
439  } else {
440  pThard = pTcalc(e, -1, -1, -1, -1, -1);
441  }
442  }
443 
444  // Find MPI veto pT if necessary.
445  if (MPIvetoMode == 1) {
446  pTMPI = (isEmt) ? pTsum / 2. : pT1;
447  }
448 
449  // Initialise other variables.
450  accepted = false;
451  nAcceptSeq = nISRveto = nFSRveto = 0;
452 
453  // Do not veto the event
454  return false;
455  }
456 
457  //--------------------------------------------------------------------------
458 
459  // ISR veto.
460 
461  inline bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
462  inline bool doVetoISREmission(int, const Event &e, int iSys) {
463  // Must be radiation from the hard system
464  if (iSys != 0) return false;
465 
466  // If we already have accepted 'vetoCount' emissions in a row, do nothing
467  if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
468 
469  // Pythia radiator after, emitted and recoiler after.
470  int iRadAft = -1, iEmt = -1, iRecAft = -1;
471  for (int i = e.size() - 1; i > 0; i--) {
472  if (showerModel == 1) {
473  // Pythia default.
474  if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
475  else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
476  else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
477  } else if (showerModel == 2) {
478  // Vincia.
479  if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
480  else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
481  else if (iRecAft == -1
482  && (e[i].status() == -41 || e[i].status() == 44)) iRecAft = i;
483  }
484  if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
485  }
486  if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
487  loggerPtr->ABORT_MSG("could not find ISR emission");
488  exit(1);
489  }
490 
491  // pTemtMode == 0: pT of emitted w.r.t. radiator
492  // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
493  // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
494  int xSR = (pTemtMode == 0) ? 0 : -1;
495  int i = (pTemtMode == 0) ? iRadAft : -1;
496  int j = (pTemtMode != 2) ? iEmt : -1;
497  int k = -1;
498  int r = (pTemtMode == 0) ? iRecAft : -1;
499  double pTemt = pTcalc(e, i, j, k, r, xSR);
500 
501  // If a Born configuration, and a photon, and QEDvetoMode=2,
502  // then don't veto photons, W, or Z harder than pThard.
503  bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
504  ? false: true;
505 
506  // Veto if pTemt > pThard.
507  if (pTemt > pThard) {
508  if(!vetoParton) {
509  // Don't veto ANY emissions afterwards.
510  nAcceptSeq = vetoCount-1;
511  } else {
512  nAcceptSeq = 0;
513  nISRveto++;
514  return true;
515  }
516  }
517 
518  // Else mark that an emission has been accepted and continue.
519  nAcceptSeq++;
520  accepted = true;
521  return false;
522  }
523 
524  //--------------------------------------------------------------------------
525 
526  // FSR veto.
527 
528  inline bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
529  inline bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
530  // Must be radiation from the hard system.
531  if (iSys != 0) return false;
532 
533  // If we already have accepted 'vetoCount' emissions in a row, do nothing.
534  if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
535 
536  // Pythia radiator (before and after), emitted and recoiler (after).
537  int iRecAft = e.size() - 1;
538  int iEmt = e.size() - 2;
539  int iRadAft = e.size() - 3;
540  int iRadBef = e[iEmt].mother1();
541  bool stop = false;
542  if (showerModel == 2) {
543  // Vincia.
544  if ( (e[iRecAft].status() != 51 && e[iRecAft].status() != 52) ||
545  e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop = true;
546  } else {
547  // Pythia default.
548  if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
549  e[iEmt].status() != 51 || e[iRadAft].status() != 51) stop = true;
550  }
551  if (stop) {
552  loggerPtr->ABORT_MSG("could not find FSR emission");
553  exit(1);
554  }
555 
556  // Behaviour based on pTemtMode:
557  // 0 - pT of emitted w.r.t. radiator before
558  // 1 - min(pT of emitted w.r.t. all incoming/outgoing)
559  // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
560  int xSR = (pTemtMode == 0) ? 1 : -1;
561  int i = (pTemtMode == 0) ? iRadBef : -1;
562  int k = (pTemtMode == 0) ? iRadAft : -1;
563  int r = (pTemtMode == 0) ? iRecAft : -1;
564 
565  // When pTemtMode is 0 or 1, iEmt has been selected.
566  double pTemt = 0.;
567  if (pTemtMode == 0 || pTemtMode == 1) {
568  // Which parton is emitted, based on emittedMode:
569  // 0 - Pythia definition of emitted
570  // 1 - Pythia definition of radiated after emission
571  // 2 - Random selection of emitted or radiated after emission
572  // 3 - Try both emitted and radiated after emission
573  int j = iRadAft;
574  if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
575 
576  for (int jLoop = 0; jLoop < 2; jLoop++) {
577  if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
578  else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
579 
580  // For emittedMode == 3, have tried iRadAft, now try iEmt
581  if (emittedMode != 3) break;
582  if (k != -1) swap(j, k); else j = iEmt;
583  }
584 
585  // If pTemtMode is 2, then try all final-state partons as emitted
586  } else if (pTemtMode == 2) {
587  pTemt = pTcalc(e, i, -1, k, r, xSR);
588 
589  }
590 
591  // If a Born configuration, and a photon, and QEDvetoMode=2,
592  // then don't veto photons, W's or Z's harder than pThard
593  bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
594  ? false: true;
595 
596  // Veto if pTemt > pThard
597  if (pTemt > pThard) {
598  if(!vetoParton) {
599  // Don't veto ANY emissions afterwards
600  nAcceptSeq = vetoCount-1;
601  } else {
602  nAcceptSeq = 0;
603  nFSRveto++;
604  return true;
605  }
606  }
607 
608  // Else mark that an emission has been accepted and continue
609  nAcceptSeq++;
610  accepted = true;
611  return false;
612  }
613 
614  //--------------------------------------------------------------------------
615 
616  // MPI veto.
617 
618  inline bool canVetoMPIEmission() {return (MPIvetoMode == 0) ? false : true;}
619  inline bool doVetoMPIEmission(int, const Event &e) {
620  if (MPIvetoMode == 1) {
621  if (e[e.size() - 1].pT() > pTMPI) return true;
622  }
623  return false;
624  }
625 
626  //--------------------------------------------------------------------------
627 
628  // Functions to return information.
629 
630  inline int getNISRveto() { return nISRveto; }
631  inline int getNFSRveto() { return nFSRveto; }
632 
633  //--------------------------------------------------------------------------
634 
635  protected:
636  double pThard, pTMPI;
637 
638  private:
639  int showerModel, nFinal, vetoMode, MPIvetoMode, QEDvetoMode, vetoCount;
640  int pThardMode, pTemtMode, emittedMode, pTdefMode;
641  bool accepted, isEmt;
642  // The number of accepted emissions (in a row)
643  // Flag for PowHeg Born or Radiation
644  int nAcceptSeq;
645  // Statistics on vetos
646  unsigned long int nISRveto, nFSRveto;
647 
648 };
649 
650 //--------------------------------------------------------------------------
651 
652 // Declare the plugin.
653 
654 PYTHIA8_PLUGIN_CLASS(UserHooks, PowhegHooks, false, false, false)
655 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
656 
657 //==========================================================================
658 
659 } // end namespace Pythia8
660 
661 #endif // end Pythia8_PowhegHooks_H
double pTpowheg(const Event &e, int i, int j, bool FSR)
Compute the POWHEG pT separation between i and j.
Definition: PowhegHooks.h:188
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:173
bool initAfterBeams()
Initialize settings, detailing merging strategy to use.
Definition: PowhegHooks.h:39
Settings * settingsPtr
Pointer to the settings database.
Definition: PhysicsBase.h:85
The Event class holds all info on the generated event.
Definition: Event.h:408
bool canVetoMPIStep()
Definition: PowhegHooks.h:362
PartonSystems * partonSystemsPtr
Pointer to information on subcollision parton locations.
Definition: PhysicsBase.h:116
Rndm * rndmPtr
Pointer to the random number generator.
Definition: PhysicsBase.h:97
Definition: Logger.h:23
Use userhooks to veto PYTHIA emissions above the POWHEG scale.
Definition: PowhegHooks.h:27
ParticleData * particleDataPtr
Pointer to the particle data table.
Definition: PhysicsBase.h:88
bool canVetoMPIEmission()
MPI veto.
Definition: PowhegHooks.h:618
double pTvincia(const Event &event, int i1, int i3, int i2)
Definition: PowhegHooks.h:128
bool doVetoMPIStep(int nMPI, const Event &e)
Definition: PowhegHooks.h:364
int getNISRveto()
Functions to return information.
Definition: PowhegHooks.h:630
void bst(double betaX, double betaY, double betaZ)
Boost (simple).
Definition: Basics.cc:433
bool doVetoFSREmission(int, const Event &e, int iSys, bool)
Definition: PowhegHooks.h:529
bool doVetoISREmission(int, const Event &e, int iSys)
Definition: PowhegHooks.h:462
int size() const
Event record size.
Definition: Event.h:459
bool canVetoISREmission()
ISR veto.
Definition: PowhegHooks.h:461
double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)
Definition: PowhegHooks.h:74
double m2(const Vec4 &v1)
The squared invariant mass of one or more four-vectors.
Definition: Basics.cc:597
Logger * loggerPtr
Pointer to logger.
Definition: PhysicsBase.h:91
bool doVetoMPIEmission(int, const Event &e)
Definition: PowhegHooks.h:619
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:72
Info * infoPtr
Definition: PhysicsBase.h:82
UserHooks is base class for user access to program execution.
Definition: UserHooks.h:32
int numberVetoMPIStep()
Up to how many MPI steps should be checked.
Definition: PowhegHooks.h:363
double flat()
Generate next random number uniformly between 0 and 1.
Definition: Basics.cc:189
double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin)
Definition: PowhegHooks.h:228
bool canVetoFSREmission()
FSR veto.
Definition: PowhegHooks.h:528
double pT(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)
Top-level wrapper for shower pTs.
Definition: PowhegHooks.h:61
Definition: Basics.h:32
PowhegHooks()
Constructors and destructor.
Definition: PowhegHooks.h:32
Definition: Settings.h:196