Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    PhotonIDSimpleAnalyzer
00004 // Class:      PhotonIDSimpleAnalyzer
00005 // 
00014 //
00015 // Original Author:  J. Stilley, A. Askew
00016 //  Editing Author:  M.B. Anderson
00017 //
00018 //         Created:  Fri May 9 11:03:51 CDT 2008
00019 // $Id: PhotonIDSimpleAnalyzer.cc,v 1.10 2010/01/14 17:29:32 nancy Exp $
00020 //
00022 //                    header file for this analyzer                  //
00024 #include "RecoEgamma/PhotonIdentification/plugins/PhotonIDSimpleAnalyzer.h"
00025 
00027 //                        CMSSW includes                             //
00029 //#include "DataFormats/EgammaCandidates/interface/PhotonIDFwd.h"
00030 //#include "DataFormats/EgammaCandidates/interface/PhotonID.h"
00031 //#include "DataFormats/EgammaCandidates/interface/PhotonIDAssociation.h"
00032 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00033 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00034 
00035 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00036 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00037 #include "DataFormats/JetReco/interface/CaloJet.h"
00038 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00039 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00040 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00041 #include "DataFormats/DetId/interface/DetId.h"
00042 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00043 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00044 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00045 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
00046 #include "FWCore/Framework/interface/ESHandle.h"
00047 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00048 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00049 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00050 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00052 //                      Root include files                           //
00054 #include "TH1F.h"
00055 #include "TH2F.h"
00056 #include "TFile.h"
00057 #include "TMath.h"
00058 #include "TTree.h"
00059 
00060 using namespace std;
00061 
00063 //                           Constructor                             //
00065 PhotonIDSimpleAnalyzer::PhotonIDSimpleAnalyzer(const edm::ParameterSet& ps)
00066 {
00067   // Read Parameters from configuration file
00068 
00069   // output filename
00070   outputFile_   = ps.getParameter<std::string>("outputFile");
00071   // Read variables that must be passed to allow a 
00072   //  supercluster to be placed in histograms as a photon.
00073   minPhotonEt_     = ps.getParameter<double>("minPhotonEt");
00074   minPhotonAbsEta_ = ps.getParameter<double>("minPhotonAbsEta");
00075   maxPhotonAbsEta_ = ps.getParameter<double>("maxPhotonAbsEta");
00076   minPhotonR9_     = ps.getParameter<double>("minPhotonR9");
00077   maxPhotonHoverE_ = ps.getParameter<double>("maxPhotonHoverE");
00078 
00079   // Read variable to that decidedes whether
00080   // a TTree of photons is created or not
00081   createPhotonTTree_ = ps.getParameter<bool>("createPhotonTTree");
00082 
00083   // open output file to store histograms
00084   rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE");
00085 }
00086 
00088 //                            Destructor                             //
00090 PhotonIDSimpleAnalyzer::~PhotonIDSimpleAnalyzer()
00091 {
00092 
00093 // do anything here that needs to be done at desctruction time
00094 // (e.g. close files, deallocate resources etc.)
00095 
00096   delete rootFile_;
00097 
00098 }
00099 
00101 //    method called once each job just before starting event loop    //
00103 void 
00104 PhotonIDSimpleAnalyzer::beginJob()
00105 {
00106 
00107   // go to *OUR* rootfile
00108   rootFile_->cd();
00109 
00110   // Book Histograms
00111   // PhotonID Histograms
00112   h_isoEcalRecHit_ = new TH1F("photonEcalIso",          "Ecal Rec Hit Isolation", 300, 0, 300);
00113   h_isoHcalRecHit_ = new TH1F("photonHcalIso",          "Hcal Rec Hit Isolation", 300, 0, 300);
00114   h_trk_pt_solid_  = new TH1F("photonTrackSolidIso",    "Sum of track pT in a cone of #DeltaR" , 300, 0, 300);
00115   h_trk_pt_hollow_ = new TH1F("photonTrackHollowIso",   "Sum of track pT in a hollow cone" ,     300, 0, 300);
00116   h_ntrk_solid_    = new TH1F("photonTrackCountSolid",  "Number of tracks in a cone of #DeltaR", 100, 0, 100);
00117   h_ntrk_hollow_   = new TH1F("photonTrackCountHollow", "Number of tracks in a hollow cone",     100, 0, 100);
00118   h_ebetagap_         = new TH1F("photonInEBEtagap",          "Ecal Barrel eta gap flag",  2, -0.5, 1.5);
00119   h_ebphigap_         = new TH1F("photonInEBEtagap",          "Ecal Barrel phi gap flag",  2, -0.5, 1.5);
00120   h_eeringGap_         = new TH1F("photonInEERinggap",          "Ecal Endcap ring gap flag",  2, -0.5, 1.5);
00121   h_eedeeGap_         = new TH1F("photonInEEDeegap",          "Ecal Endcap dee gap flag",  2, -0.5, 1.5);
00122   h_ebeeGap_       = new TH1F("photonInEEgap",          "Ecal Barrel/Endcap gap flag",  2, -0.5, 1.5);
00123   h_r9_            = new TH1F("photonR9",               "R9 = E(3x3) / E(SuperCluster)", 300, 0, 3);
00124 
00125   // Photon Histograms
00126   h_photonEt_      = new TH1F("photonEt",     "Photon E_{T}",  200,  0, 200);
00127   h_photonEta_     = new TH1F("photonEta",    "Photon #eta",   800, -4,   4);
00128   h_photonPhi_     = new TH1F("photonPhi",    "Photon #phi",   628, -1.*TMath::Pi(), TMath::Pi());
00129   h_hadoverem_     = new TH1F("photonHoverE", "Hadronic over EM", 200, 0, 1);
00130 
00131   // Photon's SuperCluster Histograms
00132   h_photonScEt_       = new TH1F("photonScEt",  "Photon SuperCluster E_{T}", 200,  0, 200);
00133   h_photonScEta_      = new TH1F("photonScEta", "Photon #eta",               800, -4,   4);
00134   h_photonScPhi_      = new TH1F("photonScPhi", "Photon #phi",628, -1.*TMath::Pi(), TMath::Pi());
00135   h_photonScEtaWidth_ = new TH1F("photonScEtaWidth","#eta-width",            100,  0,  .1);
00136 
00137   // Composite or Other Histograms
00138   h_photonInAnyGap_   = new TH1F("photonInAnyGap",     "Photon in any gap flag",  2, -0.5, 1.5);
00139   h_nPassingPho_      = new TH1F("photonLoosePhoton", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
00140   h_nPassEM_          = new TH1F("photonLooseEM", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
00141   h_nPho_             = new TH1F("photonCount",        "Number of photons passing cuts in event",  10,  0,  10);
00142 
00143   // Create a TTree of photons if set to 'True' in config file
00144   if ( createPhotonTTree_ ) {
00145     tree_PhotonAll_     = new TTree("TreePhotonAll", "Reconstructed Photon");
00146     tree_PhotonAll_->Branch("recPhoton", &recPhoton.isolationEcalRecHit, "isolationEcalRecHit/F:isolationHcalRecHit:isolationSolidTrkCone:isolationHollowTrkCone:nTrkSolidCone:nTrkHollowCone:isEBGap:isEEGap:isEBEEGap:r9:et:eta:phi:hadronicOverEm");
00147   }
00148 }
00149 
00151 //                method called to for each event                    //
00153 void
00154 PhotonIDSimpleAnalyzer::analyze(const edm::Event& evt, const edm::EventSetup& es)
00155 {
00156   
00157   using namespace std;
00158   using namespace edm;
00159   
00160   // grab photons
00161   Handle<reco::PhotonCollection> photonColl;
00162   evt.getByLabel("photons", "", photonColl);
00163 
00164   Handle<edm::ValueMap<Bool_t> > loosePhotonQual;
00165   evt.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonQual);
00166   Handle<edm::ValueMap<Bool_t> > looseEMQual;
00167   evt.getByLabel("PhotonIDProd","PhotonCutBasedIDLooseEM",looseEMQual);
00168   // grab PhotonId objects  
00169 //   Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
00170 //   evt.getByLabel("PhotonIDProd", "PhotonAssociatedID", photonIDMapColl);
00171   
00172   // create reference to the object types we are interested in
00173   const reco::PhotonCollection *photons = photonColl.product();  
00174   const edm::ValueMap<Bool_t> *phoMap = loosePhotonQual.product();
00175   const edm::ValueMap<Bool_t> *lEMMap = looseEMQual.product();
00176   int photonCounter = 0;
00177   int idxpho=0;
00178   reco::PhotonCollection::const_iterator pho;
00179   for (pho = (*photons).begin(); pho!= (*photons).end(); pho++){   
00180     
00181     edm::Ref<reco::PhotonCollection> photonref(photonColl, idxpho);
00182     //reco::PhotonIDAssociationCollection::const_iterator photonIter = phoMap->find(photonref);
00183     //const reco::PhotonIDRef &phtn = photonIter->val;
00184     //const reco::PhotonRef &pho = photonIter->key;
00185 
00186     float photonEt       = pho->et();
00187     float superClusterEt = (pho->superCluster()->energy())/(cosh(pho->superCluster()->position().eta()));
00188     Bool_t LoosePhotonQu = (*phoMap)[photonref];
00189     h_nPassingPho_->Fill(LoosePhotonQu);
00190     Bool_t LooseEMQu = (*lEMMap)[photonref];
00191     h_nPassEM_->Fill(LooseEMQu);
00192     // Only store photon candidates (SuperClusters) that pass some simple cuts
00193     bool passCuts = (              photonEt > minPhotonEt_     ) &&
00194                     (      fabs(pho->eta()) > minPhotonAbsEta_ ) &&
00195                     (      fabs(pho->eta()) < maxPhotonAbsEta_ ) &&
00196                     (             pho->r9() > minPhotonR9_     ) &&
00197                     ( pho->hadronicOverEm() < maxPhotonHoverE_ ) ;
00198 
00199     if ( passCuts )
00200     {
00202       //                fill histograms                    //
00204       // PhotonID Variables
00205       h_isoEcalRecHit_->Fill(pho->ecalRecHitSumEtConeDR04());
00206       h_isoHcalRecHit_->Fill(pho->hcalTowerSumEtConeDR04());
00207       h_trk_pt_solid_ ->Fill(pho->trkSumPtSolidConeDR04());
00208       h_trk_pt_hollow_->Fill(pho->trkSumPtHollowConeDR04());
00209       h_ntrk_solid_->   Fill(pho->nTrkSolidConeDR04());
00210       h_ntrk_hollow_->  Fill(pho->nTrkHollowConeDR04());
00211       h_ebetagap_->        Fill(pho->isEBEtaGap());
00212       h_ebphigap_->        Fill(pho->isEBPhiGap());
00213       h_eeringGap_->        Fill(pho->isEERingGap()); 
00214       h_eedeeGap_->        Fill(pho->isEEDeeGap()); 
00215       h_ebeeGap_->      Fill(pho->isEBEEGap());
00216       h_r9_->           Fill(pho->r9());
00217 
00218       // Photon Variables
00219       h_photonEt_->  Fill(photonEt);
00220       h_photonEta_-> Fill(pho->eta());
00221       h_photonPhi_-> Fill(pho->phi());
00222       h_hadoverem_-> Fill(pho->hadronicOverEm());
00223 
00224       // Photon's SuperCluster Variables
00225       // eta is with respect to detector (not physics) vertex,
00226       // thus Et and eta are different from photon.
00227       h_photonScEt_->      Fill(superClusterEt);
00228       h_photonScEta_->     Fill(pho->superCluster()->position().eta());
00229       h_photonScPhi_->     Fill(pho->superCluster()->position().phi());
00230       h_photonScEtaWidth_->Fill(pho->superCluster()->etaWidth());
00231       
00232       // It passed photon cuts, mark it
00233 
00235       //                fill TTree (optional)              //
00237       if ( createPhotonTTree_ ) {
00238         recPhoton.isolationEcalRecHit    = pho->ecalRecHitSumEtConeDR04();
00239         recPhoton.isolationHcalRecHit    = pho->hcalTowerSumEtConeDR04();
00240         recPhoton.isolationSolidTrkCone  = pho->trkSumPtSolidConeDR04();
00241         recPhoton.isolationHollowTrkCone = pho->trkSumPtHollowConeDR04();
00242         recPhoton.nTrkSolidCone          = pho->nTrkSolidConeDR04();
00243         recPhoton.nTrkHollowCone         = pho->nTrkHollowConeDR04();
00244         recPhoton.isEBEtaGap                = pho->isEBEtaGap();
00245         recPhoton.isEBPhiGap                = pho->isEBPhiGap();
00246         recPhoton.isEERingGap                = pho->isEERingGap();
00247         recPhoton.isEEDeeGap                = pho->isEEDeeGap();
00248         recPhoton.isEBEEGap              = pho->isEBEEGap();
00249         recPhoton.r9                     = pho->r9();
00250         recPhoton.et                     = pho->et();
00251         recPhoton.eta                    = pho->eta();
00252         recPhoton.phi                    = pho->phi();
00253         recPhoton.hadronicOverEm         = pho->hadronicOverEm();
00254 
00255         // Fill the tree (this records all the recPhoton.* since
00256         // tree_PhotonAll_ was set to point at that.
00257         tree_PhotonAll_->Fill();
00258       }
00259 
00260       // Record whether it was near any module gap.
00261       // Very convoluted at the moment.
00262       bool inAnyGap = pho->isEBEEGap() || (pho->isEB()&&pho->isEBEtaGap()) ||(pho->isEB()&&pho->isEBPhiGap())  || (pho->isEE()&&pho->isEERingGap()) || (pho->isEE()&&pho->isEEDeeGap()) ;
00263       if (inAnyGap) {
00264         h_photonInAnyGap_->Fill(1.0);
00265       } else {
00266         h_photonInAnyGap_->Fill(0.0);
00267       }
00268 
00269       photonCounter++;
00270     } 
00271     else
00272     {
00273       // This didn't pass photon cuts, mark it
00274    
00275     }
00276     idxpho++;
00277   } // End Loop over photons
00278   h_nPho_->Fill(photonCounter);
00279  
00280 }
00281 
00283 //    method called once each job just after ending the event loop   //
00285 void 
00286 PhotonIDSimpleAnalyzer::endJob()
00287 {
00288 
00289   // go to *OUR* root file and store histograms
00290   rootFile_->cd();
00291 
00292   // PhotonID Histograms
00293   h_isoEcalRecHit_->Write();
00294   h_isoHcalRecHit_->Write();
00295   h_trk_pt_solid_-> Write();
00296   h_trk_pt_hollow_->Write();
00297   h_ntrk_solid_->   Write();
00298   h_ntrk_hollow_->  Write();
00299   h_ebetagap_->     Write();
00300   h_ebphigap_->     Write();
00301   h_eeringGap_->     Write();
00302   h_eedeeGap_->     Write();
00303   h_ebeeGap_->   Write();
00304   h_r9_->        Write();
00305 
00306   // Photon Histograms
00307   h_photonEt_->  Write();
00308   h_photonEta_-> Write();
00309   h_photonPhi_-> Write();
00310   h_hadoverem_-> Write();
00311 
00312   // Photon's SuperCluster Histograms
00313   h_photonScEt_->      Write();
00314   h_photonScEta_->     Write();
00315   h_photonScPhi_->     Write();
00316   h_photonScEtaWidth_->Write();
00317 
00318   // Composite or Other Histograms
00319   h_photonInAnyGap_->Write();
00320   h_nPassingPho_->   Write();
00321   h_nPassEM_->       Write();
00322   h_nPho_->          Write();
00323 
00324   // Write the root file (really writes the TTree)
00325   rootFile_->Write();
00326   rootFile_->Close();
00327 
00328 }
00329 
00330 //define this as a plug-in
00331 // DEFINE_FWK_MODULE(PhotonIDSimpleAnalyzer);
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/RecoEgamma/PhotonIdentification/plugins/PhotonIDSimpleAnalyzer.cc
