Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:27

0001 // PowhegHooksBB4L.h
0002 // Copyright (C) 2017 Silvia Ferrario Ravasio, Tomas Jezo, Paolo Nason, Markus Seidel
0003 // inspired by PowhegHooks.h by Richard Corke
0004 // adjusted to work with EmissionVetoHook1 in CMSSW by Alexander Grohsjean
0005 
0006 #ifndef Pythia8_PowhegHooksBB4L_H
0007 #define Pythia8_PowhegHooksBB4L_H
0008 
0009 // Includes
0010 #include "Pythia8/Pythia.h"
0011 #include <cassert>
0012 struct {
0013   int radtype;
0014 } radtype_;
0015 
0016 namespace Pythia8 {
0017 
0018   class PowhegHooksBB4L : public UserHooks {
0019   public:
0020     //--- Constructor and destructor -------------------------------------------
0021     PowhegHooksBB4L() : nFSRvetoBB4l(0) {}
0022     ~PowhegHooksBB4L() override { std::cout << "Number of FSR vetoed in BB4l = " << nFSRvetoBB4l << std::endl; }
0023 
0024     //--- Initialization -----------------------------------------------------------------------
0025     bool initAfterBeams() override {
0026       // settings of this class
0027       vetoFSREmission = settingsPtr->flag("POWHEG:bb4l:FSREmission:veto");
0028       onlyDistance1 = settingsPtr->flag("POWHEG:bb4l:FSREmission:onlyDistance1");
0029       dryRunFSR = settingsPtr->flag("POWHEG:bb4l:FSREmission:dryRun");
0030       vetoAtPL = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoAtPL");
0031       vetoQED = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoQED");
0032       vetoPartonLevel = settingsPtr->flag("POWHEG:bb4l:PartonLevel:veto");
0033       excludeFSRConflicting = settingsPtr->flag("POWHEG:bb4l:PartonLevel:excludeFSRConflicting");
0034       debug = settingsPtr->flag("POWHEG:bb4l:DEBUG");
0035       scaleResonanceVeto = settingsPtr->flag("POWHEG:bb4l:ScaleResonance:veto");
0036       vetoDipoleFrame = settingsPtr->flag("POWHEG:bb4l:FSREmission:vetoDipoleFrame");
0037       pTpythiaVeto = settingsPtr->flag("POWHEG:bb4l:FSREmission:pTpythiaVeto");
0038       //vetoProduction = (settingsPtr->mode("POWHEG:veto")==1);
0039       pTmin = settingsPtr->parm("POWHEG:bb4l:pTminVeto");
0040       return true;
0041     }
0042 
0043     //--- PROCESS LEVEL HOOK ---------------------------------------------------
0044 
0045     // called at the LHE level
0046     inline bool canVetoProcessLevel() override { return true; }
0047     inline bool doVetoProcessLevel(Event &e) override {
0048       // extract the radtype from the event comment
0049       stringstream ss;
0050       // use eventattribute as comments not filled when using edm input
0051       //ss << infoPtr->getEventComments();
0052       ss << infoPtr->getEventAttribute("#rwgt");
0053       string temp;
0054       ss >> temp >> radtype_.radtype;
0055       assert(temp == "#rwgt");
0056 
0057       // find last top and the last anti-top in the record
0058       int i_top = -1, i_atop = -1;
0059       for (int i = 0; i < e.size(); i++) {
0060         if (e[i].id() == 6)
0061           i_top = i;
0062         if (e[i].id() == -6)
0063           i_atop = i;
0064       }
0065       if (i_top != -1)
0066         topresscale = findresscale(i_top, e);
0067       else
0068         topresscale = 1e30;
0069       if (i_top != -1)
0070         atopresscale = findresscale(i_atop, e);
0071       else
0072         atopresscale = 1e30;
0073       // initialize stuff
0074       doVetoFSRInit();
0075       // do not veto, ever
0076       return false;
0077     }
0078 
0079     //--- PARTON LEVEL HOOK ----------------------------------------------------
0080 
0081     // called after shower
0082     bool retryPartonLevel() override { return vetoPartonLevel || vetoAtPL; }
0083     inline bool canVetoPartonLevel() override { return vetoPartonLevel || vetoAtPL; }
0084     inline bool doVetoPartonLevel(const Event &e) override {
0085       if (radtype_.radtype == 2)
0086         return false;
0087       if (debug) {
0088         if (dryRunFSR && wouldVetoFsr) {
0089           double scale = getdechardness(vetoTopCharge, e);
0090           cout << "FSRdecScale = " << vetoDecScale << ", PLdecScale = " << scale << ", ratio = " << vetoDecScale / scale
0091                << endl;
0092         }
0093       }
0094       if (vetoPartonLevel) {
0095         double topdecscale = getdechardness(1, e);
0096         double atopdecscale = getdechardness(-1, e);
0097         if ((topdecscale > topresscale) || (atopdecscale > atopresscale)) {
0098           //if(dryRunFSR && ! wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0;
0099           return true;
0100         } else
0101           //if(dryRunFSR && wouldVetoFsr) mydatacontainer_.excludeEvent = excludeFSRConflicting?1:0;
0102           return false;
0103       }
0104       if (vetoAtPL) {
0105         if (dryRunFSR && wouldVetoFsr)
0106           return true;
0107         else
0108           return false;
0109       }
0110       return false;
0111     }
0112 
0113     //--- FSR EMISSION LEVEL HOOK ----------------------------------------------
0114 
0115     // FSR veto: this should be true if we want PowhegHooksBB4l veto in decay
0116     //           OR PowhegHooks veto in production. (The virtual method
0117     //           PowhegHooks::canVetoFSREmission has been replaced by
0118     //           PowhegHooksBB4L::canVetoFSREmission, so FSR veto in production
0119     //           must be handled here. ISR and MPI veto are instead still
0120     //           handled by PowhegHooks.)
0121     inline bool canVetoFSREmission() override { return vetoFSREmission; }  // || vetoProduction; }
0122     inline bool doVetoFSREmission(int sizeOld, const Event &e, int iSys, bool inResonance) override {
0123       //////////////////////////////
0124       //VETO INSIDE THE RESONANCE //
0125       //////////////////////////////
0126       if (inResonance && vetoFSREmission) {
0127         int iRecAft = e.size() - 1;
0128         int iEmt = e.size() - 2;
0129         int iRadAft = e.size() - 3;
0130         int iRadBef = e[iEmt].mother1();
0131 
0132         // find the top resonance the radiator originates from
0133         int iTop = e[iRadBef].mother1();
0134         int distance = 1;
0135         while (abs(e[iTop].id()) != 6 && iTop > 0) {
0136           iTop = e[iTop].mother1();
0137           distance++;
0138         }
0139         if (iTop == 0) {
0140           infoPtr->errorMsg(
0141               "Warning in PowhegHooksBB4L::doVetoFSREmission: emission in resonance not from top quark, not vetoing");
0142           return doVetoFSR(false, 0, 0);
0143           //return false;
0144         }
0145         int iTopCharge = (e[iTop].id() > 0) ? 1 : -1;
0146 
0147         // calculate the scale of the emission
0148         double scale;
0149         //using pythia pT definition ...
0150         if (pTpythiaVeto)
0151           scale = pTpythia(e, iRadAft, iEmt, iRecAft);
0152         //.. or using POWHEG pT definition
0153         else {
0154           Vec4 pr(e[iRadAft].p()), pe(e[iEmt].p()), pt(e[iTop].p()), prec(e[iRecAft].p()), psystem;
0155           // The computation of the POWHEG pT can be done in the top rest frame or in the diple one.
0156           // pdipole = pemt +prec +prad (after the emission)
0157           // For the first emission off the top resonance pdipole = pw +pb (before the emission) = ptop
0158           if (vetoDipoleFrame)
0159             psystem = pr + pe + prec;
0160           else
0161             psystem = pt;
0162 
0163           // gluon splitting into two partons
0164           if (e[iRadBef].id() == 21)
0165             scale = gSplittingScale(psystem, pr, pe);
0166           // quark emitting a gluon (or a photon)
0167           else if (abs(e[iRadBef].id()) == 5 && ((e[iEmt].id() == 21) && !vetoQED))
0168             scale = qSplittingScale(psystem, pr, pe);
0169           // other stuff (which we should not veto)
0170           else
0171             scale = 0;
0172         }
0173 
0174         if (iTopCharge > 0) {
0175           if (onlyDistance1) {
0176             if (debug && (distance == 1) && scale > topresscale && !wouldVetoFsr)
0177               cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id()
0178                    << "; " << scale << endl;
0179             return doVetoFSR((distance == 1) && scale > topresscale, scale, iTopCharge);
0180           } else {
0181             if (debug && scale > topresscale && !wouldVetoFsr)
0182               cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id()
0183                    << "; " << scale << endl;
0184             return doVetoFSR(scale > topresscale, scale, iTopCharge);
0185           }
0186         } else if (iTopCharge < 0) {
0187           if (onlyDistance1) {
0188             if (debug && (distance == 1) && scale > atopresscale && !wouldVetoFsr)
0189               cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id()
0190                    << "; " << scale << endl;
0191             return doVetoFSR((distance == 1) && scale > atopresscale, scale, iTopCharge);
0192           } else {
0193             if (debug && scale > topresscale && !wouldVetoFsr)
0194               cout << e[iTop].id() << ": " << e[iRadBef].id() << " > " << e[iRadAft].id() << " + " << e[iEmt].id()
0195                    << "; " << scale << endl;
0196             return doVetoFSR(scale > atopresscale, scale, iTopCharge);
0197           }
0198         } else {
0199           cout << "Bug in PohwgeHooksBB4l" << endl;
0200         }
0201       }
0202       /////////////////////////////////
0203       // VETO THE PRODUCTION PROCESS //
0204       /////////////////////////////////
0205       // covered by multiuserhook, i.e. need to turn on EV1
0206       // else if(!inResonance && vetoProduction){
0207       // return EmissionVetoHook1::doVetoFSREmission(sizeOld, e, iSys, inResonance);
0208       // }
0209 
0210       return false;
0211     }
0212 
0213     inline bool doVetoFSR(bool condition, double scale, int iTopCharge) {
0214       if (radtype_.radtype == 2)
0215         return false;
0216       if (condition) {
0217         if (!wouldVetoFsr) {
0218           wouldVetoFsr = true;
0219           vetoDecScale = scale;
0220           vetoTopCharge = iTopCharge;
0221         }
0222         if (dryRunFSR)
0223           return false;
0224         else {
0225           nFSRvetoBB4l++;
0226           return true;
0227         }
0228       } else
0229         return false;
0230     }
0231 
0232     inline void doVetoFSRInit() {
0233       wouldVetoFsr = false;
0234       vetoDecScale = -1;
0235       vetoTopCharge = 0;
0236     }
0237 
0238     //--- SCALE RESONANCE HOOK -------------------------------------------------
0239     // called before each resonance decay shower
0240     inline bool canSetResonanceScale() override { return scaleResonanceVeto; }
0241     // if the resonance is the (anti)top set the scale to:
0242     //  ---> (anti)top virtuality if radtype=2
0243     //  ---> (a)topresscale otherwise
0244     // if is not the top, set it to a big number
0245     inline double scaleResonance(int iRes, const Event &e) override {
0246       if (e[iRes].id() == 6) {
0247         if (radtype_.radtype == 2)
0248           return sqrt(e[iRes].m2Calc());
0249         else
0250           return topresscale;
0251       } else if (e[iRes].id() == -6) {
0252         if (radtype_.radtype == 2)
0253           return sqrt(e[iRes].m2Calc());
0254         else
0255           return atopresscale;
0256       } else
0257         return pow(10.0, 30.);
0258     }
0259 
0260     //--- Internal helper functions --------------------------------------------
0261 
0262     // Calculates the scale of the hardest emission from within the resonance system
0263     // translated by Markus Seidel modified by Tomas Jezo
0264     inline double findresscale(const int iRes, const Event &event) {
0265       double scale = 0.;
0266 
0267       int nDau = event[iRes].daughterList().size();
0268 
0269       if (nDau == 0) {
0270         // No resonance found, set scale to high value
0271         // Pythia will shower any MC generated resonance unrestricted
0272         scale = 1e30;
0273       } else if (nDau < 3) {
0274         // No radiating resonance found
0275         scale = pTmin;
0276       } else if (abs(event[iRes].id()) == 6) {
0277         // Find top daughters
0278         int idw = -1, idb = -1, idg = -1;
0279 
0280         for (int i = 0; i < nDau; i++) {
0281           int iDau = event[iRes].daughterList()[i];
0282           if (abs(event[iDau].id()) == 24)
0283             idw = iDau;
0284           if (abs(event[iDau].id()) == 5)
0285             idb = iDau;
0286           if (abs(event[iDau].id()) == 21)
0287             idg = iDau;
0288         }
0289 
0290         // Get daughter 4-vectors in resonance frame
0291         Vec4 pw(event[idw].p());
0292         pw.bstback(event[iRes].p());
0293 
0294         Vec4 pb(event[idb].p());
0295         pb.bstback(event[iRes].p());
0296 
0297         Vec4 pg(event[idg].p());
0298         pg.bstback(event[iRes].p());
0299 
0300         // Calculate scale
0301         scale = sqrt(2 * pg * pb * pg.e() / pb.e());
0302       } else {
0303         scale = 1e30;
0304       }
0305 
0306       return scale;
0307     }
0308 
0309     // 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`,
0310     // 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`
0311     // respectively
0312     // if exitOnExtraLegs==true, it will exit if the decay has more particles than expected, but not less
0313     inline bool match_decay(int iparticle,
0314                             const Event &e,
0315                             const vector<int> &ids,
0316                             vector<int> &positions,
0317                             vector<Vec4> &momenta,
0318                             bool exitOnExtraLegs = true) {
0319       // compare sizes
0320       if (e[iparticle].daughterList().size() != ids.size()) {
0321         if (exitOnExtraLegs && e[iparticle].daughterList().size() > ids.size()) {
0322           cout << "extra leg" << endl;
0323           exit(-1);
0324         }
0325         return false;
0326       }
0327       // compare content
0328       for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) {
0329         int di = e[iparticle].daughterList()[i];
0330         if (ids[i] != 0 && e[di].id() != ids[i])
0331           return false;
0332       }
0333       // reset the positions and momenta vectors (because they may be reused)
0334       positions.clear();
0335       momenta.clear();
0336       // construct the array of momenta
0337       for (unsigned i = 0; i < e[iparticle].daughterList().size(); i++) {
0338         int di = e[iparticle].daughterList()[i];
0339         positions.push_back(di);
0340         momenta.push_back(e[di].p());
0341       }
0342       return true;
0343     }
0344 
0345     inline double qSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2) {
0346       p1.bstback(pt);
0347       p2.bstback(pt);
0348       return sqrt(2 * p1 * p2 * p2.e() / p1.e());
0349     }
0350 
0351     inline double gSplittingScale(Vec4 pt, Vec4 p1, Vec4 p2) {
0352       p1.bstback(pt);
0353       p2.bstback(pt);
0354       return sqrt(2 * p1 * p2 * p1.e() * p2.e() / (pow(p1.e() + p2.e(), 2)));
0355     }
0356 
0357     // Routines to calculate the pT (according to pTdefMode) in a FS splitting:
0358     // i (radiator before) -> j (emitted after) k (radiator after)
0359     // For the Pythia pT definition, a recoiler (after) must be specified.
0360     // (INSPIRED BY pythia8F77_31.cc double pTpythia)
0361     inline double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch) {
0362       // Convenient shorthands for later
0363       Vec4 radVec = e[RadAfterBranch].p();
0364       Vec4 emtVec = e[EmtAfterBranch].p();
0365       Vec4 recVec = e[RecAfterBranch].p();
0366       int radID = e[RadAfterBranch].id();
0367 
0368       // Calculate virtuality of splitting
0369       Vec4 Q(radVec + emtVec);
0370       double Qsq = Q.m2Calc();
0371 
0372       // Mass term of radiator
0373       double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? pow2(particleDataPtr->m0(radID)) : 0.;
0374 
0375       // z values for FSR
0376       double z, pTnow;
0377       // Construct 2 -> 3 variables
0378       Vec4 sum = radVec + recVec + emtVec;
0379       double m2Dip = sum.m2Calc();
0380 
0381       double x1 = 2. * (sum * radVec) / m2Dip;
0382       double x3 = 2. * (sum * emtVec) / m2Dip;
0383       z = x1 / (x1 + x3);
0384       pTnow = z * (1. - z);
0385 
0386       // Virtuality
0387       pTnow *= (Qsq - m2Rad);
0388 
0389       if (pTnow < 0.) {
0390         cout << "Warning: pTpythia was negative" << endl;
0391         return -1.;
0392       } else
0393         return (sqrt(pTnow));
0394     }
0395 
0396     inline double getdechardness(int topcharge, const Event &e) {
0397       int tid = 6 * topcharge, wid = 24 * topcharge, bid = 5 * topcharge, gid = 21, wildcard = 0;
0398       // find last top in the record
0399       int i_top = -1;
0400       Vec4 p_top, p_b, p_g, p_g1, p_g2;
0401       for (int i = 0; i < e.size(); i++)
0402         if (e[i].id() == tid) {
0403           i_top = i;
0404           p_top = e[i].p();
0405         }
0406       if (i_top == -1)
0407         return -1.0;
0408 
0409       // summary of cases
0410       // 1.) t > W b
0411       //   a.) b > 3     ... error
0412       //   b.) b > b g   ... h = sqrt(2*p_g*p_b*p_g.e()/p_b.e())
0413       //   c.) b > other ... h = -1
0414       //   return h
0415       // 2.) t > W b g
0416       //   a.)   b > 3     ... error
0417       //   b.)   b > b g   ... h1 = sqrt(2*p_g*p_b*p_g.e()/p_b.e())
0418       //   c.)   b > other ... h1 = -1
0419       //   i.)   g > 3     ... error
0420       //   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))) );
0421       //   iii.) g > other ... h2 = -1
0422       //   return max(h1,h2)
0423       // 3.) else ... error
0424 
0425       vector<Vec4> momenta;
0426       vector<int> positions;
0427 
0428       // 1.) t > b W
0429       if (match_decay(i_top, e, vector<int>{wid, bid}, positions, momenta, false)) {
0430         double h;
0431         int i_b = positions[1];
0432         // a.+b.) b > 3 or b > b g
0433         if (match_decay(i_b, e, vector<int>{bid, gid}, positions, momenta))
0434           h = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
0435         // c.) b > other
0436         else
0437           h = -1;
0438         return h;
0439       }
0440       // 2.) t > b W g
0441       else if (match_decay(i_top, e, vector<int>{wid, bid, gid}, positions, momenta, false)) {
0442         double h1, h2;
0443         int i_b = positions[1], i_g = positions[2];
0444         // a.+b.) b > 3 or b > b g
0445         if (match_decay(i_b, e, vector<int>{bid, gid}, positions, momenta))
0446           h1 = qSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
0447         // c.) b > other
0448         else
0449           h1 = -1;
0450         // i.+ii.) g > 3 or g > 2
0451         if (match_decay(i_g, e, vector<int>{wildcard, wildcard}, positions, momenta))
0452           h2 = gSplittingScale(e[i_top].p(), momenta[0], momenta[1]);
0453         // c.) b > other
0454         else
0455           h2 = -1;
0456         return max(h1, h2);
0457       }
0458       // 3.) else
0459       else {
0460         cout << "getdechardness" << endl;
0461         cout << "top at position " << i_top << endl;
0462         cout << "with " << e[i_top].daughterList().size() << " daughters " << endl;
0463         for (unsigned i = 0; i < e[i_top].daughterList().size(); i++) {
0464           int di = e[i_top].daughterList()[i];
0465           cout << "with daughter " << di << ": " << e[di].id() << endl;
0466         }
0467         exit(-1);
0468       }
0469     }
0470 
0471     //--------------------------------------------------------------------------
0472 
0473     // Functions to return information
0474 
0475     //  inline int    getNFSRveto() { return nFSRveto; }
0476 
0477     //--------------------------------------------------------------------------
0478 
0479   private:
0480     // FSR emission veto flags
0481     bool vetoFSREmission, dryRunFSR, wouldVetoFsr, onlyDistance1, vetoAtPL, vetoQED;
0482     // Parton Level veto flags
0483     bool vetoPartonLevel, excludeFSRConflicting;
0484     // Scale Resonance veto flags
0485     double scaleResonanceVeto;
0486     // other flags
0487     bool debug;
0488     // internal: resonance scales
0489     double topresscale, atopresscale;
0490     // internal: inter veto communication
0491     double vetoDecScale;
0492     int vetoTopCharge;
0493     bool vetoDipoleFrame;
0494     bool pTpythiaVeto;
0495     //bool vetoProduction;
0496     double pTmin;
0497     // Statistics on vetos
0498     unsigned long int nFSRvetoBB4l;
0499   };
0500 
0501   //==========================================================================
0502 
0503 }  // end namespace Pythia8
0504 
0505 #endif  // end Pythia8_PowhegHooksBB4L_H