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