GenFilters/plugins/ComphepSingletopFilterPy8.cc

0001 // Package:    GenFilters
0002 // Class:      ComphepSingletopFilterPy8
0003 //
0004 /*
0005  class ComphepSingletopFilterPy8 ComphepSingletopFilterPy8.cc GeneratorInterface/GenFilters/src/ComphepSingletopFilterPy8.cc
0006
0007  Description: a filter to match LO/NLO in comphep-generated singletop
0008
0009  Implementation:
0010      <Notes on implementation>
0011 */
0012 //
0013 // Original Author:  Vladimir Molchanov
0014 //         Created:  Wed Mar 25 19:43:12 CET 2009
0015 // $Id: ComphepSingletopFilterPy8.cc,v 1.3 2009/12/15 10:29:32 fabiocos Exp $
0016 // Author:  Alexey Baskakov
0017 //         Created:  Oct 2014
0018 //  ComphepSingletopFilterPy8.cc,v 2.1
0019 //
0020
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/global/EDFilter.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/Utilities/interface/EDGetToken.h"
0029 #include "FWCore/Utilities/interface/InputTag.h"
0030 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0031
0032 #include <atomic>
0033 #include <cmath>
0034 #include <cstdlib>
0035 #include <sstream>
0036 #include <string>
0037 #include <vector>
0038
0039 class ComphepSingletopFilterPy8 : public edm::global::EDFilter<> {
0040 public:
0041   explicit ComphepSingletopFilterPy8(const edm::ParameterSet &);
0042
0043 private:
0044   bool filter(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0045   void endJob() override;
0046
0047   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0048   const double ptsep;
0049   mutable std::atomic<int> read22;
0050   mutable std::atomic<int> read23;
0051   mutable std::atomic<int> pass22;
0052   mutable std::atomic<int> pass23;
0053   const int hardLep;
0054 };
0055
0056 using namespace std;
0057 using namespace HepMC;
0058
0059 ComphepSingletopFilterPy8::ComphepSingletopFilterPy8(const edm::ParameterSet &iConfig)
0060     : token_(consumes<edm::HepMCProduct>(
0061           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0062       ptsep(iConfig.getParameter<double>("pTSep")),
0063       read22(0),
0064       read23(0),
0065       pass22(0),
0066       pass23(0),
0067       hardLep(23) {}  //identifies the "hard part" in Pythia8
0068
0069 void ComphepSingletopFilterPy8::endJob() {
0070   std::ostringstream ss;
0071   ss << "Proc: " << std::setw(9) << "2-->2" << std::setw(9) << "2-->3" << std::setw(9) << "Total"
0072      << "\n"
0073      << "Read: " << std::setw(9) << read22 << std::setw(9) << read23 << std::setw(9) << read22 + read23 << "\n"
0074      << "Pass: " << std::setw(9) << pass22 << std::setw(9) << pass23 << std::setw(9) << pass22 + pass23 << "\n";
0075   edm::LogAbsolute("ComphepSingletopFilterPy8") << ss.str();
0076 }
0077
0078 bool ComphepSingletopFilterPy8::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &) const {
0079   edm::Handle<edm::HepMCProduct> evt;
0080   iEvent.getByToken(token_, evt);
0081   const HepMC::GenEvent *myEvt = evt->GetEvent();
0082
0083   int id_bdec = 0, id_lJet = 0, id_b_from_top = 0, id_lep = 0;
0084   //vars for lepton top
0085   const GenParticle *gp_clep = nullptr;
0086   GenVertex *gv_hard = nullptr;
0087
0088   const GenParticle *gp_lep_FSR = nullptr;
0089   GenVertex *gv_lep = nullptr;
0090   vector<const GenParticle *> vgp_lep;
0091
0092   //vars for add b from top
0093   GenVertex *gv = nullptr;
0094   const GenParticle *gp = nullptr;
0095   vector<const GenParticle *> vgp_bsec;
0096
0097   //vars for light Jet (light q)
0098   GenVertex *gv_lJet = nullptr;
0099   const GenParticle *gplJet = nullptr;
0100   vector<const GenParticle *> vgp_lJet;
0101
0102   //vars for b from top
0103   GenVertex *gv_b_t = nullptr;
0104   const GenParticle *gp_b_t = nullptr;
0105   vector<const GenParticle *> vgp_b_t;
0106
0107   //Run through all particles in myEvt, if lepton found saves it gp_clep
0108   for (GenEvent::particle_const_iterator it = myEvt->particles_begin(); it != myEvt->particles_end(); ++it) {
0109     if (abs((*it)->status()) == hardLep) {
0110       int abs_id = abs((*it)->pdg_id());  //11=e, -11=E, 13=mu, -13=Mu, 15=l(tau), -15=L
0111       if (abs_id == 11 || abs_id == 13 || abs_id == 15) {
0112         gp_clep = *it;
0113         id_lep = (*it)->pdg_id();
0114         gv_lep = (*it)->production_vertex();
0115         // Lepton FSR
0116         while (gv_lep) {
0117           gp_lep_FSR = nullptr;
0118           for (GenVertex::particle_iterator it = gv_lep->particles_begin(children);
0119                it != gv_lep->particles_end(children);
0120                ++it) {
0121             if ((*it)->pdg_id() == id_lep) {
0122               if (!gp_lep_FSR || (*it)->momentum().perp2() > gp_lep_FSR->momentum().perp2()) {
0123                 gp_lep_FSR = *it;
0124               }
0125             }
0126           }
0127           if (gp_lep_FSR) {
0128             gv_lep = gp_lep_FSR->end_vertex();
0129             vgp_lep.push_back(gp_lep_FSR);  //vertex of
0130           } else {
0131             gv_lep = nullptr;  //exits the "while" loop
0132           }
0133         }
0134         break;
0135       }
0136     }
0137   }
0138   // Goes to lepton production vertex
0139   gv_hard = gp_clep->production_vertex();
0140
0141   if (!gp_clep) {
0142     edm::LogError("ComphepSingletopFilterPy8") << "ERROR: ComphepSingletopFilterPy8: no charged lepton" << endl;
0143     return false;
0144   }
0145
0146   //Run through lepton production_vertex
0147   for (GenVertex::particle_iterator it = gv_hard->particles_begin(children); it != gv_hard->particles_end(children);
0148        ++it) {
0149     int pdg_id = (*it)->pdg_id();  // KF in Pythia
0150     int abs_id = abs(pdg_id);
0151
0152     //selection of light quark among particles in primary vertex
0153     if (abs_id < 5) {
0154       id_lJet = (*it)->pdg_id();
0155       gv_lJet = (*it)->production_vertex();
0156
0157       while (gv_lJet) {
0158         gplJet = nullptr;
0159         for (GenVertex::particle_iterator it = gv_lJet->particles_begin(children);
0160              it != gv_lJet->particles_end(children);
0161              ++it) {
0162           if ((*it)->pdg_id() == id_lJet) {
0163             if (!gplJet || (*it)->momentum().perp2() > gplJet->momentum().perp2()) {
0164               gplJet = *it;
0165             }
0166           }
0167         }
0168
0169         if (gplJet) {
0170           gv_lJet = gplJet->end_vertex();
0171           vgp_lJet.push_back(gplJet);  //vertex of
0172         } else {
0173           gv_lJet = nullptr;  //exits the "while" loop
0174         }
0175       }
0176     }
0177
0178     if (abs(pdg_id) == 5)  // 5 = b
0179     {
0180       if (pdg_id * (gp_clep->pdg_id()) < 0) {  //b is from top
0181         id_bdec = pdg_id;
0182
0183         id_b_from_top = (*it)->pdg_id();
0184         gv_b_t = (*it)->production_vertex();
0185
0186         while (gv_b_t) {
0187           gp_b_t = nullptr;
0188           for (GenVertex::particle_iterator it = gv_b_t->particles_begin(children);
0189                it != gv_b_t->particles_end(children);
0190                ++it) {
0191             if ((*it)->pdg_id() == id_b_from_top) {
0192               if (!gp_b_t || (*it)->momentum().perp2() > gp_b_t->momentum().perp2()) {
0193                 gp_b_t = *it;
0194               }
0195             }
0196           }
0197
0198           if (gp_b_t) {
0199             gv_b_t = gp_b_t->end_vertex();
0200             vgp_b_t.push_back(gp_b_t);  //vertex of
0201           } else {
0202             gv_b_t = nullptr;  //exits the "while" loop
0203           }
0204         }
0205       } else {  //If process 2-3, then aditional b in the initial state fills
0206         vgp_bsec.push_back(*it);
0207       }
0208     }
0209   }
0210
0211   bool process22 =
0212       (vgp_bsec.empty());  //if there is no aditional b-quark in primary vexrtes, then it is tq-process (2->2)
0213
0214   if (process22) {
0215     for (GenVertex::particle_iterator it = gv_hard->particles_begin(parents); it != gv_hard->particles_end(parents);
0216          ++it)  //Among parents of lepton production vertex(primary vertex) we find b-quark.
0217     {
0218       if ((*it)->pdg_id() == id_bdec) {
0219         gv = (*it)->production_vertex();
0220         break;
0221       }
0222     }
0223     if (!gv) {
0224       edm::LogError("ComphepSingletopFilterPy8")
0225           << "ERROR: ComphepSingletopFilterPy8: HepMC inconsistency (! gv)" << endl;
0226       myEvt->print();
0227       return false;
0228     }
0229   } else {
0230     gv = vgp_bsec.back()->end_vertex();
0231   }
0232
0233   //##Correction for GV vertex, because of additional gluons in ISR
0234   bool WeFoundAdditional_b_quark = false;
0235   int loopCount = 0;
0236   while (WeFoundAdditional_b_quark != true) {
0237     ////we go through b or B quark (not from top, top parent) production vertex
0238     for (GenVertex::particle_iterator it = gv->particles_begin(children); it != gv->particles_end(children); ++it) {
0239       ////if we found b, but anti to ours, than it is additional b quark from ISR
0240       if ((*it)->pdg_id() == -id_bdec) {   //we found right vertex of ISR gluon splitting!
0241         WeFoundAdditional_b_quark = true;  //gv = (*it)->production_vertex();
0242       }
0243     }
0244     if (WeFoundAdditional_b_quark == false) {
0245       ////If we don't find add b quark, we need to go to parents vertex, to find it there
0246       for (GenVertex::particle_iterator it = gv->particles_begin(parents); it != gv->particles_end(parents); ++it) {
0247         if ((*it)->pdg_id() == id_bdec) {
0248           gv = (*it)->production_vertex();
0249         }
0250       }
0251     }
0252     loopCount++;
0253
0254     if (loopCount > 100)  //loop protection, nothing more
0255     {
0256       edm::LogError("ComphepSingletopFilterPy8")
0257           << "ERROR: ComphepSingletopFilterPy8: HepMC inconsistency (No add b vertex found)" << endl;
0258       break;
0259     }
0260   }
0261
0262   while (gv) {
0263     gp = nullptr;
0264
0265     for (GenVertex::particle_iterator it = gv->particles_begin(children); it != gv->particles_end(children); ++it) {
0266       if ((*it)->pdg_id() == -id_bdec) {
0267         if (!gp || (*it)->momentum().perp2() > gp->momentum().perp2()) {
0268           gp = *it;
0269         }
0270       }
0271     }
0272     if (gp) {
0273       gv = gp->end_vertex();
0274       vgp_bsec.push_back(gp);  //vertex of
0275     } else {
0276       gv = nullptr;  //exits the "while" loop
0277     }
0278   }
0279
0280   if (vgp_bsec.empty()) {
0281     edm::LogError("ComphepSingletopFilterPy8")
0282         << "ERROR: ComphepSingletopFilterPy8: HepMC inconsistency (vgp_bsec.size() == 0)" << endl;
0283     return false;
0284   }
0285
0286   double pt = vgp_bsec.back()->momentum().perp();
0287   // double eta = vgp_bsec.back()->momentum().eta();
0288   bool pass;
0289   if (process22) {
0290     ++read22;
0291     pass = pt < ptsep;
0292     if (pass)
0293       ++pass22;
0294   } else {
0295     ++read23;
0296     pass = ptsep <= pt;
0297     if (pass)
0298       ++pass23;
0299   }
0300   return pass;
0301 }
0302
0303 DEFINE_FWK_MODULE(ComphepSingletopFilterPy8);