PhotonIdentification/plugins/PhotonIDSimpleAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhotonIDSimpleAnalyzer
0004 // Class:      PhotonIDSimpleAnalyzer
0005 //
0006 /**\class PhotonIDSimpleAnalyzer PhotonIDSimpleAnalyzer.cc RecoEgamma/PhotonIdentification/test/PhotonIDSimpleAnalyzer.cc
0007
0008  Description: Generate various histograms for cuts and important
0009               photon ID parameters using a data sample of photons in QCD events.
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  J. Stilley, A. Askew
0016 //  Editing Author:  M.B. Anderson
0017 //
0018 //         Created:  Fri May 9 11:03:51 CDT 2008
0019 //
0020
0021 ///////////////////////////////////////////////////////////////////////
0022 //                        CMSSW includes                             //
0023 ///////////////////////////////////////////////////////////////////////
0024 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0025 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0026
0027 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0028 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0029 #include "DataFormats/JetReco/interface/CaloJet.h"
0030 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0032 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0033 #include "DataFormats/DetId/interface/DetId.h"
0034 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0035 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0036 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0037 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0038 #include "FWCore/Framework/interface/ESHandle.h"
0039 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0040 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0041 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0042 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0043 ///////////////////////////////////////////////////////////////////////
0044 //                      Root include files                           //
0045 ///////////////////////////////////////////////////////////////////////
0046 #include "TH1F.h"
0047 #include "TH2F.h"
0048 #include "TFile.h"
0049 #include "TTree.h"
0050
0051 // system include files
0052 #include <memory>
0053
0054 // user include files
0055 #include "FWCore/Framework/interface/Frameworkfwd.h"
0056 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0057
0058 #include "FWCore/Framework/interface/Event.h"
0059 #include "FWCore/Framework/interface/MakerMacros.h"
0060
0061 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0062
0063 #include <string>
0064 #include "TH1.h"
0065 #include "TTree.h"
0066
0067 class TFile;
0068
0069 //
0070 // class declaration
0071 //
0072 class PhotonIDSimpleAnalyzer : public edm::one::EDAnalyzer<> {
0073 public:
0074   explicit PhotonIDSimpleAnalyzer(const edm::ParameterSet&);
0075   ~PhotonIDSimpleAnalyzer() override;
0076
0077   void analyze(const edm::Event&, const edm::EventSetup&) override;
0078   void beginJob() override;
0079   void endJob() override;
0080
0081 private:
0082   std::string outputFile_;  // output file
0083   double minPhotonEt_;      // minimum photon Et
0084   double minPhotonAbsEta_;  // min and
0085   double maxPhotonAbsEta_;  // max abs(eta)
0086   double minPhotonR9_;      // minimum R9 = E(3x3)/E(SuperCluster)
0087   double maxPhotonHoverE_;  // maximum HCAL / ECAL
0088   bool createPhotonTTree_;  // Create a TTree of photon variables
0089
0090   // Will be used for creating TTree of photons.
0091   // These names did not have to match those from a phtn->...
0092   // but do match for clarity.
0093   struct struct_recPhoton {
0094     float isolationEcalRecHit;
0095     float isolationHcalRecHit;
0096     float isolationSolidTrkCone;
0097     float isolationHollowTrkCone;
0098     float nTrkSolidCone;
0099     float nTrkHollowCone;
0100     float isEBEtaGap;
0101     float isEBPhiGap;
0102     float isEERingGap;
0103     float isEEDeeGap;
0104     float isEBEEGap;
0105     float r9;
0106     float et;
0107     float eta;
0108     float phi;
0109     float hadronicOverEm;
0110   };
0111   struct_recPhoton recPhoton;
0112
0113   // root file to store histograms
0114   TFile* rootFile_;
0115
0116   // data members for histograms to be filled
0117
0118   // PhotonID Histograms
0119   TH1F* h_isoEcalRecHit_;
0120   TH1F* h_isoHcalRecHit_;
0121   TH1F* h_trk_pt_solid_;
0122   TH1F* h_trk_pt_hollow_;
0123   TH1F* h_ntrk_solid_;
0124   TH1F* h_ntrk_hollow_;
0125   TH1F* h_ebetagap_;
0126   TH1F* h_ebphigap_;
0127   TH1F* h_eeringGap_;
0128   TH1F* h_eedeeGap_;
0129   TH1F* h_ebeeGap_;
0130   TH1F* h_r9_;
0131
0132   // Photon Histograms
0133   TH1F* h_photonEt_;
0134   TH1F* h_photonEta_;
0135   TH1F* h_photonPhi_;
0136   TH1F* h_hadoverem_;
0137
0138   // Photon's SuperCluster Histograms
0139   TH1F* h_photonScEt_;
0140   TH1F* h_photonScEta_;
0141   TH1F* h_photonScPhi_;
0142   TH1F* h_photonScEtaWidth_;
0143
0144   // Composite or Other Histograms
0145   TH1F* h_photonInAnyGap_;
0146   TH1F* h_nPassingPho_;
0147   TH1F* h_nPassEM_;
0148   TH1F* h_nPho_;
0149
0150   // TTree
0151   TTree* tree_PhotonAll_;
0152 };
0153
0154 using namespace std;
0155
0156 ///////////////////////////////////////////////////////////////////////
0157 //                           Constructor                             //
0158 ///////////////////////////////////////////////////////////////////////
0159 PhotonIDSimpleAnalyzer::PhotonIDSimpleAnalyzer(const edm::ParameterSet& ps) {
0160   // Read Parameters from configuration file
0161
0162   // output filename
0163   outputFile_ = ps.getParameter<std::string>("outputFile");
0164   // Read variables that must be passed to allow a
0165   //  supercluster to be placed in histograms as a photon.
0166   minPhotonEt_ = ps.getParameter<double>("minPhotonEt");
0167   minPhotonAbsEta_ = ps.getParameter<double>("minPhotonAbsEta");
0168   maxPhotonAbsEta_ = ps.getParameter<double>("maxPhotonAbsEta");
0169   minPhotonR9_ = ps.getParameter<double>("minPhotonR9");
0170   maxPhotonHoverE_ = ps.getParameter<double>("maxPhotonHoverE");
0171
0172   // Read variable to that decidedes whether
0173   // a TTree of photons is created or not
0174   createPhotonTTree_ = ps.getParameter<bool>("createPhotonTTree");
0175
0176   // open output file to store histograms
0177   rootFile_ = TFile::Open(outputFile_.c_str(), "RECREATE");
0178 }
0179
0180 ///////////////////////////////////////////////////////////////////////
0181 //                            Destructor                             //
0182 ///////////////////////////////////////////////////////////////////////
0183 PhotonIDSimpleAnalyzer::~PhotonIDSimpleAnalyzer() {
0184   // do anything here that needs to be done at desctruction time
0185   // (e.g. close files, deallocate resources etc.)
0186
0187   delete rootFile_;
0188 }
0189
0190 ///////////////////////////////////////////////////////////////////////
0191 //    method called once each job just before starting event loop    //
0192 ///////////////////////////////////////////////////////////////////////
0193 void PhotonIDSimpleAnalyzer::beginJob() {
0194   // go to *OUR* rootfile
0195   rootFile_->cd();
0196
0197   // Book Histograms
0198   // PhotonID Histograms
0199   h_isoEcalRecHit_ = new TH1F("photonEcalIso", "Ecal Rec Hit Isolation", 300, 0, 300);
0200   h_isoHcalRecHit_ = new TH1F("photonHcalIso", "Hcal Rec Hit Isolation", 300, 0, 300);
0201   h_trk_pt_solid_ = new TH1F("photonTrackSolidIso", "Sum of track pT in a cone of #DeltaR", 300, 0, 300);
0202   h_trk_pt_hollow_ = new TH1F("photonTrackHollowIso", "Sum of track pT in a hollow cone", 300, 0, 300);
0203   h_ntrk_solid_ = new TH1F("photonTrackCountSolid", "Number of tracks in a cone of #DeltaR", 100, 0, 100);
0204   h_ntrk_hollow_ = new TH1F("photonTrackCountHollow", "Number of tracks in a hollow cone", 100, 0, 100);
0205   h_ebetagap_ = new TH1F("photonInEBEtagap", "Ecal Barrel eta gap flag", 2, -0.5, 1.5);
0206   h_ebphigap_ = new TH1F("photonInEBEtagap", "Ecal Barrel phi gap flag", 2, -0.5, 1.5);
0207   h_eeringGap_ = new TH1F("photonInEERinggap", "Ecal Endcap ring gap flag", 2, -0.5, 1.5);
0208   h_eedeeGap_ = new TH1F("photonInEEDeegap", "Ecal Endcap dee gap flag", 2, -0.5, 1.5);
0209   h_ebeeGap_ = new TH1F("photonInEEgap", "Ecal Barrel/Endcap gap flag", 2, -0.5, 1.5);
0210   h_r9_ = new TH1F("photonR9", "R9 = E(3x3) / E(SuperCluster)", 300, 0, 3);
0211
0212   // Photon Histograms
0213   h_photonEt_ = new TH1F("photonEt", "Photon E_{T}", 200, 0, 200);
0214   h_photonEta_ = new TH1F("photonEta", "Photon #eta", 800, -4, 4);
0215   h_photonPhi_ = new TH1F("photonPhi", "Photon #phi", 628, -1. * M_PI, M_PI);
0216   h_hadoverem_ = new TH1F("photonHoverE", "Hadronic over EM", 200, 0, 1);
0217
0218   // Photon's SuperCluster Histograms
0219   h_photonScEt_ = new TH1F("photonScEt", "Photon SuperCluster E_{T}", 200, 0, 200);
0220   h_photonScEta_ = new TH1F("photonScEta", "Photon #eta", 800, -4, 4);
0221   h_photonScPhi_ = new TH1F("photonScPhi", "Photon #phi", 628, -1. * M_PI, M_PI);
0222   h_photonScEtaWidth_ = new TH1F("photonScEtaWidth", "#eta-width", 100, 0, .1);
0223
0224   // Composite or Other Histograms
0225   h_photonInAnyGap_ = new TH1F("photonInAnyGap", "Photon in any gap flag", 2, -0.5, 1.5);
0226   h_nPassingPho_ = new TH1F("photonLoosePhoton", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
0227   h_nPassEM_ = new TH1F("photonLooseEM", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
0228   h_nPho_ = new TH1F("photonCount", "Number of photons passing cuts in event", 10, 0, 10);
0229
0230   // Create a TTree of photons if set to 'True' in config file
0231   if (createPhotonTTree_) {
0232     tree_PhotonAll_ = new TTree("TreePhotonAll", "Reconstructed Photon");
0233     tree_PhotonAll_->Branch("recPhoton",
0234                             &recPhoton.isolationEcalRecHit,
0235                             "isolationEcalRecHit/"
0236                             "F:isolationHcalRecHit:isolationSolidTrkCone:isolationHollowTrkCone:nTrkSolidCone:"
0237                             "nTrkHollowCone:isEBGap:isEEGap:isEBEEGap:r9:et:eta:phi:hadronicOverEm");
0238   }
0239 }
0240
0241 ///////////////////////////////////////////////////////////////////////
0242 //                method called to for each event                    //
0243 ///////////////////////////////////////////////////////////////////////
0244 void PhotonIDSimpleAnalyzer::analyze(const edm::Event& evt, const edm::EventSetup& es) {
0245   using namespace std;
0246   using namespace edm;
0247
0248   // grab photons
0249   Handle<reco::PhotonCollection> photonColl;
0250   evt.getByLabel("photons", "", photonColl);
0251
0252   Handle<edm::ValueMap<bool> > loosePhotonQual;
0253   evt.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonQual);
0254   Handle<edm::ValueMap<bool> > looseEMQual;
0255   evt.getByLabel("PhotonIDProd", "PhotonCutBasedIDLooseEM", looseEMQual);
0256   // grab PhotonId objects
0257   //   Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
0258   //   evt.getByLabel("PhotonIDProd", "PhotonAssociatedID", photonIDMapColl);
0259
0260   // create reference to the object types we are interested in
0261   const reco::PhotonCollection* photons = photonColl.product();
0262   const edm::ValueMap<bool>* phoMap = loosePhotonQual.product();
0263   const edm::ValueMap<bool>* lEMMap = looseEMQual.product();
0264   int photonCounter = 0;
0265   int idxpho = 0;
0266   reco::PhotonCollection::const_iterator pho;
0267   for (pho = (*photons).begin(); pho != (*photons).end(); pho++) {
0268     edm::Ref<reco::PhotonCollection> photonref(photonColl, idxpho);
0269     //reco::PhotonIDAssociationCollection::const_iterator photonIter = phoMap->find(photonref);
0270     //const reco::PhotonIDRef &phtn = photonIter->val;
0271     //const reco::PhotonRef &pho = photonIter->key;
0272
0273     float photonEt = pho->et();
0274     float superClusterEt = (pho->superCluster()->energy()) / (cosh(pho->superCluster()->position().eta()));
0275     bool LoosePhotonQu = (*phoMap)[photonref];
0276     h_nPassingPho_->Fill(LoosePhotonQu);
0277     bool LooseEMQu = (*lEMMap)[photonref];
0278     h_nPassEM_->Fill(LooseEMQu);
0279     // Only store photon candidates (SuperClusters) that pass some simple cuts
0280     bool passCuts = (photonEt > minPhotonEt_) && (fabs(pho->eta()) > minPhotonAbsEta_) &&
0281                     (fabs(pho->eta()) < maxPhotonAbsEta_) && (pho->r9() > minPhotonR9_) &&
0282                     (pho->hadronicOverEm() < maxPhotonHoverE_);
0283
0284     if (passCuts) {
0285       ///////////////////////////////////////////////////////
0286       //                fill histograms                    //
0287       ///////////////////////////////////////////////////////
0288       // PhotonID Variables
0289       h_isoEcalRecHit_->Fill(pho->ecalRecHitSumEtConeDR04());
0290       h_isoHcalRecHit_->Fill(pho->hcalTowerSumEtConeDR04());
0291       h_trk_pt_solid_->Fill(pho->trkSumPtSolidConeDR04());
0292       h_trk_pt_hollow_->Fill(pho->trkSumPtHollowConeDR04());
0293       h_ntrk_solid_->Fill(pho->nTrkSolidConeDR04());
0294       h_ntrk_hollow_->Fill(pho->nTrkHollowConeDR04());
0295       h_ebetagap_->Fill(pho->isEBEtaGap());
0296       h_ebphigap_->Fill(pho->isEBPhiGap());
0297       h_eeringGap_->Fill(pho->isEERingGap());
0298       h_eedeeGap_->Fill(pho->isEEDeeGap());
0299       h_ebeeGap_->Fill(pho->isEBEEGap());
0300       h_r9_->Fill(pho->r9());
0301
0302       // Photon Variables
0303       h_photonEt_->Fill(photonEt);
0304       h_photonEta_->Fill(pho->eta());
0305       h_photonPhi_->Fill(pho->phi());
0306       h_hadoverem_->Fill(pho->hadronicOverEm());
0307
0308       // Photon's SuperCluster Variables
0309       // eta is with respect to detector (not physics) vertex,
0310       // thus Et and eta are different from photon.
0311       h_photonScEt_->Fill(superClusterEt);
0312       h_photonScEta_->Fill(pho->superCluster()->position().eta());
0313       h_photonScPhi_->Fill(pho->superCluster()->position().phi());
0314       h_photonScEtaWidth_->Fill(pho->superCluster()->etaWidth());
0315
0316       // It passed photon cuts, mark it
0317
0318       ///////////////////////////////////////////////////////
0319       //                fill TTree (optional)              //
0320       ///////////////////////////////////////////////////////
0321       if (createPhotonTTree_) {
0322         recPhoton.isolationEcalRecHit = pho->ecalRecHitSumEtConeDR04();
0323         recPhoton.isolationHcalRecHit = pho->hcalTowerSumEtConeDR04();
0324         recPhoton.isolationSolidTrkCone = pho->trkSumPtSolidConeDR04();
0325         recPhoton.isolationHollowTrkCone = pho->trkSumPtHollowConeDR04();
0326         recPhoton.nTrkSolidCone = pho->nTrkSolidConeDR04();
0327         recPhoton.nTrkHollowCone = pho->nTrkHollowConeDR04();
0328         recPhoton.isEBEtaGap = pho->isEBEtaGap();
0329         recPhoton.isEBPhiGap = pho->isEBPhiGap();
0330         recPhoton.isEERingGap = pho->isEERingGap();
0331         recPhoton.isEEDeeGap = pho->isEEDeeGap();
0332         recPhoton.isEBEEGap = pho->isEBEEGap();
0333         recPhoton.r9 = pho->r9();
0334         recPhoton.et = pho->et();
0335         recPhoton.eta = pho->eta();
0336         recPhoton.phi = pho->phi();
0337         recPhoton.hadronicOverEm = pho->hadronicOverEm();
0338
0339         // Fill the tree (this records all the recPhoton.* since
0340         // tree_PhotonAll_ was set to point at that.
0341         tree_PhotonAll_->Fill();
0342       }
0343
0344       // Record whether it was near any module gap.
0345       // Very convoluted at the moment.
0346       bool inAnyGap = pho->isEBEEGap() || (pho->isEB() && pho->isEBEtaGap()) || (pho->isEB() && pho->isEBPhiGap()) ||
0347                       (pho->isEE() && pho->isEERingGap()) || (pho->isEE() && pho->isEEDeeGap());
0348       if (inAnyGap) {
0349         h_photonInAnyGap_->Fill(1.0);
0350       } else {
0351         h_photonInAnyGap_->Fill(0.0);
0352       }
0353
0354       photonCounter++;
0355     } else {
0356       // This didn't pass photon cuts, mark it
0357     }
0358     idxpho++;
0359   }  // End Loop over photons
0360   h_nPho_->Fill(photonCounter);
0361 }
0362
0363 ///////////////////////////////////////////////////////////////////////
0364 //    method called once each job just after ending the event loop   //
0365 ///////////////////////////////////////////////////////////////////////
0366 void PhotonIDSimpleAnalyzer::endJob() {
0367   // go to *OUR* root file and store histograms
0368   rootFile_->cd();
0369
0370   // PhotonID Histograms
0371   h_isoEcalRecHit_->Write();
0372   h_isoHcalRecHit_->Write();
0373   h_trk_pt_solid_->Write();
0374   h_trk_pt_hollow_->Write();
0375   h_ntrk_solid_->Write();
0376   h_ntrk_hollow_->Write();
0377   h_ebetagap_->Write();
0378   h_ebphigap_->Write();
0379   h_eeringGap_->Write();
0380   h_eedeeGap_->Write();
0381   h_ebeeGap_->Write();
0382   h_r9_->Write();
0383
0384   // Photon Histograms
0385   h_photonEt_->Write();
0386   h_photonEta_->Write();
0387   h_photonPhi_->Write();
0388   h_hadoverem_->Write();
0389
0390   // Photon's SuperCluster Histograms
0391   h_photonScEt_->Write();
0392   h_photonScEta_->Write();
0393   h_photonScPhi_->Write();
0394   h_photonScEtaWidth_->Write();
0395
0396   // Composite or Other Histograms
0397   h_photonInAnyGap_->Write();
0398   h_nPassingPho_->Write();
0399   h_nPassEM_->Write();
0400   h_nPho_->Write();
0401
0402   // Write the root file (really writes the TTree)
0403   rootFile_->Write();
0404   rootFile_->Close();
0405 }
0406
0407 //define this as a plug-in
0408 DEFINE_FWK_MODULE(PhotonIDSimpleAnalyzer);