CMS 3D CMS Logo

PowhegHooksBB4L.h
Go to the documentation of this file.
1 // PowhegHooksBB4L.h
2 // Copyright (C) 2017 Silvia Ferrario Ravasio, Tomas Jezo, Paolo Nason, Markus Seidel
3 // inspired by PowhegHooks.h by Richard Corke
4 // adjusted to work with EmissionVetoHook1 in CMSSW by Alexander Grohsjean
5 
6 #ifndef Pythia8_PowhegHooksBB4L_H
7 #define Pythia8_PowhegHooksBB4L_H
8 
9 // Includes
10 #include "Pythia8/Pythia.h"
11 #include <cassert>
12 struct{
13  int radtype;
14 } radtype_;
15 
16 namespace Pythia8 {
17 
18  class PowhegHooksBB4L : public UserHooks {
19 
20  public:
21 
22  //--- Constructor and destructor -------------------------------------------
24  ~PowhegHooksBB4L() override {
25  std::cout << "Number of FSR vetoed in BB4l = " << nFSRvetoBB4l << std::endl;
26  }
27 
28  //--- Initialization -----------------------------------------------------------------------
29  bool initAfterBeams() override {
30  // settings of this class
31  vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto");
32  onlyDistance1 = settingsPtr->flag("POWHEG:bb4l:FSREmission:onlyDistance1");
33  dryRunFSR = settingsPtr->flag("POWHEG:bb4l:FSREmission:dryRun");
34  vetoAtPL = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoAtPL");
35  vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED");
36  vetoPartonLevel = settingsPtr->flag("POWHEG:bb4l:PartonLevel:veto");
37  excludeFSRConflicting = settingsPtr->flag("POWHEG:bb4l:PartonLevel:excludeFSRConflicting");
38  debug = settingsPtr->flag("POWHEG:bb4l:DEBUG");
39  scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto");
40  vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame");
41  pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto");
42  //vetoProduction = (settingsPtr->mode("POWHEG:veto")==1);
43  pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto");
44  return true;
45  }
46 
47 
48  //--- PROCESS LEVEL HOOK ---------------------------------------------------
49 
50  // called at the LHE level
51  inline bool canVetoProcessLevel() override { return true; }
52  inline bool doVetoProcessLevel(Event &e) override {
53 
54  // extract the radtype from the event comment
55  stringstream ss;
56  // use eventattribute as comments not filled when using edm input
57  //ss << infoPtr->getEventComments();
58  ss << infoPtr->getEventAttribute("#rwgt");
59  string temp;
60  ss >> temp >> radtype_.radtype;
61  assert (temp == "#rwgt");
62 
63  // find last top and the last anti-top in the record
64  int i_top = -1, i_atop = -1;
65  for (int i = 0; i < e.size(); i++) {
66  if (e[i].id() == 6) i_top = i;
67  if (e[i].id() == -6) i_atop = i;
68  }
69  if (i_top != -1)
70  topresscale = findresscale(i_top, e);
71  else
72  topresscale = 1e30;
73  if (i_top != -1)
74  atopresscale = findresscale(i_atop, e);
75  else
76  atopresscale = 1e30;
77  // initialize stuff
78  doVetoFSRInit();
79  // do not veto, ever
80  return false;
81  }
82 
83  //--- PARTON LEVEL HOOK ----------------------------------------------------
84 
85  // called after shower
86  bool retryPartonLevel() override { return vetoPartonLevel || vetoAtPL; }
87  inline bool canVetoPartonLevel() override { return vetoPartonLevel || vetoAtPL; }
88  inline bool doVetoPartonLevel(const Event &e) override {
89  if(radtype_.radtype==2)
90  return false;
91  if (debug){
92  if (dryRunFSR && wouldVetoFsr) {
93  double scale = getdechardness(vetoTopCharge, e);
94  cout << "FSRdecScale = " << vetoDecScale << ", PLdecScale = " << scale << ", ratio = " << vetoDecScale/scale << endl;
95  }
96  }
97  if (vetoPartonLevel) {
98  double topdecscale = getdechardness(1, e);
99  double atopdecscale = getdechardness(-1, e);
100  if ((topdecscale > topresscale) || (atopdecscale > atopresscale)) {
101  //if(dryRunFSR && ! wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0;
102  return true;
103  }
104  else
105  //if(dryRunFSR && wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0;
106  return false;
107  }
108  if (vetoAtPL) {
109  if (dryRunFSR && wouldVetoFsr) return true;
110  else return false;
111  }
112  return false;
113  }
114 
115  //--- FSR EMISSION LEVEL HOOK ----------------------------------------------
116 
117  // FSR veto: this should be true if we want PowhegHooksBB4l veto in decay
118  // OR PowhegHooks veto in production. (The virtual method
119  // PowhegHooks::canVetoFSREmission has been replaced by
120  // PowhegHooksBB4L::canVetoFSREmission, so FSR veto in production
121  // must be handled here. ISR and MPI veto are instead still
122  // handled by PowhegHooks.)
123  inline bool canVetoFSREmission() override { return vetoFSREmission; } // || vetoProduction; }
124  inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) override {
126  //VETO INSIDE THE RESONANCE //
128  if (inResonance && vetoFSREmission) {
129  int iRecAft = e.size() - 1;
130  int iEmt = e.size() - 2;
131  int iRadAft = e.size() - 3;
132  int iRadBef = e[iEmt].mother1();
133 
134  // find the top resonance the radiator originates from
135  int iTop = e[iRadBef].mother1();
136  int distance = 1;
137  while (abs(e[iTop].id()) != 6 && iTop > 0) {
138  iTop = e[iTop].mother1();
139  distance ++;
140  }
141  if (iTop == 0) {
142  infoPtr->errorMsg("Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from top quark, not vetoing");
143  return doVetoFSR(false,0,0);
144  //return false;
145  }
146  int iTopCharge = (e[iTop].id()>0)?1:-1;
147 
148  // calculate the scale of the emission
149  double scale;
150  //using pythia pT definition ...
151  if(pTpythiaVeto)
152  scale = pTpythia(e, iRadAft, iEmt, iRecAft);
153  //.. or using POWHEG pT definition
154  else{
155  Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pt(e[iTop].p()), prec(e[iRecAft].p()), psystem;
156  // The computation of the POWHEG pT can be done in the top rest frame or in the diple one.
157  // pdipole = pemt +prec +prad (after the emission)
158  // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop
159  if(vetoDipoleFrame)
160  psystem = pr+pe+prec;
161  else
162  psystem = pt;
163 
164  // gluon splitting into two partons
165  if (e[iRadBef].id() == 21)
166  scale = gSplittingScale(psystem, pr, pe);
167  // quark emitting a gluon (or a photon)
168  else if (abs(e[iRadBef].id()) == 5 && ((e[iEmt].id() == 21) && ! vetoQED) )
169  scale = qSplittingScale(psystem, pr, pe);
170  // other stuff (which we should not veto)
171  else
172  scale = 0;
173  }
174 
175  if (iTopCharge > 0) {
176  if (onlyDistance1) {
177  if ( debug && (distance == 1) && scale > topresscale && ! wouldVetoFsr)
178  cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
179  return doVetoFSR((distance == 1) && scale > topresscale,scale,iTopCharge);
180  }
181  else {
182  if ( debug && scale > topresscale && ! wouldVetoFsr)
183  cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
184  return doVetoFSR(scale > topresscale,scale,iTopCharge);
185  }
186  }
187  else if (iTopCharge < 0){
188  if (onlyDistance1){
189  if ( debug && (distance == 1) && scale > atopresscale && ! wouldVetoFsr)
190  cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
191  return doVetoFSR((distance == 1) && scale > atopresscale,scale,iTopCharge);
192  }
193  else {
194  if ( debug && scale > topresscale && ! wouldVetoFsr)
195  cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id() << "; " << scale << endl;
196  return doVetoFSR(scale > atopresscale,scale,iTopCharge);
197  }
198  }
199  else {
200  cout << "Bug in PohwgeHooksBB4l" << endl;
201  }
202  }
204  // VETO THE PRODUCTION PROCESS //
206  // covered by multiuserhook, i.e. need to turn on EV1
207  // else if(!inResonance && vetoProduction){
208  // return EmissionVetoHook1::doVetoFSREmission(sizeOld, e, iSys, inResonance);
209  // }
210 
211  return false;
212  }
213 
214  inline bool doVetoFSR(bool condition, double scale, int iTopCharge) {
215  if(radtype_.radtype==2)
216  return false;
217  if (condition) {
218  if (!wouldVetoFsr) {
219  wouldVetoFsr = true;
221  vetoTopCharge = iTopCharge;
222  }
223  if (dryRunFSR) return false;
224  else {
225  nFSRvetoBB4l++;
226  return true;
227  }
228  }
229  else return false;
230  }
231 
232  inline void doVetoFSRInit() {
233  wouldVetoFsr = false;
234  vetoDecScale = -1;
235  vetoTopCharge = 0;
236  }
237 
238  //--- SCALE RESONANCE HOOK -------------------------------------------------
239  // called before each resonance decay shower
240  inline bool canSetResonanceScale() override { return scaleResonanceVeto; }
241  // if the resonance is the (anti)top set the scale to:
242  // ---> (anti)top virtuality if radtype=2
243  // ---> (a)topresscale otherwise
244  // if is not the top, set it to a big number
245  inline double scaleResonance(int iRes, const Event &e) override {
246  if (e[iRes].id() == 6){
247  if(radtype_.radtype == 2)
248  return sqrt(e[iRes].m2Calc());
249  else
250  return topresscale;
251  }
252  else if (e[iRes].id() == -6){
253  if(radtype_.radtype == 2)
254  return sqrt(e[iRes].m2Calc());
255  else
256  return atopresscale;
257  }
258  else
259  return pow(10.0,30.);
260  }
261 
262 
263  //--- Internal helper functions --------------------------------------------
264 
265  // Calculates the scale of the hardest emission from within the resonance system
266  // translated by Markus Seidel modified by Tomas Jezo
267  inline double findresscale( const int iRes, const Event& event) {
268  double scale = 0.;
269 
270  int nDau = event[iRes].daughterList().size();
271 
272  if (nDau == 0) {
273  // No resonance found, set scale to high value
274  // Pythia will shower any MC generated resonance unrestricted
275  scale = 1e30;
276  }
277  else if (nDau < 3) {
278  // No radiating resonance found
279  scale = pTmin;
280  }
281  else if (abs(event[iRes].id()) == 6) {
282  // Find top daughters
283  int idw = -1, idb = -1, idg = -1;
284 
285  for (int i = 0; i < nDau; i++) {
286  int iDau = event[iRes].daughterList()[i];
287  if (abs(event[iDau].id()) == 24) idw = iDau;
288  if (abs(event[iDau].id()) == 5) idb = iDau;
289  if (abs(event[iDau].id()) == 21) idg = iDau;
290  }
291 
292  // Get daughter 4-vectors in resonance frame
293  Vec4 pw(event[idw].p());
294  pw.bstback(event[iRes].p());
295 
296  Vec4 pb(event[idb].p());
297  pb.bstback(event[iRes].p());
298 
299  Vec4 pg(event[idg].p());
300  pg.bstback(event[iRes].p());
301 
302  // Calculate scale
303  scale = sqrt(2*pg*pb*pg.e()/pb.e());
304  }
305  else {
306  scale = 1e30;
307  }
308 
309  return scale;
310  }
311 
312  // The following routine will match daughters of particle `e[iparticle]` to an expected pattern specified via the list of expected particle PDG ID's `ids`,
313  // id wildcard is specified as 0 if match is obtained, the positions and the momenta of these particles are returned in vectors `positions` and `momenta`
314  // respectively
315  // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less
316  inline bool match_decay(int iparticle, const Event &e, const vector<int> &ids, vector<int> &positions, vector<Vec4> &momenta, bool exitOnExtraLegs = true){
317  // compare sizes
318  if (e[iparticle].daughterList().size() != ids.size()) {
319  if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) {
320  cout << "extra leg" << endl;
321  exit(-1);
322  }
323  return false;
324  }
325  // compare content
326  for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) {
327  int di = e[iparticle].daughterList()[i];
328  if (ids[i] != 0 && e[di].id() != ids[i])
329  return false;
330  }
331  // reset the positions and momenta vectors (because they may be reused)
332  positions.clear();
333  momenta.clear();
334  // construct the array of momenta
335  for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) {
336  int di = e[iparticle].daughterList()[i];
337  positions.push_back(di);
338  momenta.push_back(e[di].p());
339  }
340  return true;
341  }
342 
343  inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){
344  p1.bstback(pt);
345  p2.bstback(pt);
346  return sqrt( 2*p1*p2*p2.e()/p1.e() );
347  }
348 
349  inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2){
350  p1.bstback(pt);
351  p2.bstback(pt);
352  return sqrt( 2*p1*p2*p1.e()*p2.e()/(pow(p1.e()+p2.e(),2)) );
353  }
354 
355  // Routines to calculate the pT (according to pTdefMode) in a FS splitting:
356  // i (radiator before) -> j (emitted after) k (radiator after)
357  // For the Pythia pT definition, a recoiler (after) must be specified.
358  // (INSPIRED BY pythia8F77_31.cc double pTpythia)
359  inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch,
360  int RecAfterBranch)
361  {
362 
363  // Convenient shorthands for later
364  Vec4 radVec = e[RadAfterBranch].p();
365  Vec4 emtVec = e[EmtAfterBranch].p();
366  Vec4 recVec = e[RecAfterBranch].p();
367  int radID = e[RadAfterBranch].id();
368 
369  // Calculate virtuality of splitting
370  Vec4 Q(radVec + emtVec);
371  double Qsq = Q.m2Calc();
372 
373 
374  // Mass term of radiator
375  double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
376  pow2(particleDataPtr->m0(radID)) : 0.;
377 
378  // z values for FSR
379  double z, pTnow;
380  // Construct 2 -> 3 variables
381  Vec4 sum = radVec + recVec + emtVec;
382  double m2Dip = sum.m2Calc();
383 
384  double x1 = 2. * (sum * radVec) / m2Dip;
385  double x3 = 2. * (sum * emtVec) / m2Dip;
386  z = x1 / (x1 + x3);
387  pTnow = z * (1. - z);
388 
389 
390  // Virtuality
391  pTnow *= (Qsq - m2Rad);
392 
393  if (pTnow < 0.) {
394  cout << "Warning: pTpythia was negative" << endl;
395  return -1.;
396  }
397  else
398  return(sqrt(pTnow));
399  }
400 
401  inline double getdechardness(int topcharge, const Event &e){
402  int tid = 6*topcharge, wid = 24*topcharge, bid = 5*topcharge, gid = 21, wildcard = 0;
403  // find last top in the record
404  int i_top = -1;
405  Vec4 p_top, p_b, p_g, p_g1, p_g2;
406  for (int i = 0; i < e.size(); i++)
407  if (e[i].id() == tid) {
408  i_top = i;
409  p_top = e[i].p();
410  }
411  if (i_top == -1) return -1.0;
412 
413  // summary of cases
414  // 1.) t > W b
415  // a.) b > 3 ... error
416  // b.) b > b g ... h = sqrt(2*p_g*p_b*p_g.e()/p_b.e())
417  // c.) b > other ... h = -1
418  // return h
419  // 2.) t > W b g
420  // a.) b > 3 ... error
421  // b.) b > b g ... h1 = sqrt(2*p_g*p_b*p_g.e()/p_b.e())
422  // c.) b > other ... h1 = -1
423  // i.) g > 3 ... error
424  // ii.) g > 2 ... h2 = sqrt(2*p_g1*p_g2*p_g1.e()*p_g2.e()/(pow(p_g1.e(),2)+pow(p_g2.e(),2))) );
425  // iii.) g > other ... h2 = -1
426  // return max(h1,h2)
427  // 3.) else ... error
428 
429  vector<Vec4> momenta;
430  vector<int> positions;
431 
432  // 1.) t > b W
433  if ( match_decay(i_top, e, vector<int> {wid, bid}, positions, momenta, false) ) {
434  double h;
435  int i_b = positions[1];
436  // a.+b.) b > 3 or b > b g
437  if ( match_decay(i_b, e, vector<int> {bid, gid}, positions, momenta) )
438  h = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
439  // c.) b > other
440  else
441  h = -1;
442  return h;
443  }
444  // 2.) t > b W g
445  else if ( match_decay(i_top, e, vector<int> {wid, bid, gid}, positions, momenta, false) ) {
446  double h1, h2;
447  int i_b = positions[1], i_g = positions[2];
448  // a.+b.) b > 3 or b > b g
449  if ( match_decay(i_b, e, vector<int> {bid, gid}, positions, momenta) )
450  h1 = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
451  // c.) b > other
452  else
453  h1 = -1;
454  // i.+ii.) g > 3 or g > 2
455  if ( match_decay(i_g, e, vector<int> {wildcard, wildcard}, positions, momenta) )
456  h2 = gSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
457  // c.) b > other
458  else
459  h2 = -1;
460  return max(h1, h2);
461  }
462  // 3.) else
463  else {
464  cout << "getdechardness" << endl;
465  cout << "top at position " << i_top << endl;
466  cout << "with " << e[i_top].daughterList().size() << " daughters " << endl;
467  for (unsigned i = 0; i < e[i_top].daughterList().size(); i++) {
468  int di = e[i_top].daughterList()[i];
469  cout << "with daughter " << di << ": " << e[di].id() << endl;
470  }
471  exit(-1);
472  }
473  }
474 
475 
476 
477  //--------------------------------------------------------------------------
478 
479  // Functions to return information
480 
481  // inline int getNFSRveto() { return nFSRveto; }
482 
483  //--------------------------------------------------------------------------
484 
485  private:
486  // FSR emission veto flags
488  // Parton Level veto flags
490  // Scale Resonance veto flags
492  // other flags
493  bool debug;
494  // internal: resonance scales
496  // internal: inter veto communication
497  double vetoDecScale;
501  //bool vetoProduction;
502  double pTmin;
503  // Statistics on vetos
504  unsigned long int nFSRvetoBB4l;
505 
506  };
507 
508  //==========================================================================
509 
510 } // end namespace Pythia8
511 
512 #endif // end Pythia8_PowhegHooksBB4L_H
size
Write out results.
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
struct @706 radtype_
ExtVec< T, 4 > Vec4
Definition: ExtVec.h:62
double scaleResonance(int iRes, const Event &e) override
double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch)
double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2)
bool initAfterBeams() override
double getdechardness(int topcharge, const Event &e)
bool canVetoProcessLevel() override
bool doVetoFSR(bool condition, double scale, int iTopCharge)
bool canVetoFSREmission() override
T sqrt(T t)
Definition: SSEVec.h:18
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double p2[4]
Definition: TauolaWrapper.h:90
double findresscale(const int iRes, const Event &event)
double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2)
bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) override
bool canVetoPartonLevel() override
int radtype
bool doVetoProcessLevel(Event &e) override
bool retryPartonLevel() override
bool canSetResonanceScale() override
bool doVetoPartonLevel(const Event &e) override
double p1[4]
Definition: TauolaWrapper.h:89
bool match_decay(int iparticle, const Event &e, const vector< int > &ids, vector< int > &positions, vector< Vec4 > &momenta, bool exitOnExtraLegs=true)
unsigned long int nFSRvetoBB4l
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: event.py:1