Skims/src/WZInterestingEventSelector.cc

0001 // -*- C++ -*-
0002 //
0003 //
0004 //
0005 //
0006
0007 // system include files
0008 #include <memory>
0009 #include <iostream>
0010 #include <fstream>
0011 #include <sstream>
0012 #include <string>
0013 #include <cstdio>
0014 #include <cstdlib>
0015 // user include files
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/one/EDFilter.h"
0018
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024
0025 #include "FWCore/Framework/interface/LuminosityBlock.h"
0026 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0027 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0028 #include "DataFormats/METReco/interface/CaloMETFwd.h"
0029 #include "DataFormats/METReco/interface/CaloMET.h"
0030 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0031 #include "DataFormats/METReco/interface/PFMETFwd.h"
0032 #include "DataFormats/METReco/interface/PFMET.h"
0033 #include "DataFormats/METReco/interface/PFMETCollection.h"
0034 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0035 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0036
0037 #include "TLorentzVector.h"
0038 //
0039 // class declaration
0040 //
0041
0042 using namespace reco;
0043
0044 class WZInterestingEventSelector : public edm::one::EDFilter<> {
0045 public:
0046   struct event {
0047     long run;
0048     long event;
0049     long ls;
0050     int nEle;
0051     float maxPt;
0052     float maxPtEleEta;
0053     float maxPtElePhi;
0054     float invMass;
0055     float met;
0056     float metPhi;
0057   };
0058
0059   explicit WZInterestingEventSelector(const edm::ParameterSet&);
0060   ~WZInterestingEventSelector() override;
0061
0062 private:
0063   bool filter(edm::Event&, const edm::EventSetup&) override;
0064   void endJob() override;
0065   bool electronSelection(const GsfElectron* eleRef, const math::XYZPoint& bspotPosition);
0066   // ----------member data ---------------------------
0067
0068   //std::vector<event> interestingEvents_;
0069
0070   //Pt cut
0071   float ptCut_;
0072   int missHitCut_;
0073
0074   //EB ID+ISO cuts
0075   float eb_trIsoCut_;
0076   float eb_ecalIsoCut_;
0077   float eb_hcalIsoCut_;
0078   float eb_hoeCut_;
0079   float eb_seeCut_;
0080
0081   //EE ID+ISO cuts
0082   float ee_trIsoCut_;
0083   float ee_ecalIsoCut_;
0084   float ee_hcalIsoCut_;
0085   float ee_hoeCut_;
0086   float ee_seeCut_;
0087
0088   //met Cut
0089   float metCut_;
0090
0091   //invMass Cut
0092   float invMassCut_;
0093
0094   edm::InputTag electronCollection_;
0095   edm::InputTag pfMetCollection_;
0096   edm::InputTag offlineBSCollection_;
0097 };
0098
0099 //
0100 // constructors and destructor
0101 //
0102 WZInterestingEventSelector::WZInterestingEventSelector(const edm::ParameterSet& iConfig) {
0103   ptCut_ = iConfig.getParameter<double>("ptCut");
0104   missHitCut_ = iConfig.getParameter<int>("missHitsCut");
0105
0106   eb_trIsoCut_ = iConfig.getParameter<double>("eb_trIsoCut");
0107   eb_ecalIsoCut_ = iConfig.getParameter<double>("eb_ecalIsoCut");
0108   eb_hcalIsoCut_ = iConfig.getParameter<double>("eb_hcalIsoCut");
0109   eb_hoeCut_ = iConfig.getParameter<double>("eb_hoeCut");
0110   eb_seeCut_ = iConfig.getParameter<double>("eb_seeCut");
0111
0112   ee_trIsoCut_ = iConfig.getParameter<double>("ee_trIsoCut");
0113   ee_ecalIsoCut_ = iConfig.getParameter<double>("ee_ecalIsoCut");
0114   ee_hcalIsoCut_ = iConfig.getParameter<double>("ee_hcalIsoCut");
0115   ee_hoeCut_ = iConfig.getParameter<double>("ee_hoeCut");
0116   ee_seeCut_ = iConfig.getParameter<double>("ee_seeCut");
0117
0118   metCut_ = iConfig.getParameter<double>("metCut");
0119   invMassCut_ = iConfig.getParameter<double>("invMassCut");
0120
0121   electronCollection_ =
0122       iConfig.getUntrackedParameter<edm::InputTag>("electronCollection", edm::InputTag("gsfElectrons"));
0123   pfMetCollection_ = iConfig.getUntrackedParameter<edm::InputTag>("pfMetCollection", edm::InputTag("pfMet"));
0124   offlineBSCollection_ =
0125       iConfig.getUntrackedParameter<edm::InputTag>("offlineBSCollection", edm::InputTag("offlineBeamSpot"));
0126 }
0127
0128 WZInterestingEventSelector::~WZInterestingEventSelector() {
0129   // do anything here that needs to be done at desctruction time
0130   // (e.g. close files, deallocate resources etc.)
0131 }
0132
0133 //
0134 // member functions
0135 //
0136
0137 // ------------ method called on each new Event  ------------
0138
0139 void WZInterestingEventSelector::endJob() {
0140   //   if (interestingEvents_.size()<1)
0141   //     return;
0142
0143   //   std::ostringstream oss;
0144   //   for (unsigned int iEvent=0;iEvent<interestingEvents_.size();++iEvent)
0145   //     {
0146   //       oss << "==================================" << std::endl;
0147   //       oss << "Run: " << interestingEvents_[iEvent].run << " Event: " <<  interestingEvents_[iEvent].event << " LS: " <<  interestingEvents_[iEvent].ls << std::endl;
0148   //       oss << "nGoodEle: " << interestingEvents_[iEvent].nEle << " maxPt " << interestingEvents_[iEvent].maxPt <<  " maxPtEta " << interestingEvents_[iEvent].maxPtEleEta << " maxPtPhi " << interestingEvents_[iEvent].maxPtElePhi << std::endl;
0149   //       oss << "invMass " << interestingEvents_[iEvent].invMass << " met " << interestingEvents_[iEvent].met << " metPhi " << interestingEvents_[iEvent].metPhi << std::endl;
0150   //     }
0151   //   std::string mailText;
0152   //   mailText = oss.str();
0153
0154   //   std::ofstream outputTxt;
0155   //   outputTxt.open("interestingEvents.txt");
0156   //   outputTxt << mailText;
0157   //   outputTxt.close();
0158
0159   //Sending email
0160   //   std::ostringstream subject;
0161   //   subject << "Interesting events in Run#" << interestingEvents_[0].run;
0162
0163   //   std::ostringstream command;
0164   //   command << "cat interestingEvents.txt | mail -s \"" << subject.str() << "\" Paolo.Meridiani@cern.ch";
0165
0166   //   std::string commandStr = command.str();
0167   //   char* pch = (char*)malloc( sizeof( char ) *(commandStr.length() +1) );
0168   //   string::traits_type::copy( pch, commandStr.c_str(), commandStr.length() +1 );
0169   //   int i=system(pch);
0170 }
0171
0172 bool WZInterestingEventSelector::electronSelection(const GsfElectron* eleRef, const math::XYZPoint& bspotPosition) {
0173   //   if (eleRef->trackerDrivenSeed() && !eleRef->ecalDrivenSeed())
0174   //     return false;
0175
0176   //   if (eleRef->ecalDrivenSeed())
0177   //     {
0178
0179   if (eleRef->pt() < ptCut_)
0180     return false;
0181
0182   if (eleRef->isEB()) {
0183     if (eleRef->dr03TkSumPt() / eleRef->pt() > eb_trIsoCut_)
0184       return false;
0185     if (eleRef->dr03EcalRecHitSumEt() / eleRef->pt() > eb_ecalIsoCut_)
0186       return false;
0187     if (eleRef->dr03HcalTowerSumEt() / eleRef->pt() > eb_hcalIsoCut_)
0188       return false;
0189     if (eleRef->sigmaIetaIeta() > eb_seeCut_)
0190       return false;
0191     if (eleRef->hcalOverEcal() > eb_hoeCut_)
0192       return false;
0193   } else if (eleRef->isEE()) {
0194     if (eleRef->dr03TkSumPt() / eleRef->pt() > ee_trIsoCut_)
0195       return false;
0196     if (eleRef->dr03EcalRecHitSumEt() / eleRef->pt() > ee_ecalIsoCut_)
0197       return false;
0198     if (eleRef->dr03HcalTowerSumEt() / eleRef->pt() > ee_hcalIsoCut_)
0199       return false;
0200     if (eleRef->sigmaIetaIeta() > ee_seeCut_)
0201       return false;
0202     if (eleRef->hcalOverEcal() > ee_hoeCut_)
0203       return false;
0204   }
0205
0206   if (eleRef->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) > missHitCut_)
0207     return false;
0208
0209   return true;
0210 }
0211
0212 bool WZInterestingEventSelector::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0213   // using namespace edm;
0214   edm::Handle<reco::GsfElectronCollection> gsfElectrons;
0215   iEvent.getByLabel(electronCollection_, gsfElectrons);
0216
0217   //   edm::Handle<reco::CaloMETCollection> caloMET;
0218   //   iEvent.getByLabel(edm::InputTag("met"), caloMET);
0219
0220   edm::Handle<reco::PFMETCollection> pfMET;
0221   iEvent.getByLabel(pfMetCollection_, pfMET);
0222
0223   edm::Handle<reco::BeamSpot> pBeamSpot;
0224   iEvent.getByLabel(offlineBSCollection_, pBeamSpot);
0225
0226   const reco::BeamSpot* bspot = pBeamSpot.product();
0227   const math::XYZPoint& bspotPosition = bspot->position();
0228
0229   std::vector<const GsfElectron*> goodElectrons;
0230   float ptMax = -999.;
0231   const GsfElectron* ptMaxEle = nullptr;
0232   for (reco::GsfElectronCollection::const_iterator myEle = gsfElectrons->begin(); myEle != gsfElectrons->end();
0233        ++myEle) {
0234     //Apply a minimal isolated electron selection
0235     if (!electronSelection(&(*myEle), bspotPosition))
0236       continue;
0237     goodElectrons.push_back(&(*myEle));
0238     if (myEle->pt() > ptMax) {
0239       ptMax = myEle->pt();
0240       ptMaxEle = &(*myEle);
0241     }
0242   }
0243
0244   float maxInv = -999.;
0245   TLorentzVector v1;
0246   if (ptMaxEle)
0247     v1.SetPtEtaPhiM(ptMaxEle->pt(), ptMaxEle->eta(), ptMaxEle->phi(), 0);
0248   if (goodElectrons.size() > 1) {
0249     for (unsigned int iEle = 0; iEle < goodElectrons.size(); ++iEle)
0250       if (goodElectrons[iEle] != ptMaxEle && (goodElectrons[iEle]->charge() * ptMaxEle->charge() == -1)) {
0251         TLorentzVector v2;
0252         v2.SetPtEtaPhiM(goodElectrons[iEle]->pt(), goodElectrons[iEle]->eta(), goodElectrons[iEle]->phi(), 0.);
0253         if ((v1 + v2).M() > maxInv)
0254           maxInv = (v1 + v2).M();
0255       }
0256   }
0257
0258   //Z filt: Retain event if more then 1 good ele and invMass above threshold (zee)
0259   if (goodElectrons.size() > 1 && maxInv > invMassCut_) {
0260     //interestingEvents_.push_back(thisEvent);
0261     return true;
0262   }
0263
0264   //W filt: Retain event also event with at least 1 good ele and some met
0265   if (!goodElectrons.empty() && (pfMET->begin()->et() > metCut_)) {
0266     //interestingEvents_.push_back(thisEvent);
0267     return true;
0268   }
0269
0270   return false;
0271 }
0272
0273 //define this as a plug-in
0274 DEFINE_FWK_MODULE(WZInterestingEventSelector);