CMS 3D CMS Logo

EmissionVetoHook1.cc
Go to the documentation of this file.
1 #include "Pythia8/Pythia.h"
2 
3 using namespace Pythia8;
4 
6 
9  throw edm::Exception(edm::errors::Configuration, "Pythia8Interface") << "EmissionVeto: " << message << std::endl;
10 }
11 
12 //--------------------------------------------------------------------------
13 
14 // Routines to calculate the pT (according to pTdefMode) in a splitting:
15 // ISR: i (radiator after) -> j (emitted after) k (radiator before)
16 // FSR: i (radiator before) -> j (emitted after) k (radiator after)
17 // For the Pythia pT definition, a recoiler (after) must be specified.
18 
19 // Compute the Pythia pT separation. Based on pTLund function in History.cc
21  const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR) {
22  // Convenient shorthands for later
23  Pythia8::Vec4 radVec = e[RadAfterBranch].p();
24  Pythia8::Vec4 emtVec = e[EmtAfterBranch].p();
25  Pythia8::Vec4 recVec = e[RecAfterBranch].p();
26  int radID = e[RadAfterBranch].id();
27 
28  // Calculate virtuality of splitting
29  double sign = (FSR) ? 1. : -1.;
30  Pythia8::Vec4 Q(radVec + sign * emtVec);
31  double Qsq = sign * Q.m2Calc();
32 
33  // Mass term of radiator
34  double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? Pythia8::pow2(particleDataPtr->m0(radID)) : 0.;
35 
36  // z values for FSR and ISR
37  double z, pTnow;
38  if (FSR) {
39  // Construct 2 -> 3 variables
40  Pythia8::Vec4 sum = radVec + recVec + emtVec;
41  double m2Dip = sum.m2Calc();
42  double x1 = 2. * (sum * radVec) / m2Dip;
43  double x3 = 2. * (sum * emtVec) / m2Dip;
44  z = x1 / (x1 + x3);
45  pTnow = z * (1. - z);
46 
47  } else {
48  // Construct dipoles before/after splitting
49  Pythia8::Vec4 qBR(radVec - emtVec + recVec);
50  Pythia8::Vec4 qAR(radVec + recVec);
51  z = qBR.m2Calc() / qAR.m2Calc();
52  pTnow = (1. - z);
53  }
54 
55  // Virtuality with correct sign
56  pTnow *= (Qsq - sign * m2Rad);
57 
58  // Can get negative pT for massive splittings
59  if (pTnow < 0.) {
60  std::cout << "Warning: pTpythia was negative" << std::endl;
61  return -1.;
62  }
63 
64 #ifdef DBGOUTPUT
65  std::cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " << EmtAfterBranch << ", rec = " << RecAfterBranch
66  << ", pTnow = " << sqrt(pTnow) << std::endl;
67 #endif
68 
69  // Return pT
70  return sqrt(pTnow);
71 }
72 
73 // Compute the POWHEG pT separation between i and j
74 double EmissionVetoHook1::pTpowheg(const Pythia8::Event &e, int i, int j, bool FSR) {
75  // pT value for FSR and ISR
76  double pTnow = 0.;
77  if (FSR) {
78  // POWHEG d_ij (in CM frame). Note that the incoming beams have not
79  // been updated in the parton systems pointer yet (i.e. prior to any
80  // potential recoil).
81  int iInA = partonSystemsPtr->getInA(0);
82  int iInB = partonSystemsPtr->getInB(0);
83  double betaZ = -(e[iInA].pz() + e[iInB].pz()) / (e[iInA].e() + e[iInB].e());
84  Pythia8::Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
85  iVecBst.bst(0., 0., betaZ);
86  jVecBst.bst(0., 0., betaZ);
87  pTnow = sqrt((iVecBst + jVecBst).m2Calc() * iVecBst.e() * jVecBst.e() / Pythia8::pow2(iVecBst.e() + jVecBst.e()));
88 
89  } else {
90  // POWHEG pT_ISR is just kinematic pT
91  pTnow = e[j].pT();
92  }
93 
94  // Check result
95  if (pTnow < 0.) {
96  std::cout << "Warning: pTpowheg was negative" << std::endl;
97  return -1.;
98  }
99 
100 #ifdef DBGOUTPUT
101  std::cout << "pTpowheg: i = " << i << ", j = " << j << ", pTnow = " << pTnow << std::endl;
102 #endif
103 
104  return pTnow;
105 }
106 
107 // Calculate pT for a splitting based on pTdefMode.
108 // If j is -1, all final-state partons are tried.
109 // If i, k, r and xSR are -1, then all incoming and outgoing
110 // partons are tried.
111 // xSR set to 0 means ISR, while xSR set to 1 means FSR
112 double EmissionVetoHook1::pTcalc(const Pythia8::Event &e, int i, int j, int k, int r, int xSRin) {
113  // Loop over ISR and FSR if necessary
114  double pTemt = -1., pTnow;
115  int xSR1 = (xSRin == -1) ? 0 : xSRin;
116  int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
117  for (int xSR = xSR1; xSR < xSR2; xSR++) {
118  // FSR flag
119  bool FSR = (xSR == 0) ? false : true;
120 
121  // If all necessary arguments have been given, then directly calculate.
122  // POWHEG ISR and FSR, need i and j.
123  if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
124  pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
125 
126  // Pythia ISR, need i, j and r.
127  } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
128  pTemt = pTpythia(e, i, j, r, FSR);
129 
130  // Pythia FSR, need k, j and r.
131  } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
132  pTemt = pTpythia(e, k, j, r, FSR);
133 
134  // Otherwise need to try all possible combintations.
135  } else {
136  // Start by finding incoming legs to the hard system after
137  // branching (radiator after branching, i for ISR).
138  // Use partonSystemsPtr to find incoming just prior to the
139  // branching and track mothers.
140  int iInA = partonSystemsPtr->getInA(0);
141  int iInB = partonSystemsPtr->getInB(0);
142  while (e[iInA].mother1() != 1) {
143  iInA = e[iInA].mother1();
144  }
145  while (e[iInB].mother1() != 2) {
146  iInB = e[iInB].mother1();
147  }
148 
149  // If we do not have j, then try all final-state partons
150  int jNow = (j > 0) ? j : 0;
151  int jMax = (j > 0) ? j + 1 : e.size();
152  for (; jNow < jMax; jNow++) {
153  // Final-state only
154  if (!e[jNow].isFinal())
155  continue;
156  // Exclude photons (and W/Z!)
157  if (QEDvetoMode == 0 && e[jNow].colType() == 0)
158  continue;
159 
160  // POWHEG
161  if (pTdefMode == 0 || pTdefMode == 1) {
162  // ISR - only done once as just kinematical pT
163  if (!FSR) {
164  pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
165  if (pTnow > 0.)
166  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
167 
168  // FSR - try all outgoing partons from system before branching
169  // as i. Note that for the hard system, there is no
170  // "before branching" information.
171  } else {
172  int outSize = partonSystemsPtr->sizeOut(0);
173  for (int iMem = 0; iMem < outSize; iMem++) {
174  int iNow = partonSystemsPtr->getOut(0, iMem);
175 
176  // if i != jNow and no carbon copies
177  if (iNow == jNow)
178  continue;
179  // Exlude photons (and W/Z!)
180  if (QEDvetoMode == 0 && e[iNow].colType() == 0)
181  continue;
182  if (jNow == e[iNow].daughter1() && jNow == e[iNow].daughter2())
183  continue;
184 
185  pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) ? false : FSR);
186  if (pTnow > 0.)
187  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
188  } // for (iMem)
189 
190  } // if (!FSR)
191 
192  // Pythia
193  } else if (pTdefMode == 2) {
194  // ISR - other incoming as recoiler
195  if (!FSR) {
196  pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
197  if (pTnow > 0.)
198  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
199  pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
200  if (pTnow > 0.)
201  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
202 
203  // FSR - try all final-state coloured partons as radiator
204  // after emission (k).
205  } else {
206  for (int kNow = 0; kNow < e.size(); kNow++) {
207  if (kNow == jNow || !e[kNow].isFinal())
208  continue;
209  if (QEDvetoMode == 0 && e[kNow].colType() == 0)
210  continue;
211 
212  // For this kNow, need to have a recoiler.
213  // Try two incoming.
214  pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
215  if (pTnow > 0.)
216  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
217  pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
218  if (pTnow > 0.)
219  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
220 
221  // Try all other outgoing.
222  for (int rNow = 0; rNow < e.size(); rNow++) {
223  if (rNow == kNow || rNow == jNow || !e[rNow].isFinal())
224  continue;
225  if (QEDvetoMode == 0 && e[rNow].colType() == 0)
226  continue;
227  pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
228  if (pTnow > 0.)
229  pTemt = (pTemt < 0) ? pTnow : std::min(pTemt, pTnow);
230  } // for (rNow)
231 
232  } // for (kNow)
233  } // if (!FSR)
234  } // if (pTdefMode)
235  } // for (j)
236  }
237  } // for (xSR)
238 
239 #ifdef DBGOUTPUT
240  std::cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k << ", r = " << r << ", xSR = " << xSRin
241  << ", pTemt = " << pTemt << std::endl;
242 #endif
243 
244  return pTemt;
245 }
246 
247 //--------------------------------------------------------------------------
248 
249 // Extraction of pThard based on the incoming event.
250 // Assume that all the final-state particles are in a continuous block
251 // at the end of the event and the final entry is the POWHEG emission.
252 // If there is no POWHEG emission, then pThard is set to QRen.
253 bool EmissionVetoHook1::doVetoMPIStep(int nMPI, const Pythia8::Event &e) {
254  if (nFinalMode == 3 && pThardMode != 0)
255  fatalEmissionVeto(std::string(
256  "When nFinalMode is set to 3, ptHardMode should be set to 0, since the emission variables in doVetoMPIStep are "
257  "not set correctly case when there are three possible particle Born particle counts."));
258 
259  // Extra check on nMPI
260  if (nMPI > 1)
261  return false;
262 
263  // Find if there is a POWHEG emission. Go backwards through the
264  // event record until there is a non-final particle. Also sum pT and
265  // find pT_1 for possible MPI vetoing
266  // Flag if POWHEG radiation is present and index at the same time
267  int count = 0, inonfinal = 0;
268  double pT1 = 0., pTsum = 0.;
269  isEmt = false;
270  int iEmt = -1;
271 
272  for (int i = e.size() - 1; i > 0; i--) {
273  inonfinal = i;
274  if (e[i].isFinal()) {
275  count++;
276  pT1 = e[i].pT();
277  pTsum += e[i].pT();
278  // the following was added for bbbar4l and will be triggered by specifying nfinalmode == 2
279  // the solution provided by Tomas may not be process independent but should work for hvq and bb4l
280  // if the particle is a light quark or a gluon and
281  // comes for a light quark or a gluon (not a resonance, not a top)
282  // then it is the POWHEG emission (hvq) or the POWHEG emission in production (bb4l)
283  if (nFinalMode == 2) {
284  if ((abs(e[i].id()) < 6 || e[i].id() == 21) &&
285  (abs(e[e[i].mother1()].id()) < 6 || e[e[i].mother1()].id() == 21)) {
286  isEmt = true;
287  iEmt = i;
288  }
289  }
290  } else
291  break;
292  }
293 
294  nFinal = nFinalExt;
295  if (nFinal < 0 || nFinalMode == 1) { // nFinal is not specified from external, try to find out
296  int first = -1, myid;
297  int last = -1;
298  for (int ip = 2; ip < e.size(); ip++) {
299  myid = e[ip].id();
300  if (abs(myid) < 6 || abs(myid) == 21)
301  continue;
302  first = ip;
303  break;
304  }
305  if (first < 0)
306  fatalEmissionVeto(std::string("signal particles not found"));
307  for (int ip = first; ip < e.size(); ip++) {
308  myid = e[ip].id();
309  if (abs(myid) < 6 || abs(myid) == 21)
310  continue;
311  last = ip;
312  }
313  nFinal = last - inonfinal;
314  }
315 
316  // don't perform a cross check in case of nfinalmode == 2
317  if (nFinalMode != 2) {
318  // Extra check that we have the correct final state
319  // In POWHEG WGamma, both w+0/1 jets and w+gamma+0/1 jets are generated at the same time, which leads to three different possible numbers of particles
320  // Normally, there would be only two possible numbers of particles for X and X+1 jet
321  if (nFinalMode == 3) {
322  if (count != nFinal && count != nFinal + 1 && count != nFinal - 1)
323  fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
324  } else {
325  if (count != nFinal && count != nFinal + 1)
326  fatalEmissionVeto(std::string("Wrong number of final state particles in event"));
327  }
328  // Flag if POWHEG radiation present and index
329  if (count == nFinal + 1)
330  isEmt = true;
331  if (isEmt)
332  iEmt = e.size() - 1;
333  }
334 
335  // If there is no radiation or if pThardMode is 0 then set pThard to QRen.
336  if (!isEmt || pThardMode == 0) {
337  pThard = infoPtr->QRen();
338 
339  // If pThardMode is 1 then the pT of the POWHEG emission is checked against
340  // all other incoming and outgoing partons, with the minimal value taken
341  } else if (pThardMode == 1) {
342  pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
343 
344  // If pThardMode is 2, then the pT of all final-state partons is checked
345  // against all other incoming and outgoing partons, with the minimal value
346  // taken
347  } else if (pThardMode == 2) {
348  pThard = pTcalc(e, -1, -1, -1, -1, -1);
349  }
350 
351  // Find MPI veto pT if necessary
352  if (MPIvetoOn) {
353  pTMPI = (isEmt) ? pTsum / 2. : pT1;
354  }
355 
356  if (Verbosity)
357  std::cout << "doVetoMPIStep: QFac = " << infoPtr->QFac() << ", QRen = " << infoPtr->QRen()
358  << ", pThard = " << pThard << std::endl
359  << std::endl;
360 
361  // Initialise other variables
362  accepted = false;
363  nAcceptSeq = 0; // should not reset nISRveto = nFSRveto = nFSRvetoBB4l here
364 
365  // Do not veto the event
366  return false;
367 }
368 
369 //--------------------------------------------------------------------------
370 
371 // ISR veto
372 
373 bool EmissionVetoHook1::doVetoISREmission(int, const Pythia8::Event &e, int iSys) {
374  // Must be radiation from the hard system
375  if (iSys != 0)
376  return false;
377 
378  // If we already have accepted 'vetoCount' emissions in a row, do nothing
379  if (vetoOn && nAcceptSeq >= vetoCount)
380  return false;
381 
382  // Pythia radiator after, emitted and recoiler after.
383  int iRadAft = -1, iEmt = -1, iRecAft = -1;
384  for (int i = e.size() - 1; i > 0; i--) {
385  if (iRadAft == -1 && e[i].status() == -41)
386  iRadAft = i;
387  else if (iEmt == -1 && e[i].status() == 43)
388  iEmt = i;
389  else if (iRecAft == -1 && e[i].status() == -42)
390  iRecAft = i;
391  if (iRadAft != -1 && iEmt != -1 && iRecAft != -1)
392  break;
393  }
394  if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
395  e.list();
396  fatalEmissionVeto(std::string("Couldn't find Pythia ISR emission"));
397  }
398 
399  // pTemtMode == 0: pT of emitted w.r.t. radiator
400  // pTemtMode == 1: std::min(pT of emitted w.r.t. all incoming/outgoing)
401  // pTemtMode == 2: std::min(pT of all outgoing w.r.t. all incoming/outgoing)
402  int xSR = (pTemtMode == 0) ? 0 : -1;
403  int i = (pTemtMode == 0) ? iRadAft : -1;
404  int j = (pTemtMode != 2) ? iEmt : -1;
405  int k = -1;
406  int r = (pTemtMode == 0) ? iRecAft : -1;
407  double pTemt = pTcalc(e, i, j, k, r, xSR);
408 
409 #ifdef DBGOUTPUT
410  std::cout << "doVetoISREmission: pTemt = " << pTemt << std::endl << std::endl;
411 #endif
412 
413  // If a Born configuration, and a colorless FS, and QEDvetoMode=2,
414  // then don't veto photons, W, or Z harder than pThard
415  bool vetoParton = (!isEmt && e[iEmt].colType() == 0 && QEDvetoMode == 2) ? false : true;
416 
417  // Veto if pTemt > pThard
418  if (pTemt > pThard) {
419  if (!vetoParton) {
420  // Don't veto ANY emissions afterwards
421  nAcceptSeq = vetoCount - 1;
422  } else {
423  nAcceptSeq = 0;
424  nISRveto++;
425  return true;
426  }
427  }
428 
429  // Else mark that an emission has been accepted and continue
430  nAcceptSeq++;
431  accepted = true;
432  return false;
433 }
434 
435 //--------------------------------------------------------------------------
436 
437 // FSR veto
438 
439 bool EmissionVetoHook1::doVetoFSREmission(int, const Pythia8::Event &e, int iSys, bool inResonance) {
440  // Must be radiation from the hard system
441  if (iSys != 0)
442  return false;
443 
444  // only use for outside resonance vetos in combination with bb4l:FSREmission:veto
445  if (inResonance && settingsPtr->flag("POWHEG:bb4l:FSREmission:veto") == 1)
446  return false;
447 
448  // If we already have accepted 'vetoCount' emissions in a row, do nothing
449  if (vetoOn && nAcceptSeq >= vetoCount)
450  return false;
451 
452  // Pythia radiator (before and after), emitted and recoiler (after)
453  int iRecAft = e.size() - 1;
454  int iEmt = e.size() - 2;
455  int iRadAft = e.size() - 3;
456  int iRadBef = e[iEmt].mother1();
457  if ((e[iRecAft].status() != 52 && e[iRecAft].status() != -53) || e[iEmt].status() != 51 ||
458  e[iRadAft].status() != 51) {
459  e.list();
460  fatalEmissionVeto(std::string("Couldn't find Pythia FSR emission"));
461  }
462 
463  // Behaviour based on pTemtMode:
464  // 0 - pT of emitted w.r.t. radiator before
465  // 1 - std::min(pT of emitted w.r.t. all incoming/outgoing)
466  // 2 - std::min(pT of all outgoing w.r.t. all incoming/outgoing)
467  int xSR = (pTemtMode == 0) ? 1 : -1;
468  int i = (pTemtMode == 0) ? iRadBef : -1;
469  int k = (pTemtMode == 0) ? iRadAft : -1;
470  int r = (pTemtMode == 0) ? iRecAft : -1;
471 
472  // When pTemtMode is 0 or 1, iEmt has been selected
473  double pTemt = -1.;
474  if (pTemtMode == 0 || pTemtMode == 1) {
475  // Which parton is emitted, based on emittedMode:
476  // 0 - Pythia definition of emitted
477  // 1 - Pythia definition of radiated after emission
478  // 2 - Random selection of emitted or radiated after emission
479  // 3 - Try both emitted and radiated after emission
480  int j = iRadAft;
481  if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5))
482  j++;
483 
484  for (int jLoop = 0; jLoop < 2; jLoop++) {
485  if (jLoop == 0)
486  pTemt = pTcalc(e, i, j, k, r, xSR);
487  else if (jLoop == 1)
488  pTemt = std::min(pTemt, pTcalc(e, i, j, k, r, xSR));
489 
490  // For emittedMode == 3, have tried iRadAft, now try iEmt
491  if (emittedMode != 3)
492  break;
493  if (k != -1)
494  std::swap(j, k);
495  else
496  j = iEmt;
497  }
498 
499  // If pTemtMode is 2, then try all final-state partons as emitted
500  } else if (pTemtMode == 2) {
501  pTemt = pTcalc(e, i, -1, k, r, xSR);
502  }
503 
504 #ifdef DBGOUTPUT
505  std::cout << "doVetoFSREmission: pTemt = " << pTemt << std::endl << std::endl;
506 #endif
507 
508  // If a Born configuration, and a colorless FS, and QEDvetoMode=2,
509  // then don't veto photons, W, or Z harder than pThard
510  bool vetoParton = (!isEmt && e[iEmt].colType() == 0 && QEDvetoMode == 2) ? false : true;
511 
512  // Veto if pTemt > pThard
513  if (pTemt > pThard) {
514  if (!vetoParton) {
515  // Don't veto ANY emissions afterwards
516  nAcceptSeq = vetoCount - 1;
517  } else {
518  nAcceptSeq = 0;
519  nFSRveto++;
520  return true;
521  }
522  }
523 
524  // Else mark that an emission has been accepted and continue
525  nAcceptSeq++;
526  accepted = true;
527  return false;
528 }
529 
530 //--------------------------------------------------------------------------
531 
532 // MPI veto
533 
534 bool EmissionVetoHook1::doVetoMPIEmission(int, const Pythia8::Event &e) {
535  if (MPIvetoOn) {
536  if (e[e.size() - 1].pT() > pTMPI) {
537 #ifdef DBGOUTPUT
538  std::cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() << ", pTMPI = " << pTMPI << std::endl
539  << std::endl;
540 #endif
541  return true;
542  }
543  }
544  return false;
545 }
bool doVetoFSREmission(int, const Pythia8::Event &e, int iSys, bool) override
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:112
bool doVetoISREmission(int, const Pythia8::Event &e, int iSys) override
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ExtVec< T, 4 > Vec4
Definition: ExtVec.h:60
double pTpowheg(const Pythia8::Event &e, int i, int j, bool FSR)
bool doVetoMPIStep(int nMPI, const Pythia8::Event &e) override
double pTcalc(const Pythia8::Event &e, int i, int j, int k, int r, int xSRin)
double pTpythia(const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)
constexpr int pow2(int x)
Definition: conifer.h:28
void fatalEmissionVeto(std::string message)
bool doVetoMPIEmission(int, const Pythia8::Event &e) override