GenFilters/plugins/PythiaFilterEMJetHeep.cc

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 ptMaxHadron = 0;
0204     float setCharged = 0;
0205     unsigned int Ncharged = 0;
0206
0207     for (HepMC::GenEvent::particle_const_iterator pp = myGenEvent->particles_begin(); pp != myGenEvent->particles_end();
0208          ++pp) {
0209       if ((*pp)->status() != 1)
0210         continue;  // select just stable particles
0211       int pid = (*pp)->pdg_id();
0212       int apid = std::abs(pid);
0213       if (apid == 310)
0214         apid = 22;
0215       if (apid > 11 && apid < 20)
0216         continue;  //get rid of muons and neutrinos
0217
0218       double pt_ = (*pp)->momentum().perp();
0219       double eta_ = (*pp)->momentum().eta();
0220       double phi_ = (*pp)->momentum().phi();
0221
0222       float dr = deltaR(eta_, phi_, eta, phi);
0223       if (dr <= cone_iso)
0224         setCone_iso += pt_;
0225
0226       bool charged = false;
0227       if (apid == 211 || apid == 321 || apid == 2212 || apid == 11)
0228         charged = true;
0229       if (dr <= cone_clust) {
0230         setCone_clust += pt_;
0231         if (apid == 22 || apid == 11)
0232           setEM += pt_;
0233         if (apid > 100 && pt_ > ptMaxHadron)
0234           ptMaxHadron = pt_;
0235         if (charged && pt_ > 1) {
0236           Ncharged++;
0237           setCharged += pt_;
0238         }
0239       }
0240     }
0241     setCone_iso -= setCone_clust;
0242
0243     if (pt / setCone_clust < 0.15)
0244       continue;
0245
0246     // Selection #3: min. pT of all particles in the proto-cluster
0247     if (EB && setCone_clust < fracConePtMin_EB * minEventPt)
0248       continue;
0249     if (EE && setCone_clust < fracConePtMin_EE * minEventPt)
0250       continue;
0251     if (debug)
0252       std::cout << " Selection 3 passed " << std::endl;
0253
0254     // Selections #4: min/max pT fractions of e/gamma in the proto-cluster from the total pT
0255     if (EB && setEM / setCone_clust < fracEmPtMin_EB)
0256       continue;
0257     if (EE && setEM / setCone_clust < fracEmPtMin_EE)
0258       continue;
0259     if (debug)
0260       std::cout << " Selection 4 passed " << std::endl;
0261
0262     // Selection 5: max. track pT fractions and number of tracks in the proto-cluster
0263     if ((EB && setCharged / setCone_clust > fracTrkPtMax_EB) || Ncharged > ntrkMax_EB)
0264       continue;
0265     if ((EE && setCharged / setCone_clust > fracTrkPtMax_EE) || Ncharged > ntrkMax_EE)
0266       continue;
0267     if (debug)
0268       std::cout << " Selection 5 passed " << std::endl;
0269
0270     // Selection #6: max. pT of the hadron in the proto-cluster
0271     if (EB && ptMaxHadron > ptHdMax_EB)
0272       continue;
0273     if (EE && ptMaxHadron > ptHdMax_EE)
0274       continue;
0275     if (debug)
0276       std::cout << " Selection 6 passed " << std::endl;
0277
0278     // Selection #7: max. pT fraction in the dR=cone_iso around the proto-cluster
0279     if (EB && setCone_iso / setCone_clust > isoConeMax_EB)
0280       continue;
0281     if (EE && setCone_iso / setCone_clust > isoConeMax_EE)
0282       continue;
0283     if (debug)
0284       std::cout << " Selection 7 passed " << std::endl;
0285
0286     if (debug) {
0287       std::cout << "(*itSeed)->pdg_id() = " << (*itSeed)->pdg_id() << std::endl;
0288       std::cout << "pt = " << pt << std::endl;
0289       std::cout << "setCone_clust = " << setCone_clust << std::endl;
0290       std::cout << "setEM = " << setEM << std::endl;
0291       std::cout << "ptMaxHadron = " << ptMaxHadron << std::endl;
0292       std::cout << "setCone_iso = " << setCone_iso << std::endl;
0293       std::cout << "setCone_iso / setCone_clust = " << setCone_iso / setCone_clust << std::endl;
0294     }
0295
0296     // Selection #8: any two objects should be separated by dr > drMin
0297     bool same = false;
0298     std::list<const HepMC::GenParticle*>::iterator itPart;
0299     for (itPart = candidates.begin(); itPart != candidates.end(); ++itPart) {
0300       float eta_ = (*itPart)->momentum().eta();
0301       float phi_ = (*itPart)->momentum().phi();
0302       float dr = deltaR(eta_, phi_, eta, phi);
0303       if (dr < drMin)
0304         same = true;
0305     }
0306     if (same)
0307       continue;
0308     if (debug)
0309       std::cout << " Selection 8 passed " << std::endl;
0310
0311     candidates.push_back(*itSeed);
0312   }
0313
0314   if (candidates.size() >= nPartMin) {
0315     accepted = true;
0316   }
0317
0318   if (debug)
0319     std::cout << " Proccess ID " << myGenEvent->signal_process_id() << std::endl;
0320
0321   return accepted;
0322 }
0323
0324 DEFINE_FWK_MODULE(PythiaFilterEMJetHeep);