CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PhotonIDSimpleAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: PhotonIDSimpleAnalyzer
4 // Class: PhotonIDSimpleAnalyzer
5 //
14 //
15 // Original Author: J. Stilley, A. Askew
16 // Editing Author: M.B. Anderson
17 //
18 // Created: Fri May 9 11:03:51 CDT 2008
19 // $Id: PhotonIDSimpleAnalyzer.cc,v 1.10 2010/01/14 17:29:32 nancy Exp $
20 //
22 // header file for this analyzer //
24 #include "RecoEgamma/PhotonIdentification/plugins/PhotonIDSimpleAnalyzer.h"
25 
27 // CMSSW includes //
29 //#include "DataFormats/EgammaCandidates/interface/PhotonIDFwd.h"
30 //#include "DataFormats/EgammaCandidates/interface/PhotonID.h"
31 //#include "DataFormats/EgammaCandidates/interface/PhotonIDAssociation.h"
32 #include "DataFormats/EgammaCandidates/interface/Photon.h"
33 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
34 
35 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
36 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
37 #include "DataFormats/JetReco/interface/CaloJet.h"
38 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
39 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
40 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
41 #include "DataFormats/DetId/interface/DetId.h"
42 #include "DataFormats/EcalDetId/interface/EBDetId.h"
43 #include "Geometry/Records/interface/IdealGeometryRecord.h"
44 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
45 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
46 #include "FWCore/Framework/interface/ESHandle.h"
47 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
48 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
49 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
50 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
52 // Root include files //
54 #include "TH1F.h"
55 #include "TH2F.h"
56 #include "TFile.h"
57 #include "TMath.h"
58 #include "TTree.h"
59 
60 using namespace std;
61 
63 // Constructor //
65 PhotonIDSimpleAnalyzer::PhotonIDSimpleAnalyzer(const edm::ParameterSet& ps)
66 {
67  // Read Parameters from configuration file
68 
69  // output filename
70  outputFile_ = ps.getParameter<std::string>("outputFile");
71  // Read variables that must be passed to allow a
72  // supercluster to be placed in histograms as a photon.
73  minPhotonEt_ = ps.getParameter<double>("minPhotonEt");
74  minPhotonAbsEta_ = ps.getParameter<double>("minPhotonAbsEta");
75  maxPhotonAbsEta_ = ps.getParameter<double>("maxPhotonAbsEta");
76  minPhotonR9_ = ps.getParameter<double>("minPhotonR9");
77  maxPhotonHoverE_ = ps.getParameter<double>("maxPhotonHoverE");
78 
79  // Read variable to that decidedes whether
80  // a TTree of photons is created or not
81  createPhotonTTree_ = ps.getParameter<bool>("createPhotonTTree");
82 
83  // open output file to store histograms
84  rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE");
85 }
86 
88 // Destructor //
90 PhotonIDSimpleAnalyzer::~PhotonIDSimpleAnalyzer()
91 {
92 
93 // do anything here that needs to be done at desctruction time
94 // (e.g. close files, deallocate resources etc.)
95 
96  delete rootFile_;
97 
98 }
99 
101 // method called once each job just before starting event loop //
103 void
104 PhotonIDSimpleAnalyzer::beginJob()
105 {
106 
107  // go to *OUR* rootfile
108  rootFile_->cd();
109 
110  // Book Histograms
111  // PhotonID Histograms
112  h_isoEcalRecHit_ = new TH1F("photonEcalIso", "Ecal Rec Hit Isolation", 300, 0, 300);
113  h_isoHcalRecHit_ = new TH1F("photonHcalIso", "Hcal Rec Hit Isolation", 300, 0, 300);
114  h_trk_pt_solid_ = new TH1F("photonTrackSolidIso", "Sum of track pT in a cone of #DeltaR" , 300, 0, 300);
115  h_trk_pt_hollow_ = new TH1F("photonTrackHollowIso", "Sum of track pT in a hollow cone" , 300, 0, 300);
116  h_ntrk_solid_ = new TH1F("photonTrackCountSolid", "Number of tracks in a cone of #DeltaR", 100, 0, 100);
117  h_ntrk_hollow_ = new TH1F("photonTrackCountHollow", "Number of tracks in a hollow cone", 100, 0, 100);
118  h_ebetagap_ = new TH1F("photonInEBEtagap", "Ecal Barrel eta gap flag", 2, -0.5, 1.5);
119  h_ebphigap_ = new TH1F("photonInEBEtagap", "Ecal Barrel phi gap flag", 2, -0.5, 1.5);
120  h_eeringGap_ = new TH1F("photonInEERinggap", "Ecal Endcap ring gap flag", 2, -0.5, 1.5);
121  h_eedeeGap_ = new TH1F("photonInEEDeegap", "Ecal Endcap dee gap flag", 2, -0.5, 1.5);
122  h_ebeeGap_ = new TH1F("photonInEEgap", "Ecal Barrel/Endcap gap flag", 2, -0.5, 1.5);
123  h_r9_ = new TH1F("photonR9", "R9 = E(3x3) / E(SuperCluster)", 300, 0, 3);
124 
125  // Photon Histograms
126  h_photonEt_ = new TH1F("photonEt", "Photon E_{T}", 200, 0, 200);
127  h_photonEta_ = new TH1F("photonEta", "Photon #eta", 800, -4, 4);
128  h_photonPhi_ = new TH1F("photonPhi", "Photon #phi", 628, -1.*TMath::Pi(), TMath::Pi());
129  h_hadoverem_ = new TH1F("photonHoverE", "Hadronic over EM", 200, 0, 1);
130 
131  // Photon's SuperCluster Histograms
132  h_photonScEt_ = new TH1F("photonScEt", "Photon SuperCluster E_{T}", 200, 0, 200);
133  h_photonScEta_ = new TH1F("photonScEta", "Photon #eta", 800, -4, 4);
134  h_photonScPhi_ = new TH1F("photonScPhi", "Photon #phi",628, -1.*TMath::Pi(), TMath::Pi());
135  h_photonScEtaWidth_ = new TH1F("photonScEtaWidth","#eta-width", 100, 0, .1);
136 
137  // Composite or Other Histograms
138  h_photonInAnyGap_ = new TH1F("photonInAnyGap", "Photon in any gap flag", 2, -0.5, 1.5);
139  h_nPassingPho_ = new TH1F("photonLoosePhoton", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
140  h_nPassEM_ = new TH1F("photonLooseEM", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
141  h_nPho_ = new TH1F("photonCount", "Number of photons passing cuts in event", 10, 0, 10);
142 
143  // Create a TTree of photons if set to 'True' in config file
144  if ( createPhotonTTree_ ) {
145  tree_PhotonAll_ = new TTree("TreePhotonAll", "Reconstructed Photon");
146  tree_PhotonAll_->Branch("recPhoton", &recPhoton.isolationEcalRecHit, "isolationEcalRecHit/F:isolationHcalRecHit:isolationSolidTrkCone:isolationHollowTrkCone:nTrkSolidCone:nTrkHollowCone:isEBGap:isEEGap:isEBEEGap:r9:et:eta:phi:hadronicOverEm");
147  }
148 }
149 
151 // method called to for each event //
153 void
154 PhotonIDSimpleAnalyzer::analyze(const edm::Event& evt, const edm::EventSetup& es)
155 {
156 
157  using namespace std;
158  using namespace edm;
159 
160  // grab photons
161  Handle<reco::PhotonCollection> photonColl;
162  evt.getByLabel("photons", "", photonColl);
163 
164  Handle<edm::ValueMap<Bool_t> > loosePhotonQual;
165  evt.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonQual);
166  Handle<edm::ValueMap<Bool_t> > looseEMQual;
167  evt.getByLabel("PhotonIDProd","PhotonCutBasedIDLooseEM",looseEMQual);
168  // grab PhotonId objects
169 // Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
170 // evt.getByLabel("PhotonIDProd", "PhotonAssociatedID", photonIDMapColl);
171 
172  // create reference to the object types we are interested in
173  const reco::PhotonCollection *photons = photonColl.product();
174  const edm::ValueMap<Bool_t> *phoMap = loosePhotonQual.product();
175  const edm::ValueMap<Bool_t> *lEMMap = looseEMQual.product();
176  int photonCounter = 0;
177  int idxpho=0;
178  reco::PhotonCollection::const_iterator pho;
179  for (pho = (*photons).begin(); pho!= (*photons).end(); pho++){
180 
181  edm::Ref<reco::PhotonCollection> photonref(photonColl, idxpho);
182  //reco::PhotonIDAssociationCollection::const_iterator photonIter = phoMap->find(photonref);
183  //const reco::PhotonIDRef &phtn = photonIter->val;
184  //const reco::PhotonRef &pho = photonIter->key;
185 
186  float photonEt = pho->et();
187  float superClusterEt = (pho->superCluster()->energy())/(cosh(pho->superCluster()->position().eta()));
188  Bool_t LoosePhotonQu = (*phoMap)[photonref];
189  h_nPassingPho_->Fill(LoosePhotonQu);
190  Bool_t LooseEMQu = (*lEMMap)[photonref];
191  h_nPassEM_->Fill(LooseEMQu);
192  // Only store photon candidates (SuperClusters) that pass some simple cuts
193  bool passCuts = ( photonEt > minPhotonEt_ ) &&
194  ( fabs(pho->eta()) > minPhotonAbsEta_ ) &&
195  ( fabs(pho->eta()) < maxPhotonAbsEta_ ) &&
196  ( pho->r9() > minPhotonR9_ ) &&
197  ( pho->hadronicOverEm() < maxPhotonHoverE_ ) ;
198 
199  if ( passCuts )
200  {
202  // fill histograms //
204  // PhotonID Variables
205  h_isoEcalRecHit_->Fill(pho->ecalRecHitSumEtConeDR04());
206  h_isoHcalRecHit_->Fill(pho->hcalTowerSumEtConeDR04());
207  h_trk_pt_solid_ ->Fill(pho->trkSumPtSolidConeDR04());
208  h_trk_pt_hollow_->Fill(pho->trkSumPtHollowConeDR04());
209  h_ntrk_solid_-> Fill(pho->nTrkSolidConeDR04());
210  h_ntrk_hollow_-> Fill(pho->nTrkHollowConeDR04());
211  h_ebetagap_-> Fill(pho->isEBEtaGap());
212  h_ebphigap_-> Fill(pho->isEBPhiGap());
213  h_eeringGap_-> Fill(pho->isEERingGap());
214  h_eedeeGap_-> Fill(pho->isEEDeeGap());
215  h_ebeeGap_-> Fill(pho->isEBEEGap());
216  h_r9_-> Fill(pho->r9());
217 
218  // Photon Variables
219  h_photonEt_-> Fill(photonEt);
220  h_photonEta_-> Fill(pho->eta());
221  h_photonPhi_-> Fill(pho->phi());
222  h_hadoverem_-> Fill(pho->hadronicOverEm());
223 
224  // Photon's SuperCluster Variables
225  // eta is with respect to detector (not physics) vertex,
226  // thus Et and eta are different from photon.
227  h_photonScEt_-> Fill(superClusterEt);
228  h_photonScEta_-> Fill(pho->superCluster()->position().eta());
229  h_photonScPhi_-> Fill(pho->superCluster()->position().phi());
230  h_photonScEtaWidth_->Fill(pho->superCluster()->etaWidth());
231 
232  // It passed photon cuts, mark it
233 
235  // fill TTree (optional) //
237  if ( createPhotonTTree_ ) {
238  recPhoton.isolationEcalRecHit = pho->ecalRecHitSumEtConeDR04();
239  recPhoton.isolationHcalRecHit = pho->hcalTowerSumEtConeDR04();
240  recPhoton.isolationSolidTrkCone = pho->trkSumPtSolidConeDR04();
241  recPhoton.isolationHollowTrkCone = pho->trkSumPtHollowConeDR04();
242  recPhoton.nTrkSolidCone = pho->nTrkSolidConeDR04();
243  recPhoton.nTrkHollowCone = pho->nTrkHollowConeDR04();
244  recPhoton.isEBEtaGap = pho->isEBEtaGap();
245  recPhoton.isEBPhiGap = pho->isEBPhiGap();
246  recPhoton.isEERingGap = pho->isEERingGap();
247  recPhoton.isEEDeeGap = pho->isEEDeeGap();
248  recPhoton.isEBEEGap = pho->isEBEEGap();
249  recPhoton.r9 = pho->r9();
250  recPhoton.et = pho->et();
251  recPhoton.eta = pho->eta();
252  recPhoton.phi = pho->phi();
253  recPhoton.hadronicOverEm = pho->hadronicOverEm();
254 
255  // Fill the tree (this records all the recPhoton.* since
256  // tree_PhotonAll_ was set to point at that.
257  tree_PhotonAll_->Fill();
258  }
259 
260  // Record whether it was near any module gap.
261  // Very convoluted at the moment.
262  bool inAnyGap = pho->isEBEEGap() || (pho->isEB()&&pho->isEBEtaGap()) ||(pho->isEB()&&pho->isEBPhiGap()) || (pho->isEE()&&pho->isEERingGap()) || (pho->isEE()&&pho->isEEDeeGap()) ;
263  if (inAnyGap) {
264  h_photonInAnyGap_->Fill(1.0);
265  } else {
266  h_photonInAnyGap_->Fill(0.0);
267  }
268 
269  photonCounter++;
270  }
271  else
272  {
273  // This didn't pass photon cuts, mark it
274 
275  }
276  idxpho++;
277  } // End Loop over photons
278  h_nPho_->Fill(photonCounter);
279 
280 }
281 
283 // method called once each job just after ending the event loop //
285 void
286 PhotonIDSimpleAnalyzer::endJob()
287 {
288 
289  // go to *OUR* root file and store histograms
290  rootFile_->cd();
291 
292  // PhotonID Histograms
293  h_isoEcalRecHit_->Write();
294  h_isoHcalRecHit_->Write();
295  h_trk_pt_solid_-> Write();
296  h_trk_pt_hollow_->Write();
297  h_ntrk_solid_-> Write();
298  h_ntrk_hollow_-> Write();
299  h_ebetagap_-> Write();
300  h_ebphigap_-> Write();
301  h_eeringGap_-> Write();
302  h_eedeeGap_-> Write();
303  h_ebeeGap_-> Write();
304  h_r9_-> Write();
305 
306  // Photon Histograms
307  h_photonEt_-> Write();
308  h_photonEta_-> Write();
309  h_photonPhi_-> Write();
310  h_hadoverem_-> Write();
311 
312  // Photon's SuperCluster Histograms
313  h_photonScEt_-> Write();
314  h_photonScEta_-> Write();
315  h_photonScPhi_-> Write();
316  h_photonScEtaWidth_->Write();
317 
318  // Composite or Other Histograms
319  h_photonInAnyGap_->Write();
320  h_nPassingPho_-> Write();
321  h_nPassEM_-> Write();
322  h_nPho_-> Write();
323 
324  // Write the root file (really writes the TTree)
325  rootFile_->Write();
326  rootFile_->Close();
327 
328 }
329 
330 //define this as a plug-in
331 // DEFINE_FWK_MODULE(PhotonIDSimpleAnalyzer);
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
PhotonIDSimpleAnalyzer::PhotonIDSimpleAnalyzer
PhotonIDSimpleAnalyzer(const edm::ParameterSet &)
Definition: PhotonIDSimpleAnalyzer.cc:65
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
PhotonIDSimpleAnalyzer::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: PhotonIDSimpleAnalyzer.cc:154
reco::PhotonCollection
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
Generated for CMSSW Reference Manual by doxygen 1.8.5