ElectronTools/test/ElectronIsoAnalyzer.cc

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/one/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::one::EDAnalyzer<edm::one::SharedResources> {
0046 public:
0047   explicit ElectronIsoAnalyzer(const edm::ParameterSet&);
0048   ~ElectronIsoAnalyzer() override;
0049   //
0050
0051 private:
0052   void beginJob() override;
0053   void analyze(const edm::Event&, const edm::EventSetup&) override;
0054   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   usesResource(TFileService::kSharedResource);
0105
0106   verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0107   tokenGsfElectrons_ = consumes<reco::GsfElectronCollection>(iConfig.getParameter<edm::InputTag>("Electrons"));
0108   tokensIsoValElectrons_ =
0109       edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag> >("IsoValElectrons"),
0110                             [this](edm::InputTag const& tag) { return consumes<edm::ValueMap<double> >(tag); });
0111   deltaR_ = iConfig.getParameter<std::string>("deltaR");
0112   std::string eaTarget = iConfig.getParameter<std::string>("effectiveAreaTarget");
0113   rhoIsoToken_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rhoIsoInputTag"));
0114
0115   if (eaTarget == "NoCorr")
0116     effAreaTarget_ = ElectronEffectiveArea::kEleEANoCorr;
0117   else if (eaTarget == "Data2011")
0118     effAreaTarget_ = ElectronEffectiveArea::kEleEAData2011;  // default for HZZ
0119   else if (eaTarget == "Data2012")
0120     effAreaTarget_ = ElectronEffectiveArea::kEleEAData2012;  // default for HWW
0121   else if (eaTarget == "Summer11MC")
0122     effAreaTarget_ = ElectronEffectiveArea::kEleEASummer11MC;
0123   else if (eaTarget == "Fall11MC")
0124     effAreaTarget_ = ElectronEffectiveArea::kEleEAFall11MC;
0125   else
0126     throw cms::Exception("Configuration") << "Unknown effective area " << eaTarget << "\n";
0127
0128   if (deltaR_ == "03") {
0129     effAreaGammaPlusNeutralHad_ = ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03;
0130   } else if (deltaR_ == "04") {
0131     effAreaGammaPlusNeutralHad_ = ElectronEffectiveArea::kEleGammaAndNeutralHadronIso04;
0132   } else
0133     throw cms::Exception("Configuration") << "Unsupported deltaR " << deltaR_ << "\n";
0134
0135   edm::Service<TFileService> fs;
0136   chargedBarrel_ = fs->make<TH1F>("chargedBarrel", ";Sum pT/pT", 100, 0, 4);
0137   photonBarrel_ = fs->make<TH1F>("photonBarrel", ";Sum pT/pT", 100, 0, 4);
0138   neutralBarrel_ = fs->make<TH1F>("neutralBarrel", ";Sum pT/pT", 100, 0, 4);
0139
0140   chargedEndcaps_ = fs->make<TH1F>("chargedEndcaps", ";Sum pT/pT", 100, 0, 4);
0141   photonEndcaps_ = fs->make<TH1F>("photonEndcaps", ";Sum pT/pT", 100, 0, 4);
0142   neutralEndcaps_ = fs->make<TH1F>("neutralEndcaps", ";Sum pT/pT", 100, 0, 4);
0143
0144   sumBarrel_ = fs->make<TH1F>("allbarrel", ";Sum pT/pT", 100, 0, 4);
0145   sumEndcaps_ = fs->make<TH1F>("allendcaps", ";Sum pT/pT", 100, 0, 4);
0146
0147   sumCorrBarrel_ = fs->make<TH1F>("allcorrbarrel", ";Sum pT/pT", 100, 0, 4);
0148   sumCorrEndcaps_ = fs->make<TH1F>("allcorrendcaps", ";Sum pT/pT", 100, 0, 4);
0149
0150   missHitsBarrel_ = fs->make<TH1F>("missHitsBarrel_", "", 10, 0, 10);
0151   missHitsEndcap_ = fs->make<TH1F>("missHitsEndcap_", "", 10, 0, 10);
0152 }
0153
0154 ElectronIsoAnalyzer::~ElectronIsoAnalyzer() {
0155   // do anything here that needs to be done at desctruction time
0156   // (e.g. close files, deallocate resources etc.)
0157 }
0158
0159 //
0160 // member functions
0161 //
0162
0163 // ------------ method called to for each event  ------------
0164 void ElectronIsoAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0165   edm::Handle<reco::GsfElectronCollection> theEGammaCollection;
0166   iEvent.getByToken(tokenGsfElectrons_, theEGammaCollection);
0167   const reco::GsfElectronCollection theEGamma = *(theEGammaCollection.product());
0168
0169   // rho for isolation
0170   edm::Handle<double> rhoIso_h;
0171   iEvent.getByToken(rhoIsoToken_, rhoIso_h);
0172   double rhoIso = *(rhoIso_h.product());
0173
0174   unsigned nTypes = 3;
0175   IsoValues electronIsoValues(nTypes);
0176
0177   for (size_t j = 0; j < tokensIsoValElectrons_.size(); ++j) {
0178     iEvent.getByToken(tokensIsoValElectrons_[j], electronIsoValues[j]);
0179   }
0180
0181   unsigned nele = theEGammaCollection->size();
0182
0183   for (unsigned iele = 0; iele < nele; ++iele) {
0184     reco::GsfElectronRef myElectronRef(theEGammaCollection, iele);
0185
0186     const IsoValues* myIsoValues = &electronIsoValues;
0187
0188     double charged = (*(*myIsoValues)[0])[myElectronRef];
0189     double photon = (*(*myIsoValues)[1])[myElectronRef];
0190     double neutral = (*(*myIsoValues)[2])[myElectronRef];
0191
0192     float abseta = fabs(myElectronRef->superCluster()->eta());
0193
0194     float eff_area_phnh =
0195         ElectronEffectiveArea::GetElectronEffectiveArea(effAreaGammaPlusNeutralHad_, abseta, effAreaTarget_);
0196
0197     double myRho = std::max<double>(0., rhoIso);
0198
0199     float myPfIsoPuCorr = charged + std::max<float>(0.f, (photon + neutral) - eff_area_phnh * myRho);
0200
0201     if (verbose_) {
0202       std::cout << " run " << iEvent.id().run() << " lumi " << iEvent.id().luminosityBlock() << " event "
0203                 << iEvent.id().event();
0204       std::cout << " pt " << myElectronRef->pt() << " eta " << myElectronRef->eta() << " phi " << myElectronRef->phi()
0205                 << " charge " << myElectronRef->charge() << " : " << std::endl;
0206       ;
0207
0208       // print values also from alternate code
0209       std::cout << " ChargedIso " << charged << std::endl;
0210       std::cout << " PhotonIso " << photon << std::endl;
0211       std::cout << " NeutralHadron Iso " << neutral << std::endl;
0212     }
0213     if (myElectronRef->isEB()) {
0214       chargedBarrel_->Fill(charged / myElectronRef->pt());
0215       photonBarrel_->Fill(photon / myElectronRef->pt());
0216       neutralBarrel_->Fill(neutral / myElectronRef->pt());
0217       sumBarrel_->Fill((charged + photon + neutral) / myElectronRef->pt());
0218       sumCorrBarrel_->Fill(myPfIsoPuCorr / myElectronRef->pt());
0219       missHitsBarrel_->Fill(
0220           myElectronRef->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
0221
0222     } else {
0223       chargedEndcaps_->Fill(charged / myElectronRef->pt());
0224       photonEndcaps_->Fill(photon / myElectronRef->pt());
0225       neutralEndcaps_->Fill(neutral / myElectronRef->pt());
0226       sumEndcaps_->Fill((charged + photon + neutral) / myElectronRef->pt());
0227       sumCorrEndcaps_->Fill(myPfIsoPuCorr / myElectronRef->pt());
0228       missHitsEndcap_->Fill(
0229           myElectronRef->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS));
0230     }
0231   }
0232 }
0233 // ------------ method called once each job just before starting event loop  ------------
0234 void ElectronIsoAnalyzer::beginJob() { ev = 0; }
0235
0236 // ------------ method called once each job just after ending the event loop  ------------
0237 void ElectronIsoAnalyzer::endJob() { std::cout << " endJob:: #events " << ev << std::endl; }
0238
0239 //define this as a plug-in
0240 DEFINE_FWK_MODULE(ElectronIsoAnalyzer);