ParticleFlow/test/PFIsoReaderDemo.cc

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/transform.h"
0007 #include "DataFormats/Common/interface/ValueMap.h"
0008 #include "CommonTools/ParticleFlow/test/PFIsoReaderDemo.h"
0009 #include "DataFormats/Math/interface/deltaPhi.h"
0010 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0011 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0012 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0013 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0014
0015 PFIsoReaderDemo::PFIsoReaderDemo(const edm::ParameterSet& iConfig) {
0016   usesResource(TFileService::kSharedResource);
0017   inputTagGsfElectrons_ = iConfig.getParameter<edm::InputTag>("Electrons");
0018   tokenGsfElectrons_ = consumes<reco::GsfElectronCollection>(inputTagGsfElectrons_);
0019   inputTagPhotons_ = iConfig.getParameter<edm::InputTag>("Photons");
0020   tokenPhotons_ = consumes<reco::PhotonCollection>(inputTagPhotons_);
0021
0022   //   // not needed at the moment
0023   //   inputTagPFCandidateMap_ = iConfig.getParameter< edm::InputTag>("PFCandidateMap");
0024
0025   //   inputTagIsoDepElectrons_ = iConfig.getParameter< std::vector<edm::InputTag> >("IsoDepElectron");
0026   tokensIsoDepElectrons_ = edm::vector_transform(
0027       iConfig.getParameter<std::vector<edm::InputTag> >("IsoDepElectron"),
0028       [this](edm::InputTag const& tag) { return consumes<edm::ValueMap<reco::IsoDeposit> >(tag); });
0029   //   inputTagIsoDepPhotons_ = iConfig.getParameter< std::vector<edm::InputTag> >("IsoDepPhoton");
0030   tokensIsoDepPhotons_ = edm::vector_transform(
0031       iConfig.getParameter<std::vector<edm::InputTag> >("IsoDepPhoton"),
0032       [this](edm::InputTag const& tag) { return consumes<edm::ValueMap<reco::IsoDeposit> >(tag); });
0033   // No longer needed. e/g recommendation (04/04/12)
0034   //  inputTagIsoValElectronsNoPFId_ = iConfig.getParameter< std::vector<edm::InputTag> >("IsoValElectronNoPF");
0035   //   inputTagIsoValElectronsPFId_   = iConfig.getParameter< std::vector<edm::InputTag> >("IsoValElectronPF");
0036   tokensIsoValElectronsPFId_ =
0037       edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag> >("IsoValElectronPF"),
0038                             [this](edm::InputTag const& tag) { return consumes<edm::ValueMap<double> >(tag); });
0039   //   inputTagIsoValPhotonsPFId_   = iConfig.getParameter< std::vector<edm::InputTag> >("IsoValPhoton");
0040   tokensIsoValPhotonsPFId_ =
0041       edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag> >("IsoValPhoton"),
0042                             [this](edm::InputTag const& tag) { return consumes<edm::ValueMap<double> >(tag); });
0043
0044   printElectrons_ = iConfig.getParameter<bool>("PrintElectrons");
0045   printPhotons_ = iConfig.getParameter<bool>("PrintPhotons");
0046
0047   // Control plots
0048   TFileDirectory dir = fileservice_->mkdir("PF ISO");
0049   chargedBarrelElectrons_ = dir.make<TH1F>("chargedBarrelElectrons", ";Sum pT/pT", 100, 0, 4);
0050   photonBarrelElectrons_ = dir.make<TH1F>("photonBarrelElectrons", ";Sum pT/pT", 100, 0, 4);
0051   neutralBarrelElectrons_ = dir.make<TH1F>("neutralBarrelElectrons", ";Sum pT/pT", 100, 0, 4);
0052
0053   chargedEndcapsElectrons_ = dir.make<TH1F>("chargedEndcapsElectrons", ";Sum pT/pT", 100, 0, 4);
0054   photonEndcapsElectrons_ = dir.make<TH1F>("photonEndcapsElectrons", ";Sum pT/pT", 100, 0, 4);
0055   neutralEndcapsElectrons_ = dir.make<TH1F>("neutralEndcapsElectrons", ";Sum pT/pT", 100, 0, 4);
0056
0057   sumBarrelElectrons_ = dir.make<TH1F>("allbarrelElectrons", ";Sum pT/pT", 100, 0, 4);
0058   sumEndcapsElectrons_ = dir.make<TH1F>("allendcapsElectrons", ";Sum pT/pT", 100, 0, 4);
0059
0060   chargedBarrelPhotons_ = dir.make<TH1F>("chargedBarrelPhotons", ";Sum pT/pT", 100, 0, 4);
0061   photonBarrelPhotons_ = dir.make<TH1F>("photonBarrelPhotons", ";Sum pT/pT", 100, 0, 4);
0062   neutralBarrelPhotons_ = dir.make<TH1F>("neutralBarrelPhotons", ";Sum pT/pT", 100, 0, 4);
0063
0064   chargedEndcapsPhotons_ = dir.make<TH1F>("chargedEndcapsPhotons", ";Sum pT/pT", 100, 0, 4);
0065   photonEndcapsPhotons_ = dir.make<TH1F>("photonEndcapsPhotons", ";Sum pT/pT", 100, 0, 4);
0066   neutralEndcapsPhotons_ = dir.make<TH1F>("neutralEndcapsPhotons", ";Sum pT/pT", 100, 0, 4);
0067
0068   sumBarrelPhotons_ = dir.make<TH1F>("allbarrelPhotons", ";Sum pT/pT", 100, 0, 4);
0069   sumEndcapsPhotons_ = dir.make<TH1F>("allendcapsPhotons", ";Sum pT/pT", 100, 0, 4);
0070 }
0071
0072 PFIsoReaderDemo::~PFIsoReaderDemo() { ; }
0073
0074 void PFIsoReaderDemo::beginRun(edm::Run const&, edm::EventSetup const&) { ; }
0075
0076 void PFIsoReaderDemo::analyze(const edm::Event& iEvent, const edm::EventSetup& c) {
0077   edm::Handle<reco::GsfElectronCollection> gsfElectronH;
0078   bool found = iEvent.getByToken(tokenGsfElectrons_, gsfElectronH);
0079   if (!found) {
0080     std::ostringstream err;
0081     err << " cannot get GsfElectrons: " << inputTagGsfElectrons_ << std::endl;
0082     edm::LogError("PFIsoReaderDemo") << err.str();
0083     throw cms::Exception("MissingProduct", err.str());
0084   }
0085
0086   edm::Handle<reco::PhotonCollection> photonH;
0087   found = iEvent.getByToken(tokenPhotons_, photonH);
0088   if (!found) {
0089     std::ostringstream err;
0090     err << " cannot get Photons: " << inputTagPhotons_ << std::endl;
0091     edm::LogError("PFIsoReaderDemo") << err.str();
0092     throw cms::Exception("MissingProduct", err.str());
0093   }
0094
0095   // get the iso deposits. 3 (charged hadrons, photons, neutral hadrons)
0096   unsigned nTypes = 3;
0097   IsoDepositMaps electronIsoDep(nTypes);
0098
0099   for (size_t j = 0; j < tokensIsoDepElectrons_.size(); ++j) {
0100     iEvent.getByToken(tokensIsoDepElectrons_[j], electronIsoDep[j]);
0101   }
0102
0103   IsoDepositMaps photonIsoDep(nTypes);
0104   for (size_t j = 0; j < tokensIsoDepPhotons_.size(); ++j) {
0105     iEvent.getByToken(tokensIsoDepPhotons_[j], photonIsoDep[j]);
0106   }
0107
0108   IsoDepositVals electronIsoValPFId(nTypes);
0109   IsoDepositVals photonIsoValPFId(nTypes);
0110   // just renaming
0111   const IsoDepositVals* electronIsoVals = &electronIsoValPFId;
0112
0113   for (size_t j = 0; j < tokensIsoValElectronsPFId_.size(); ++j) {
0114     iEvent.getByToken(tokensIsoValElectronsPFId_[j], electronIsoValPFId[j]);
0115   }
0116
0117   for (size_t j = 0; j < tokensIsoValPhotonsPFId_.size(); ++j) {
0118     iEvent.getByToken(tokensIsoValPhotonsPFId_[j], photonIsoValPFId[j]);
0119   }
0120
0121   // Electrons - from reco
0122   if (printElectrons_) {
0123     unsigned nele = gsfElectronH->size();
0124
0125     for (unsigned iele = 0; iele < nele; ++iele) {
0126       reco::GsfElectronRef myElectronRef(gsfElectronH, iele);
0127
0128       double charged = (*(*electronIsoVals)[0])[myElectronRef];
0129       double photon = (*(*electronIsoVals)[1])[myElectronRef];
0130       double neutral = (*(*electronIsoVals)[2])[myElectronRef];
0131
0132       std::cout << "Electron: "
0133                 << " run " << iEvent.id().run() << " lumi " << iEvent.id().luminosityBlock() << " event "
0134                 << iEvent.id().event();
0135       std::cout << " pt " << myElectronRef->pt() << " eta " << myElectronRef->eta() << " phi " << myElectronRef->phi()
0136                 << " charge " << myElectronRef->charge() << " : ";
0137       std::cout << " ChargedIso " << charged;
0138       std::cout << " PhotonIso " << photon;
0139       std::cout << " NeutralHadron Iso " << neutral << std::endl;
0140
0141       if (myElectronRef->isEB()) {
0142         chargedBarrelElectrons_->Fill(charged / myElectronRef->pt());
0143         photonBarrelElectrons_->Fill(photon / myElectronRef->pt());
0144         neutralBarrelElectrons_->Fill(neutral / myElectronRef->pt());
0145         sumBarrelElectrons_->Fill((charged + photon + neutral) / myElectronRef->pt());
0146       } else {
0147         chargedEndcapsElectrons_->Fill(charged / myElectronRef->pt());
0148         photonEndcapsElectrons_->Fill(photon / myElectronRef->pt());
0149         neutralEndcapsElectrons_->Fill(neutral / myElectronRef->pt());
0150         sumEndcapsElectrons_->Fill((charged + photon + neutral) / myElectronRef->pt());
0151       }
0152     }
0153   }
0154
0155   // Photons - from reco
0156   const IsoDepositVals* photonIsoVals = &photonIsoValPFId;
0157
0158   if (printPhotons_) {
0159     unsigned npho = photonH->size();
0160
0161     for (unsigned ipho = 0; ipho < npho; ++ipho) {
0162       reco::PhotonRef myPhotonRef(photonH, ipho);
0163
0164       double charged = (*(*photonIsoVals)[0])[myPhotonRef];
0165       double photon = (*(*photonIsoVals)[1])[myPhotonRef];
0166       double neutral = (*(*photonIsoVals)[2])[myPhotonRef];
0167
0168       std::cout << "Photon: "
0169                 << " run " << iEvent.id().run() << " lumi " << iEvent.id().luminosityBlock() << " event "
0170                 << iEvent.id().event();
0171       std::cout << " pt " << myPhotonRef->pt() << " eta " << myPhotonRef->eta() << " phi " << myPhotonRef->phi()
0172                 << " charge " << myPhotonRef->charge() << " : ";
0173       std::cout << " ChargedIso " << charged;
0174       std::cout << " PhotonIso " << photon;
0175       std::cout << " NeutralHadron Iso " << neutral << std::endl;
0176
0177       if (myPhotonRef->isEB()) {
0178         chargedBarrelPhotons_->Fill(charged / myPhotonRef->pt());
0179         photonBarrelPhotons_->Fill(photon / myPhotonRef->pt());
0180         neutralBarrelPhotons_->Fill(neutral / myPhotonRef->pt());
0181         sumBarrelPhotons_->Fill((charged + photon + neutral) / myPhotonRef->pt());
0182       } else {
0183         chargedEndcapsPhotons_->Fill(charged / myPhotonRef->pt());
0184         photonEndcapsPhotons_->Fill(photon / myPhotonRef->pt());
0185         neutralEndcapsPhotons_->Fill(neutral / myPhotonRef->pt());
0186         sumEndcapsPhotons_->Fill((charged + photon + neutral) / myPhotonRef->pt());
0187       }
0188     }
0189   }
0190 }
0191
0192 DEFINE_FWK_MODULE(PFIsoReaderDemo);