Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:29:58

0001 /** \class PythiaFilterEMJetHeep
0002  *
0003  *  PythiaFilterEMJetHeep filter implements generator-level preselections
0004  *  of events with for studying background to high-energetic electrons
0005  *
0006  * \author Dmitry Bandurin (KSU), Jeremy Werner (Princeton)
0007  *
0008  ************************************************************/
0009 
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/global/EDFilter.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/EDGetToken.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0019 
0020 #include <cmath>
0021 #include <cstdlib>
0022 #include <iostream>
0023 #include <list>
0024 #include <string>
0025 
0026 class PythiaFilterEMJetHeep : public edm::global::EDFilter<> {
0027 public:
0028   double deltaR(double eta0, double phi0, double eta, double phi) const;
0029 
0030   explicit PythiaFilterEMJetHeep(const edm::ParameterSet&);
0031 
0032   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0033   void beginJob() override;
0034 
0035 private:
0036   const edm::EDGetTokenT<edm::HepMCProduct> token_;
0037 
0038   double minEventPt;
0039   const double etaMax;
0040   const double cone_clust;
0041   const double cone_iso;
0042   const unsigned int nPartMin;
0043   const double drMin;
0044 
0045   double ptSeedMin_EB;
0046   double fracConePtMin_EB;
0047   double ptHdMax_EB;
0048   double fracEmPtMin_EB;
0049   double fracTrkPtMax_EB;
0050   unsigned int ntrkMax_EB;
0051   double isoConeMax_EB;
0052 
0053   double ptSeedMin_EE;
0054   double fracConePtMin_EE;
0055   double ptHdMax_EE;
0056   double fracEmPtMin_EE;
0057   double fracTrkPtMax_EE;
0058   unsigned int ntrkMax_EE;
0059   double isoConeMax_EE;
0060 
0061   const bool minbias;
0062 
0063   const bool debug;
0064 };
0065 
0066 double PythiaFilterEMJetHeep::deltaR(double eta0, double phi0, double eta, double phi) const {
0067   double dphi = phi - phi0;
0068   if (dphi > M_PI)
0069     dphi -= 2 * M_PI;
0070   else if (dphi <= -M_PI)
0071     dphi += 2 * M_PI;
0072   return std::sqrt(dphi * dphi + (eta - eta0) * (eta - eta0));
0073 }
0074 
0075 namespace {
0076   struct ParticlePtGreater {
0077     bool operator()(const HepMC::GenParticle* a, const HepMC::GenParticle* b) const {
0078       return a->momentum().perp() > b->momentum().perp();
0079     }
0080   };
0081 }  // namespace
0082 
0083 PythiaFilterEMJetHeep::PythiaFilterEMJetHeep(const edm::ParameterSet& iConfig)
0084     : token_(consumes<edm::HepMCProduct>(
0085           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0086       minEventPt(iConfig.getUntrackedParameter<double>("MinEventPt", 40.)),
0087       etaMax(iConfig.getUntrackedParameter<double>("MaxEta", 2.8)),
0088       cone_clust(iConfig.getUntrackedParameter<double>("ConeClust", 0.10)),
0089       cone_iso(iConfig.getUntrackedParameter<double>("ConeIso", 0.50)),
0090       nPartMin(iConfig.getUntrackedParameter<unsigned int>("NumPartMin", 2)),
0091       drMin(iConfig.getUntrackedParameter<double>("dRMin", 0.4)),
0092       minbias(iConfig.getUntrackedParameter<bool>("Minbias", false)),
0093       debug(iConfig.getUntrackedParameter<bool>("Debug", true)) {}
0094 
0095 void PythiaFilterEMJetHeep::beginJob() {
0096   // parametarizations of presel. criteria:
0097 
0098   ptSeedMin_EB = 6.0 + (minEventPt - 80.) * 0.035;
0099   ptSeedMin_EE = 5.5 + (minEventPt - 80.) * 0.033;
0100   fracConePtMin_EB = 0.60 - (minEventPt - 80.) * 0.0009;
0101   fracConePtMin_EE = fracConePtMin_EB;
0102   fracEmPtMin_EB = 0.30 + (minEventPt - 80.) * 0.0017;
0103   if (minEventPt >= 225)
0104     fracEmPtMin_EB = 0.40 + (minEventPt - 230.) * 0.00063;
0105   fracEmPtMin_EE = 0.30 + (minEventPt - 80.) * 0.0033;
0106   if (minEventPt >= 165)
0107     fracEmPtMin_EE = 0.55 + (minEventPt - 170.) * 0.0005;
0108   fracTrkPtMax_EB = 0.80 - (minEventPt - 80.) * 0.001;
0109   fracTrkPtMax_EE = 0.70 - (minEventPt - 80.) * 0.001;
0110   isoConeMax_EB = 0.35;
0111   isoConeMax_EE = 0.40;
0112   ptHdMax_EB = 40;
0113   ptHdMax_EB = 45;
0114   ntrkMax_EB = 4;
0115   ntrkMax_EE = 4;
0116 
0117   if (minbias) {
0118     minEventPt = 1.0;
0119     ptSeedMin_EB = 1.5;
0120     ptSeedMin_EE = 1.5;
0121     fracConePtMin_EB = 0.20;
0122     fracConePtMin_EE = fracConePtMin_EB;
0123     fracEmPtMin_EB = 0.20;
0124     fracEmPtMin_EE = 0.20;
0125     fracTrkPtMax_EB = 0.80;
0126     fracTrkPtMax_EE = 0.80;
0127     isoConeMax_EB = 1.0;
0128     isoConeMax_EE = 1.0;
0129     ptHdMax_EB = 7;
0130     ptHdMax_EB = 7;
0131     ntrkMax_EB = 2;
0132     ntrkMax_EE = 2;
0133   }
0134 }
0135 
0136 bool PythiaFilterEMJetHeep::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0137   bool accepted = false;
0138   edm::Handle<edm::HepMCProduct> evt;
0139   iEvent.getByToken(token_, evt);
0140 
0141   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0142 
0143   std::list<const HepMC::GenParticle*> seeds;
0144   std::list<const HepMC::GenParticle*> candidates;
0145 
0146   // collect the seeds
0147   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0148        ++p) {
0149     double pt = (*p)->momentum().perp();
0150     double eta = (*p)->momentum().eta();
0151 
0152     int pdgid = (*p)->pdg_id();
0153     if (!(pdgid == 22 || std::abs(pdgid) == 11 || std::abs(pdgid) == 211 || std::abs(pdgid) == 321))
0154       continue;
0155     //if ( !(pdgid==22 || std::abs(pdgid)==11) ) continue;
0156 
0157     //  Selection #1: seed particle ETA
0158     if (std::abs(eta) > etaMax)
0159       continue;
0160     bool EB = false;
0161     bool EE = false;
0162     if (std::fabs(eta) < 1.5)
0163       EB = true;
0164     else if (std::fabs(eta) >= 1.5 && std::fabs(eta) < etaMax)
0165       EE = true;
0166     if (debug)
0167       std::cout << " Selection 1 passed " << std::endl;
0168 
0169     //  Selection #2: seed particle pT
0170     if (EB && pt < ptSeedMin_EB)
0171       continue;
0172     if (EE && pt < ptSeedMin_EE)
0173       continue;
0174     if (debug)
0175       std::cout << " Selection 2 passed " << std::endl;
0176 
0177     seeds.push_back(*p);
0178   }
0179 
0180   if (seeds.size() < nPartMin)
0181     return false;
0182 
0183   seeds.sort(ParticlePtGreater());
0184 
0185   // select the proto-clusters
0186   std::list<const HepMC::GenParticle*>::iterator itSeed;
0187 
0188   for (itSeed = seeds.begin(); itSeed != seeds.end(); ++itSeed) {
0189     double pt = (*itSeed)->momentum().perp();
0190     double eta = (*itSeed)->momentum().eta();
0191     double phi = (*itSeed)->momentum().phi();
0192 
0193     bool EB = false;
0194     bool EE = false;
0195     if (std::fabs(eta) < 1.5)
0196       EB = true;
0197     else if (std::fabs(eta) >= 1.5 && std::fabs(eta) < etaMax)
0198       EE = true;
0199 
0200     float setCone_iso = 0;
0201     float setCone_clust = 0;
0202     float setEM = 0;
0203     float setHAD = 0;
0204     float ptMaxHadron = 0;
0205     float setCharged = 0;
0206     unsigned int Ncharged = 0;
0207 
0208     for (HepMC::GenEvent::particle_const_iterator pp = myGenEvent->particles_begin(); pp != myGenEvent->particles_end();
0209          ++pp) {
0210       if ((*pp)->status() != 1)
0211         continue;  // select just stable particles
0212       int pid = (*pp)->pdg_id();
0213       int apid = std::abs(pid);
0214       if (apid == 310)
0215         apid = 22;
0216       if (apid > 11 && apid < 20)
0217         continue;  //get rid of muons and neutrinos
0218 
0219       double pt_ = (*pp)->momentum().perp();
0220       double eta_ = (*pp)->momentum().eta();
0221       double phi_ = (*pp)->momentum().phi();
0222 
0223       float dr = deltaR(eta_, phi_, eta, phi);
0224       if (dr <= cone_iso)
0225         setCone_iso += pt_;
0226 
0227       bool charged = false;
0228       if (apid == 211 || apid == 321 || apid == 2212 || apid == 11)
0229         charged = true;
0230       if (dr <= cone_clust) {
0231         setCone_clust += pt_;
0232         if (apid == 22 || apid == 11)
0233           setEM += pt_;
0234         if (apid > 100)
0235           setHAD += pt_;
0236         if (apid > 100 && pt_ > ptMaxHadron)
0237           ptMaxHadron = pt_;
0238         if (charged && pt_ > 1) {
0239           Ncharged++;
0240           setCharged += pt_;
0241         }
0242       }
0243     }
0244     setCone_iso -= setCone_clust;
0245 
0246     if (pt / setCone_clust < 0.15)
0247       continue;
0248 
0249     // Selection #3: min. pT of all particles in the proto-cluster
0250     if (EB && setCone_clust < fracConePtMin_EB * minEventPt)
0251       continue;
0252     if (EE && setCone_clust < fracConePtMin_EE * minEventPt)
0253       continue;
0254     if (debug)
0255       std::cout << " Selection 3 passed " << std::endl;
0256 
0257     // Selections #4: min/max pT fractions of e/gamma in the proto-cluster from the total pT
0258     if (EB && setEM / setCone_clust < fracEmPtMin_EB)
0259       continue;
0260     if (EE && setEM / setCone_clust < fracEmPtMin_EE)
0261       continue;
0262     if (debug)
0263       std::cout << " Selection 4 passed " << std::endl;
0264 
0265     // Selection 5: max. track pT fractions and number of tracks in the proto-cluster
0266     if ((EB && setCharged / setCone_clust > fracTrkPtMax_EB) || Ncharged > ntrkMax_EB)
0267       continue;
0268     if ((EE && setCharged / setCone_clust > fracTrkPtMax_EE) || Ncharged > ntrkMax_EE)
0269       continue;
0270     if (debug)
0271       std::cout << " Selection 5 passed " << std::endl;
0272 
0273     // Selection #6: max. pT of the hadron in the proto-cluster
0274     if (EB && ptMaxHadron > ptHdMax_EB)
0275       continue;
0276     if (EE && ptMaxHadron > ptHdMax_EE)
0277       continue;
0278     if (debug)
0279       std::cout << " Selection 6 passed " << std::endl;
0280 
0281     // Selection #7: max. pT fraction in the dR=cone_iso around the proto-cluster
0282     if (EB && setCone_iso / setCone_clust > isoConeMax_EB)
0283       continue;
0284     if (EE && setCone_iso / setCone_clust > isoConeMax_EE)
0285       continue;
0286     if (debug)
0287       std::cout << " Selection 7 passed " << std::endl;
0288 
0289     if (debug) {
0290       std::cout << "(*itSeed)->pdg_id() = " << (*itSeed)->pdg_id() << std::endl;
0291       std::cout << "pt = " << pt << std::endl;
0292       std::cout << "setCone_clust = " << setCone_clust << std::endl;
0293       std::cout << "setEM = " << setEM << std::endl;
0294       std::cout << "ptMaxHadron = " << ptMaxHadron << std::endl;
0295       std::cout << "setCone_iso = " << setCone_iso << std::endl;
0296       std::cout << "setCone_iso / setCone_clust = " << setCone_iso / setCone_clust << std::endl;
0297     }
0298 
0299     // Selection #8: any two objects should be separated by dr > drMin
0300     bool same = false;
0301     std::list<const HepMC::GenParticle*>::iterator itPart;
0302     for (itPart = candidates.begin(); itPart != candidates.end(); ++itPart) {
0303       float eta_ = (*itPart)->momentum().eta();
0304       float phi_ = (*itPart)->momentum().phi();
0305       float dr = deltaR(eta_, phi_, eta, phi);
0306       if (dr < drMin)
0307         same = true;
0308     }
0309     if (same)
0310       continue;
0311     if (debug)
0312       std::cout << " Selection 8 passed " << std::endl;
0313 
0314     candidates.push_back(*itSeed);
0315   }
0316 
0317   if (candidates.size() >= nPartMin) {
0318     accepted = true;
0319   }
0320 
0321   if (debug)
0322     std::cout << " Proccess ID " << myGenEvent->signal_process_id() << std::endl;
0323 
0324   return accepted;
0325 }
0326 
0327 DEFINE_FWK_MODULE(PythiaFilterEMJetHeep);