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