Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:14:03

0001 // -*- C++ -*-
0002 //
0003 // Package:    ElectronIsoAnalyzer
0004 // Class:      ElectronIsoAnalyzer
0005 //
0006 /**\class ElectronIsoAnalyzer
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Daniele Benedetti
0015 
0016 // system include files
0017 #include <memory>
0018 
0019 // user include files
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/EDAnalyzer.h"
0022 
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Utilities/interface/transform.h"
0026 
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/ServiceRegistry/interface/Service.h"
0029 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0030 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0031 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0032 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0033 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0034 #include "DataFormats/Common/interface/ValueMap.h"
0035 #include "EgammaAnalysis/ElectronTools/interface/ElectronEffectiveArea.h"
0036 
0037 #include <iostream>
0038 #include <string>
0039 #include <map>
0040 
0041 //
0042 // class decleration
0043 //
0044 
0045 class ElectronIsoAnalyzer : public edm::EDAnalyzer {
0046 public:
0047   explicit ElectronIsoAnalyzer(const edm::ParameterSet&);
0048   ~ElectronIsoAnalyzer();
0049   //
0050 
0051 private:
0052   virtual void beginJob() override;
0053   virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
0054   virtual void endJob() override;
0055 
0056   edm::ParameterSet conf_;
0057 
0058   unsigned int ev;
0059   // ----------member data ---------------------------
0060   bool verbose_;
0061   edm::EDGetTokenT<reco::GsfElectronCollection> tokenGsfElectrons_;
0062   edm::EDGetTokenT<double> rhoIsoToken_;
0063   std::vector<edm::EDGetTokenT<edm::ValueMap<double> > > tokensIsoValElectrons_;
0064   typedef std::vector<edm::Handle<edm::ValueMap<double> > > IsoValues;
0065   ElectronEffectiveArea::ElectronEffectiveAreaTarget effAreaTarget_;
0066   ElectronEffectiveArea::ElectronEffectiveAreaType effAreaGammaPlusNeutralHad_;
0067   std::string rho_;
0068   std::string deltaR_;
0069 
0070   // Control histos
0071   TH1F* chargedBarrel_;
0072   TH1F* photonBarrel_;
0073   TH1F* neutralBarrel_;
0074 
0075   TH1F* chargedEndcaps_;
0076   TH1F* photonEndcaps_;
0077   TH1F* neutralEndcaps_;
0078 
0079   TH1F* sumBarrel_;
0080   TH1F* sumEndcaps_;
0081 
0082   TH1F* missHitsBarrel_;
0083   TH1F* missHitsEndcap_;
0084 
0085   TH1F* sumCorrBarrel_;
0086   TH1F* sumCorrEndcaps_;
0087 };
0088 
0089 //
0090 // constants, enums and typedefs
0091 //
0092 
0093 //
0094 // static data member definitions
0095 //
0096 
0097 //
0098 // constructors and destructor
0099 //
0100 ElectronIsoAnalyzer::ElectronIsoAnalyzer(const edm::ParameterSet& iConfig)
0101     : conf_(iConfig)
0102 
0103 {
0104   verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0105   tokenGsfElectrons_ = consumes<reco::GsfElectronCollection>(iConfig.getParameter<edm::InputTag>("Electrons"));
0106   tokensIsoValElectrons_ =
0107       edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag> >("IsoValElectrons"),
0108                             [this](edm::InputTag const& tag) { return consumes<edm::ValueMap<double> >(tag); });
0109   deltaR_ = iConfig.getParameter<std::string>("deltaR");
0110   std::string eaTarget = iConfig.getParameter<std::string>("effectiveAreaTarget");
0111   rhoIsoToken_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rhoIsoInputTag"));
0112 
0113   if (eaTarget == "NoCorr")
0114     effAreaTarget_ = ElectronEffectiveArea::kEleEANoCorr;
0115   else if (eaTarget == "Data2011")
0116     effAreaTarget_ = ElectronEffectiveArea::kEleEAData2011;  // default for HZZ
0117   else if (eaTarget == "Data2012")
0118     effAreaTarget_ = ElectronEffectiveArea::kEleEAData2012;  // default for HWW
0119   else if (eaTarget == "Summer11MC")
0120     effAreaTarget_ = ElectronEffectiveArea::kEleEASummer11MC;
0121   else if (eaTarget == "Fall11MC")
0122     effAreaTarget_ = ElectronEffectiveArea::kEleEAFall11MC;
0123   else
0124     throw cms::Exception("Configuration") << "Unknown effective area " << eaTarget << "\n";
0125 
0126   if (deltaR_ == "03") {
0127     effAreaGammaPlusNeutralHad_ = ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03;
0128   } else if (deltaR_ == "04") {
0129     effAreaGammaPlusNeutralHad_ = ElectronEffectiveArea::kEleGammaAndNeutralHadronIso04;
0130   } else
0131     throw cms::Exception("Configuration") << "Unsupported deltaR " << deltaR_ << "\n";
0132 
0133   edm::Service<TFileService> fs;
0134   chargedBarrel_ = fs->make<TH1F>("chargedBarrel", ";Sum pT/pT", 100, 0, 4);
0135   photonBarrel_ = fs->make<TH1F>("photonBarrel", ";Sum pT/pT", 100, 0, 4);
0136   neutralBarrel_ = fs->make<TH1F>("neutralBarrel", ";Sum pT/pT", 100, 0, 4);
0137 
0138   chargedEndcaps_ = fs->make<TH1F>("chargedEndcaps", ";Sum pT/pT", 100, 0, 4);
0139   photonEndcaps_ = fs->make<TH1F>("photonEndcaps", ";Sum pT/pT", 100, 0, 4);
0140   neutralEndcaps_ = fs->make<TH1F>("neutralEndcaps", ";Sum pT/pT", 100, 0, 4);
0141 
0142   sumBarrel_ = fs->make<TH1F>("allbarrel", ";Sum pT/pT", 100, 0, 4);
0143   sumEndcaps_ = fs->make<TH1F>("allendcaps", ";Sum pT/pT", 100, 0, 4);
0144 
0145   sumCorrBarrel_ = fs->make<TH1F>("allcorrbarrel", ";Sum pT/pT", 100, 0, 4);
0146   sumCorrEndcaps_ = fs->make<TH1F>("allcorrendcaps", ";Sum pT/pT", 100, 0, 4);
0147 
0148   missHitsBarrel_ = fs->make<TH1F>("missHitsBarrel_", "", 10, 0, 10);
0149   missHitsEndcap_ = fs->make<TH1F>("missHitsEndcap_", "", 10, 0, 10);
0150 }
0151 
0152 ElectronIsoAnalyzer::~ElectronIsoAnalyzer() {
0153   // do anything here that needs to be done at desctruction time
0154   // (e.g. close files, deallocate resources etc.)
0155 }
0156 
0157 //
0158 // member functions
0159 //
0160 
0161 // ------------ method called to for each event  ------------
0162 void ElectronIsoAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0163   edm::Handle<reco::GsfElectronCollection> theEGammaCollection;
0164   iEvent.getByToken(tokenGsfElectrons_, theEGammaCollection);
0165   const reco::GsfElectronCollection theEGamma = *(theEGammaCollection.product());
0166 
0167   // rho for isolation
0168   edm::Handle<double> rhoIso_h;
0169   iEvent.getByToken(rhoIsoToken_, rhoIso_h);
0170   double rhoIso = *(rhoIso_h.product());
0171 
0172   unsigned nTypes = 3;
0173   IsoValues electronIsoValues(nTypes);
0174 
0175   for (size_t j = 0; j < tokensIsoValElectrons_.size(); ++j) {
0176     iEvent.getByToken(tokensIsoValElectrons_[j], electronIsoValues[j]);
0177   }
0178 
0179   unsigned nele = theEGammaCollection->size();
0180 
0181   for (unsigned iele = 0; iele < nele; ++iele) {
0182     reco::GsfElectronRef myElectronRef(theEGammaCollection, iele);
0183 
0184     const IsoValues* myIsoValues = &electronIsoValues;
0185 
0186     double charged = (*(*myIsoValues)[0])[myElectronRef];
0187     double photon = (*(*myIsoValues)[1])[myElectronRef];
0188     double neutral = (*(*myIsoValues)[2])[myElectronRef];
0189 
0190     float abseta = fabs(myElectronRef->superCluster()->eta());
0191 
0192     float eff_area_phnh =
0193         ElectronEffectiveArea::GetElectronEffectiveArea(effAreaGammaPlusNeutralHad_, abseta, effAreaTarget_);
0194 
0195     double myRho = std::max<double>(0., rhoIso);
0196 
0197     float myPfIsoPuCorr = charged + std::max<float>(0.f, (photon + neutral) - eff_area_phnh * myRho);
0198 
0199     if (verbose_) {
0200       std::cout << " run " << iEvent.id().run() << " lumi " << iEvent.id().luminosityBlock() << " event "
0201                 << iEvent.id().event();
0202       std::cout << " pt " << myElectronRef->pt() << " eta " << myElectronRef->eta() << " phi " << myElectronRef->phi()
0203                 << " charge " << myElectronRef->charge() << " : " << std::endl;
0204       ;
0205 
0206       // print values also from alternate code
0207       std::cout << " ChargedIso " << charged << std::endl;
0208       std::cout << " PhotonIso " << photon << std::endl;
0209       std::cout << " NeutralHadron Iso " << neutral << std::endl;
0210     }
0211     if (myElectronRef->isEB()) {
0212       chargedBarrel_->Fill(charged / myElectronRef->pt());
0213       photonBarrel_->Fill(photon / myElectronRef->pt());
0214       neutralBarrel_->Fill(neutral / myElectronRef->pt());
0215       sumBarrel_->Fill((charged + photon + neutral) / myElectronRef->pt());
0216       sumCorrBarrel_->Fill(myPfIsoPuCorr / myElectronRef->pt());
0217       missHitsBarrel_->Fill(
0218           myElectronRef->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
0219 
0220     } else {
0221       chargedEndcaps_->Fill(charged / myElectronRef->pt());
0222       photonEndcaps_->Fill(photon / myElectronRef->pt());
0223       neutralEndcaps_->Fill(neutral / myElectronRef->pt());
0224       sumEndcaps_->Fill((charged + photon + neutral) / myElectronRef->pt());
0225       sumCorrEndcaps_->Fill(myPfIsoPuCorr / myElectronRef->pt());
0226       missHitsEndcap_->Fill(
0227           myElectronRef->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
0228     }
0229   }
0230 }
0231 // ------------ method called once each job just before starting event loop  ------------
0232 void ElectronIsoAnalyzer::beginJob() { ev = 0; }
0233 
0234 // ------------ method called once each job just after ending the event loop  ------------
0235 void ElectronIsoAnalyzer::endJob() { std::cout << " endJob:: #events " << ev << std::endl; }
0236 
0237 //define this as a plug-in
0238 DEFINE_FWK_MODULE(ElectronIsoAnalyzer);