StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
PowhegHooks.h
1 // PowhegHooks.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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.
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 merging for details.
13 
14 #ifndef Pythia8_PowhegHooks_H
15 #define Pythia8_PowhegHooks_H
16 
17 // Includes
18 #include "Pythia8/Pythia.h"
19 
20 namespace Pythia8 {
21 
22 //==========================================================================
23 
24 // Use userhooks to veto PYTHIA emissions above the POWHEG scale.
25 
26 class PowhegHooks : public UserHooks {
27 
28 public:
29 
30  // Constructor and destructor.
31  PowhegHooks() {}
32  ~PowhegHooks() {}
33 
34 //--------------------------------------------------------------------------
35 
36  // Initialize settings, detailing merging strategy to use.
37  bool initAfterBeams() {
38  nFinal = settingsPtr->mode("POWHEG:nFinal");
39  vetoMode = settingsPtr->mode("POWHEG:veto");
40  vetoCount = settingsPtr->mode("POWHEG:vetoCount");
41  pThardMode = settingsPtr->mode("POWHEG:pThard");
42  pTemtMode = settingsPtr->mode("POWHEG:pTemt");
43  emittedMode = settingsPtr->mode("POWHEG:emitted");
44  pTdefMode = settingsPtr->mode("POWHEG:pTdef");
45  MPIvetoMode = settingsPtr->mode("POWHEG:MPIveto");
46  QEDvetoMode = settingsPtr->mode("POWHEG:QEDveto");
47  return true;
48  }
49 
50 //--------------------------------------------------------------------------
51 
52  // Routines to calculate the pT (according to pTdefMode) in a splitting:
53  // ISR: i (radiator after) -> j (emitted after) k (radiator before)
54  // FSR: i (radiator before) -> j (emitted after) k (radiator after)
55  // For the Pythia pT definition, a recoiler (after) must be specified.
56 
57  // Compute the Pythia pT separation. Based on pTLund function in History.cc
58  inline double pTpythia(const Event &e, int RadAfterBranch,
59  int EmtAfterBranch, int RecAfterBranch, bool FSR) {
60 
61  // Convenient shorthands for later
62  Vec4 radVec = e[RadAfterBranch].p();
63  Vec4 emtVec = e[EmtAfterBranch].p();
64  Vec4 recVec = e[RecAfterBranch].p();
65  int radID = e[RadAfterBranch].id();
66 
67  // Calculate virtuality of splitting
68  double sign = (FSR) ? 1. : -1.;
69  Vec4 Q(radVec + sign * emtVec);
70  double Qsq = sign * Q.m2Calc();
71 
72  // Mass term of radiator
73  double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
74  pow2(particleDataPtr->m0(radID)) : 0.;
75 
76  // z values for FSR and ISR
77  double z, pTnow;
78  if (FSR) {
79  // Construct 2 -> 3 variables
80  Vec4 sum = radVec + recVec + emtVec;
81  double m2Dip = sum.m2Calc();
82  double x1 = 2. * (sum * radVec) / m2Dip;
83  double x3 = 2. * (sum * emtVec) / m2Dip;
84  z = x1 / (x1 + x3);
85  pTnow = z * (1. - z);
86 
87  } else {
88  // Construct dipoles before/after splitting
89  Vec4 qBR(radVec - emtVec + recVec);
90  Vec4 qAR(radVec + recVec);
91  z = qBR.m2Calc() / qAR.m2Calc();
92  pTnow = (1. - z);
93  }
94 
95  // Virtuality with correct sign
96  pTnow *= (Qsq - sign * m2Rad);
97 
98  // Can get negative pT for massive splittings
99  if (pTnow < 0.) {
100  cout << "Warning: pTpythia was negative" << endl;
101  return -1.;
102  }
103 
104  // Return pT
105  return sqrt(pTnow);
106  }
107 
108  // Compute the POWHEG pT separation between i and j
109  inline double pTpowheg(const Event &e, int i, int j, bool FSR) {
110 
111  // pT value for FSR and ISR
112  double pTnow = 0.;
113  if (FSR) {
114  // POWHEG d_ij (in CM frame). Note that the incoming beams have not
115  // been updated in the parton systems pointer yet (i.e. prior to any
116  // potential recoil).
117  int iInA = partonSystemsPtr->getInA(0);
118  int iInB = partonSystemsPtr->getInB(0);
119  double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
120  ( e[iInA].e() + e[iInB].e() );
121  Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
122  iVecBst.bst(0., 0., betaZ);
123  jVecBst.bst(0., 0., betaZ);
124  pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
125  iVecBst.e() * jVecBst.e() /
126  pow2(iVecBst.e() + jVecBst.e()) );
127 
128  } else {
129  // POWHEG pT_ISR is just kinematic pT
130  pTnow = e[j].pT();
131  }
132 
133  // Check result
134  if (pTnow < 0.) {
135  cout << "Warning: pTpowheg was negative" << endl;
136  return -1.;
137  }
138 
139  return pTnow;
140  }
141 
142  // Calculate pT for a splitting based on pTdefMode.
143  // If j is -1, all final-state partons are tried.
144  // If i, k, r and xSR are -1, then all incoming and outgoing
145  // partons are tried.
146  // xSR set to 0 means ISR, while xSR set to 1 means FSR
147  inline double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) {
148 
149  // Loop over ISR and FSR if necessary
150  double pTemt = -1., pTnow;
151  int xSR1 = (xSRin == -1) ? 0 : xSRin;
152  int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
153  for (int xSR = xSR1; xSR < xSR2; xSR++) {
154  // FSR flag
155  bool FSR = (xSR == 0) ? false : true;
156 
157  // If all necessary arguments have been given, then directly calculate.
158  // POWHEG ISR and FSR, need i and j.
159  if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
160  pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
161 
162  // Pythia ISR, need i, j and r.
163  } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
164  pTemt = pTpythia(e, i, j, r, FSR);
165 
166  // Pythia FSR, need k, j and r.
167  } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
168  pTemt = pTpythia(e, k, j, r, FSR);
169 
170  // Otherwise need to try all possible combinations.
171  } else {
172  // Start by finding incoming legs to the hard system after
173  // branching (radiator after branching, i for ISR).
174  // Use partonSystemsPtr to find incoming just prior to the
175  // branching and track mothers.
176  int iInA = partonSystemsPtr->getInA(0);
177  int iInB = partonSystemsPtr->getInB(0);
178  while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
179  while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
180 
181  // If we do not have j, then try all final-state partons
182  int jNow = (j > 0) ? j : 0;
183  int jMax = (j > 0) ? j + 1 : e.size();
184  for (; jNow < jMax; jNow++) {
185 
186  // Final-state only
187  if (!e[jNow].isFinal()) continue;
188  // Exclude photons (and W/Z!)
189  if (QEDvetoMode==0 && e[jNow].colType() == 0) continue;
190 
191  // POWHEG
192  if (pTdefMode == 0 || pTdefMode == 1) {
193 
194  // ISR - only done once as just kinematical pT
195  if (!FSR) {
196  pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
197  if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
198 
199  // FSR - try all outgoing partons from system before branching
200  // as i. Note that for the hard system, there is no
201  // "before branching" information.
202  } else {
203 
204  int outSize = partonSystemsPtr->sizeOut(0);
205  for (int iMem = 0; iMem < outSize; iMem++) {
206  int iNow = partonSystemsPtr->getOut(0, iMem);
207 
208  // i != jNow and no carbon copies
209  if (iNow == jNow ) continue;
210  // Exclude photons (and W/Z!)
211  if (QEDvetoMode==0 && e[iNow].colType() == 0) continue;
212  if (jNow == e[iNow].daughter1()
213  && jNow == e[iNow].daughter2()) continue;
214 
215  pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0)
216  ? false : FSR);
217  if (pTnow > 0.) pTemt = (pTemt < 0)
218  ? pTnow : min(pTemt, pTnow);
219  }
220  // for (iMem)
221  }
222  // if (!FSR)
223  // Pythia
224  } else if (pTdefMode == 2) {
225 
226  // ISR - other incoming as recoiler
227  if (!FSR) {
228  pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
229  if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
230  pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
231  if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
232 
233  // FSR - try all final-state coloured partons as radiator
234  // after emission (k).
235  } else {
236  for (int kNow = 0; kNow < e.size(); kNow++) {
237  if (kNow == jNow || !e[kNow].isFinal()) continue;
238  if (QEDvetoMode==0 && e[kNow].colType() == 0) continue;
239 
240  // For this kNow, need to have a recoiler.
241  // Try two incoming.
242  pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
243  if (pTnow > 0.) pTemt = (pTemt < 0)
244  ? pTnow : min(pTemt, pTnow);
245  pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
246  if (pTnow > 0.) pTemt = (pTemt < 0)
247  ? pTnow : min(pTemt, pTnow);
248 
249  // Try all other outgoing.
250  for (int rNow = 0; rNow < e.size(); rNow++) {
251  if (rNow == kNow || rNow == jNow ||
252  !e[rNow].isFinal()) continue;
253  if(QEDvetoMode==0 && e[rNow].colType() == 0) continue;
254  pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
255  if (pTnow > 0.) pTemt = (pTemt < 0)
256  ? pTnow : min(pTemt, pTnow);
257  }
258  // for (rNow)
259  }
260  // for (kNow)
261  }
262  // if (!FSR)
263  }
264  // if (pTdefMode)
265  }
266  // for (j)
267  }
268  }
269  // for (xSR)
270 
271  return pTemt;
272  }
273 
274 //--------------------------------------------------------------------------
275 
276  // Extraction of pThard based on the incoming event.
277  // Assume that all the final-state particles are in a continuous block
278  // at the end of the event and the final entry is the POWHEG emission.
279  // If there is no POWHEG emission, then pThard is set to SCALUP.
280 
281  inline bool canVetoMPIStep() { return true; }
282  inline int numberVetoMPIStep() { return 1; }
283  inline bool doVetoMPIStep(int nMPI, const Event &e) {
284  // Extra check on nMPI
285  if (nMPI > 1) return false;
286 
287  // Find if there is a POWHEG emission. Go backwards through the
288  // event record until there is a non-final particle. Also sum pT and
289  // find pT_1 for possible MPI vetoing
290  int count = 0;
291  double pT1 = 0., pTsum = 0.;
292  for (int i = e.size() - 1; i > 0; i--) {
293  if (e[i].isFinal()) {
294  count++;
295  pT1 = e[i].pT();
296  pTsum += e[i].pT();
297  } else break;
298  }
299  // Extra check that we have the correct final state
300  if (count != nFinal && count != nFinal + 1) {
301  cout << "Error: wrong number of final state particles in event" << endl;
302  exit(1);
303  }
304  // Flag if POWHEG radiation present and index
305  isEmt = (count == nFinal) ? false : true;
306  int iEmt = (isEmt) ? e.size() - 1 : -1;
307 
308  // If there is no radiation or if pThardMode is 0 then set pThard = SCALUP.
309  if (!isEmt || pThardMode == 0) {
310  pThard = infoPtr->scalup();
311 
312  // If pThardMode is 1 then the pT of the POWHEG emission is checked against
313  // all other incoming and outgoing partons, with the minimal value taken
314  } else if (pThardMode == 1) {
315  pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
316 
317  // If pThardMode is 2, then the pT of all final-state partons is checked
318  // against all other incoming and outgoing partons, with the minimal value
319  // taken
320  } else if (pThardMode == 2) {
321  pThard = pTcalc(e, -1, -1, -1, -1, -1);
322  }
323 
324  // Find MPI veto pT if necessary
325  if (MPIvetoMode == 1) {
326  pTMPI = (isEmt) ? pTsum / 2. : pT1;
327  }
328 
329  // Initialise other variables
330  accepted = false;
331  nAcceptSeq = nISRveto = nFSRveto = 0;
332 
333  // Do not veto the event
334  return false;
335  }
336 
337 //--------------------------------------------------------------------------
338 
339  // ISR veto
340 
341  inline bool canVetoISREmission() { return (vetoMode == 0) ? false : true; }
342  inline bool doVetoISREmission(int, const Event &e, int iSys) {
343  // Must be radiation from the hard system
344  if (iSys != 0) return false;
345 
346  // If we already have accepted 'vetoCount' emissions in a row, do nothing
347  if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
348 
349  // Pythia radiator after, emitted and recoiler after.
350  int iRadAft = -1, iEmt = -1, iRecAft = -1;
351  for (int i = e.size() - 1; i > 0; i--) {
352  if (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
353  else if (iEmt == -1 && e[i].status() == 43) iEmt = i;
354  else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
355  if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
356  }
357  if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
358  e.list();
359  cout << "Error: couldn't find Pythia ISR emission" << endl;
360  exit(1);
361  }
362 
363  // pTemtMode == 0: pT of emitted w.r.t. radiator
364  // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
365  // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
366  int xSR = (pTemtMode == 0) ? 0 : -1;
367  int i = (pTemtMode == 0) ? iRadAft : -1;
368  int j = (pTemtMode != 2) ? iEmt : -1;
369  int k = -1;
370  int r = (pTemtMode == 0) ? iRecAft : -1;
371  double pTemt = pTcalc(e, i, j, k, r, xSR);
372 
373  // If a Born configuration, and a photon, and QEDvetoMode=2,
374  // then don't veto photons, W, or Z harder than pThard
375  bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
376  ? false: true;
377 
378  // Veto if pTemt > pThard
379  if (pTemt > pThard) {
380  if(!vetoParton) {
381  // Don't veto ANY emissions afterwards
382  nAcceptSeq = vetoCount-1;
383  } else {
384  nAcceptSeq = 0;
385  nISRveto++;
386  return true;
387  }
388  }
389 
390  // Else mark that an emission has been accepted and continue
391  nAcceptSeq++;
392  accepted = true;
393  return false;
394  }
395 
396 //--------------------------------------------------------------------------
397 
398  // FSR veto
399 
400  inline bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; }
401  inline bool doVetoFSREmission(int, const Event &e, int iSys, bool) {
402  // Must be radiation from the hard system
403  if (iSys != 0) return false;
404 
405  // If we already have accepted 'vetoCount' emissions in a row, do nothing
406  if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false;
407 
408  // Pythia radiator (before and after), emitted and recoiler (after)
409  int iRecAft = e.size() - 1;
410  int iEmt = e.size() - 2;
411  int iRadAft = e.size() - 3;
412  int iRadBef = e[iEmt].mother1();
413  if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
414  e[iEmt].status() != 51 || e[iRadAft].status() != 51) {
415  e.list();
416  cout << "Error: couldn't find Pythia FSR emission" << endl;
417  exit(1);
418  }
419 
420  // Behaviour based on pTemtMode:
421  // 0 - pT of emitted w.r.t. radiator before
422  // 1 - min(pT of emitted w.r.t. all incoming/outgoing)
423  // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
424  int xSR = (pTemtMode == 0) ? 1 : -1;
425  int i = (pTemtMode == 0) ? iRadBef : -1;
426  int k = (pTemtMode == 0) ? iRadAft : -1;
427  int r = (pTemtMode == 0) ? iRecAft : -1;
428 
429  // When pTemtMode is 0 or 1, iEmt has been selected
430  double pTemt = 0.;
431  if (pTemtMode == 0 || pTemtMode == 1) {
432  // Which parton is emitted, based on emittedMode:
433  // 0 - Pythia definition of emitted
434  // 1 - Pythia definition of radiated after emission
435  // 2 - Random selection of emitted or radiated after emission
436  // 3 - Try both emitted and radiated after emission
437  int j = iRadAft;
438  if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
439 
440  for (int jLoop = 0; jLoop < 2; jLoop++) {
441  if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
442  else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
443 
444  // For emittedMode == 3, have tried iRadAft, now try iEmt
445  if (emittedMode != 3) break;
446  if (k != -1) swap(j, k); else j = iEmt;
447  }
448 
449  // If pTemtMode is 2, then try all final-state partons as emitted
450  } else if (pTemtMode == 2) {
451  pTemt = pTcalc(e, i, -1, k, r, xSR);
452 
453  }
454 
455  // If a Born configuration, and a photon, and QEDvetoMode=2,
456  // then don't veto photons, W's or Z's harder than pThard
457  bool vetoParton = (!isEmt && e[iEmt].colType()==0 && QEDvetoMode==2)
458  ? false: true;
459 
460  // Veto if pTemt > pThard
461  if (pTemt > pThard) {
462  if(!vetoParton) {
463  // Don't veto ANY emissions afterwards
464  nAcceptSeq = vetoCount-1;
465  } else {
466  nAcceptSeq = 0;
467  nFSRveto++;
468  return true;
469  }
470  }
471 
472  // Else mark that an emission has been accepted and continue
473  nAcceptSeq++;
474  accepted = true;
475  return false;
476  }
477 
478 //--------------------------------------------------------------------------
479 
480  // MPI veto
481 
482  inline bool canVetoMPIEmission() {return (MPIvetoMode == 0) ? false : true;}
483  inline bool doVetoMPIEmission(int, const Event &e) {
484  if (MPIvetoMode == 1) {
485  if (e[e.size() - 1].pT() > pTMPI) return true;
486  }
487  return false;
488  }
489 
490 //--------------------------------------------------------------------------
491 
492  // Functions to return information
493 
494  inline int getNISRveto() { return nISRveto; }
495  inline int getNFSRveto() { return nFSRveto; }
496 
497 //--------------------------------------------------------------------------
498 
499 private:
500  int nFinal, vetoMode, vetoCount, pThardMode, pTemtMode,
501  emittedMode, pTdefMode, MPIvetoMode, QEDvetoMode;
502  double pThard, pTMPI;
503  bool accepted, isEmt;
504  // The number of accepted emissions (in a row)
505  // Flag for PowHeg Born or Radiation
506  int nAcceptSeq;
507  // Statistics on vetos
508  unsigned long int nISRveto, nFSRveto;
509 
510 };
511 
512 //==========================================================================
513 
514 } // end namespace Pythia8
515 
516 #endif // end Pythia8_PowhegHooks_H
Definition: AgUStep.h:26