CMS 3D CMS Logo

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 //
20 
22 // CMSSW includes //
26 
44 // Root include files //
46 #include "TH1F.h"
47 #include "TH2F.h"
48 #include "TFile.h"
49 #include "TTree.h"
50 
51 // system include files
52 #include <memory>
53 
54 // user include files
57 
60 
62 
63 #include <string>
64 #include "TH1.h"
65 #include "TTree.h"
66 
67 class TFile;
68 
69 //
70 // class declaration
71 //
73 public:
75  ~PhotonIDSimpleAnalyzer() override;
76 
77  void analyze(const edm::Event&, const edm::EventSetup&) override;
78  void beginJob() override;
79  void endJob() override;
80 
81 private:
82  std::string outputFile_; // output file
83  double minPhotonEt_; // minimum photon Et
84  double minPhotonAbsEta_; // min and
85  double maxPhotonAbsEta_; // max abs(eta)
86  double minPhotonR9_; // minimum R9 = E(3x3)/E(SuperCluster)
87  double maxPhotonHoverE_; // maximum HCAL / ECAL
88  bool createPhotonTTree_; // Create a TTree of photon variables
89 
90  // Will be used for creating TTree of photons.
91  // These names did not have to match those from a phtn->...
92  // but do match for clarity.
100  float isEBEtaGap;
101  float isEBPhiGap;
102  float isEERingGap;
103  float isEEDeeGap;
104  float isEBEEGap;
105  float r9;
106  float et;
107  float eta;
108  float phi;
110  };
112 
113  // root file to store histograms
114  TFile* rootFile_;
115 
116  // data members for histograms to be filled
117 
118  // PhotonID Histograms
125  TH1F* h_ebetagap_;
126  TH1F* h_ebphigap_;
128  TH1F* h_eedeeGap_;
129  TH1F* h_ebeeGap_;
130  TH1F* h_r9_;
131 
132  // Photon Histograms
133  TH1F* h_photonEt_;
137 
138  // Photon's SuperCluster Histograms
143 
144  // Composite or Other Histograms
147  TH1F* h_nPassEM_;
148  TH1F* h_nPho_;
149 
150  // TTree
152 };
153 
154 using namespace std;
155 
157 // Constructor //
160  // Read Parameters from configuration file
161 
162  // output filename
163  outputFile_ = ps.getParameter<std::string>("outputFile");
164  // Read variables that must be passed to allow a
165  // supercluster to be placed in histograms as a photon.
166  minPhotonEt_ = ps.getParameter<double>("minPhotonEt");
167  minPhotonAbsEta_ = ps.getParameter<double>("minPhotonAbsEta");
168  maxPhotonAbsEta_ = ps.getParameter<double>("maxPhotonAbsEta");
169  minPhotonR9_ = ps.getParameter<double>("minPhotonR9");
170  maxPhotonHoverE_ = ps.getParameter<double>("maxPhotonHoverE");
171 
172  // Read variable to that decidedes whether
173  // a TTree of photons is created or not
174  createPhotonTTree_ = ps.getParameter<bool>("createPhotonTTree");
175 
176  // open output file to store histograms
177  rootFile_ = TFile::Open(outputFile_.c_str(), "RECREATE");
178 }
179 
181 // Destructor //
184  // do anything here that needs to be done at desctruction time
185  // (e.g. close files, deallocate resources etc.)
186 
187  delete rootFile_;
188 }
189 
191 // method called once each job just before starting event loop //
194  // go to *OUR* rootfile
195  rootFile_->cd();
196 
197  // Book Histograms
198  // PhotonID Histograms
199  h_isoEcalRecHit_ = new TH1F("photonEcalIso", "Ecal Rec Hit Isolation", 300, 0, 300);
200  h_isoHcalRecHit_ = new TH1F("photonHcalIso", "Hcal Rec Hit Isolation", 300, 0, 300);
201  h_trk_pt_solid_ = new TH1F("photonTrackSolidIso", "Sum of track pT in a cone of #DeltaR", 300, 0, 300);
202  h_trk_pt_hollow_ = new TH1F("photonTrackHollowIso", "Sum of track pT in a hollow cone", 300, 0, 300);
203  h_ntrk_solid_ = new TH1F("photonTrackCountSolid", "Number of tracks in a cone of #DeltaR", 100, 0, 100);
204  h_ntrk_hollow_ = new TH1F("photonTrackCountHollow", "Number of tracks in a hollow cone", 100, 0, 100);
205  h_ebetagap_ = new TH1F("photonInEBEtagap", "Ecal Barrel eta gap flag", 2, -0.5, 1.5);
206  h_ebphigap_ = new TH1F("photonInEBEtagap", "Ecal Barrel phi gap flag", 2, -0.5, 1.5);
207  h_eeringGap_ = new TH1F("photonInEERinggap", "Ecal Endcap ring gap flag", 2, -0.5, 1.5);
208  h_eedeeGap_ = new TH1F("photonInEEDeegap", "Ecal Endcap dee gap flag", 2, -0.5, 1.5);
209  h_ebeeGap_ = new TH1F("photonInEEgap", "Ecal Barrel/Endcap gap flag", 2, -0.5, 1.5);
210  h_r9_ = new TH1F("photonR9", "R9 = E(3x3) / E(SuperCluster)", 300, 0, 3);
211 
212  // Photon Histograms
213  h_photonEt_ = new TH1F("photonEt", "Photon E_{T}", 200, 0, 200);
214  h_photonEta_ = new TH1F("photonEta", "Photon #eta", 800, -4, 4);
215  h_photonPhi_ = new TH1F("photonPhi", "Photon #phi", 628, -1. * M_PI, M_PI);
216  h_hadoverem_ = new TH1F("photonHoverE", "Hadronic over EM", 200, 0, 1);
217 
218  // Photon's SuperCluster Histograms
219  h_photonScEt_ = new TH1F("photonScEt", "Photon SuperCluster E_{T}", 200, 0, 200);
220  h_photonScEta_ = new TH1F("photonScEta", "Photon #eta", 800, -4, 4);
221  h_photonScPhi_ = new TH1F("photonScPhi", "Photon #phi", 628, -1. * M_PI, M_PI);
222  h_photonScEtaWidth_ = new TH1F("photonScEtaWidth", "#eta-width", 100, 0, .1);
223 
224  // Composite or Other Histograms
225  h_photonInAnyGap_ = new TH1F("photonInAnyGap", "Photon in any gap flag", 2, -0.5, 1.5);
226  h_nPassingPho_ = new TH1F("photonLoosePhoton", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
227  h_nPassEM_ = new TH1F("photonLooseEM", "Total number photons (0=NotPassing, 1=Passing)", 2, -0.5, 1.5);
228  h_nPho_ = new TH1F("photonCount", "Number of photons passing cuts in event", 10, 0, 10);
229 
230  // Create a TTree of photons if set to 'True' in config file
231  if (createPhotonTTree_) {
232  tree_PhotonAll_ = new TTree("TreePhotonAll", "Reconstructed Photon");
233  tree_PhotonAll_->Branch("recPhoton",
234  &recPhoton.isolationEcalRecHit,
235  "isolationEcalRecHit/"
236  "F:isolationHcalRecHit:isolationSolidTrkCone:isolationHollowTrkCone:nTrkSolidCone:"
237  "nTrkHollowCone:isEBGap:isEEGap:isEBEEGap:r9:et:eta:phi:hadronicOverEm");
238  }
239 }
240 
242 // method called to for each event //
245  using namespace std;
246  using namespace edm;
247 
248  // grab photons
250  evt.getByLabel("photons", "", photonColl);
251 
252  Handle<edm::ValueMap<bool> > loosePhotonQual;
253  evt.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonQual);
254  Handle<edm::ValueMap<bool> > looseEMQual;
255  evt.getByLabel("PhotonIDProd", "PhotonCutBasedIDLooseEM", looseEMQual);
256  // grab PhotonId objects
257  // Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
258  // evt.getByLabel("PhotonIDProd", "PhotonAssociatedID", photonIDMapColl);
259 
260  // create reference to the object types we are interested in
261  const reco::PhotonCollection* photons = photonColl.product();
262  const edm::ValueMap<bool>* phoMap = loosePhotonQual.product();
263  const edm::ValueMap<bool>* lEMMap = looseEMQual.product();
264  int photonCounter = 0;
265  int idxpho = 0;
266  reco::PhotonCollection::const_iterator pho;
267  for (pho = (*photons).begin(); pho != (*photons).end(); pho++) {
268  edm::Ref<reco::PhotonCollection> photonref(photonColl, idxpho);
269  //reco::PhotonIDAssociationCollection::const_iterator photonIter = phoMap->find(photonref);
270  //const reco::PhotonIDRef &phtn = photonIter->val;
271  //const reco::PhotonRef &pho = photonIter->key;
272 
273  float photonEt = pho->et();
274  float superClusterEt = (pho->superCluster()->energy()) / (cosh(pho->superCluster()->position().eta()));
275  bool LoosePhotonQu = (*phoMap)[photonref];
276  h_nPassingPho_->Fill(LoosePhotonQu);
277  bool LooseEMQu = (*lEMMap)[photonref];
278  h_nPassEM_->Fill(LooseEMQu);
279  // Only store photon candidates (SuperClusters) that pass some simple cuts
280  bool passCuts = (photonEt > minPhotonEt_) && (fabs(pho->eta()) > minPhotonAbsEta_) &&
281  (fabs(pho->eta()) < maxPhotonAbsEta_) && (pho->r9() > minPhotonR9_) &&
282  (pho->hadronicOverEm() < maxPhotonHoverE_);
283 
284  if (passCuts) {
286  // fill histograms //
288  // PhotonID Variables
289  h_isoEcalRecHit_->Fill(pho->ecalRecHitSumEtConeDR04());
290  h_isoHcalRecHit_->Fill(pho->hcalTowerSumEtConeDR04());
291  h_trk_pt_solid_->Fill(pho->trkSumPtSolidConeDR04());
292  h_trk_pt_hollow_->Fill(pho->trkSumPtHollowConeDR04());
293  h_ntrk_solid_->Fill(pho->nTrkSolidConeDR04());
294  h_ntrk_hollow_->Fill(pho->nTrkHollowConeDR04());
295  h_ebetagap_->Fill(pho->isEBEtaGap());
296  h_ebphigap_->Fill(pho->isEBPhiGap());
297  h_eeringGap_->Fill(pho->isEERingGap());
298  h_eedeeGap_->Fill(pho->isEEDeeGap());
299  h_ebeeGap_->Fill(pho->isEBEEGap());
300  h_r9_->Fill(pho->r9());
301 
302  // Photon Variables
303  h_photonEt_->Fill(photonEt);
304  h_photonEta_->Fill(pho->eta());
305  h_photonPhi_->Fill(pho->phi());
306  h_hadoverem_->Fill(pho->hadronicOverEm());
307 
308  // Photon's SuperCluster Variables
309  // eta is with respect to detector (not physics) vertex,
310  // thus Et and eta are different from photon.
311  h_photonScEt_->Fill(superClusterEt);
312  h_photonScEta_->Fill(pho->superCluster()->position().eta());
313  h_photonScPhi_->Fill(pho->superCluster()->position().phi());
314  h_photonScEtaWidth_->Fill(pho->superCluster()->etaWidth());
315 
316  // It passed photon cuts, mark it
317 
319  // fill TTree (optional) //
321  if (createPhotonTTree_) {
322  recPhoton.isolationEcalRecHit = pho->ecalRecHitSumEtConeDR04();
323  recPhoton.isolationHcalRecHit = pho->hcalTowerSumEtConeDR04();
324  recPhoton.isolationSolidTrkCone = pho->trkSumPtSolidConeDR04();
325  recPhoton.isolationHollowTrkCone = pho->trkSumPtHollowConeDR04();
326  recPhoton.nTrkSolidCone = pho->nTrkSolidConeDR04();
327  recPhoton.nTrkHollowCone = pho->nTrkHollowConeDR04();
328  recPhoton.isEBEtaGap = pho->isEBEtaGap();
329  recPhoton.isEBPhiGap = pho->isEBPhiGap();
330  recPhoton.isEERingGap = pho->isEERingGap();
331  recPhoton.isEEDeeGap = pho->isEEDeeGap();
332  recPhoton.isEBEEGap = pho->isEBEEGap();
333  recPhoton.r9 = pho->r9();
334  recPhoton.et = pho->et();
335  recPhoton.eta = pho->eta();
336  recPhoton.phi = pho->phi();
337  recPhoton.hadronicOverEm = pho->hadronicOverEm();
338 
339  // Fill the tree (this records all the recPhoton.* since
340  // tree_PhotonAll_ was set to point at that.
341  tree_PhotonAll_->Fill();
342  }
343 
344  // Record whether it was near any module gap.
345  // Very convoluted at the moment.
346  bool inAnyGap = pho->isEBEEGap() || (pho->isEB() && pho->isEBEtaGap()) || (pho->isEB() && pho->isEBPhiGap()) ||
347  (pho->isEE() && pho->isEERingGap()) || (pho->isEE() && pho->isEEDeeGap());
348  if (inAnyGap) {
349  h_photonInAnyGap_->Fill(1.0);
350  } else {
351  h_photonInAnyGap_->Fill(0.0);
352  }
353 
354  photonCounter++;
355  } else {
356  // This didn't pass photon cuts, mark it
357  }
358  idxpho++;
359  } // End Loop over photons
360  h_nPho_->Fill(photonCounter);
361 }
362 
364 // method called once each job just after ending the event loop //
367  // go to *OUR* root file and store histograms
368  rootFile_->cd();
369 
370  // PhotonID Histograms
371  h_isoEcalRecHit_->Write();
372  h_isoHcalRecHit_->Write();
373  h_trk_pt_solid_->Write();
374  h_trk_pt_hollow_->Write();
375  h_ntrk_solid_->Write();
376  h_ntrk_hollow_->Write();
377  h_ebetagap_->Write();
378  h_ebphigap_->Write();
379  h_eeringGap_->Write();
380  h_eedeeGap_->Write();
381  h_ebeeGap_->Write();
382  h_r9_->Write();
383 
384  // Photon Histograms
385  h_photonEt_->Write();
386  h_photonEta_->Write();
387  h_photonPhi_->Write();
388  h_hadoverem_->Write();
389 
390  // Photon's SuperCluster Histograms
391  h_photonScEt_->Write();
392  h_photonScEta_->Write();
393  h_photonScPhi_->Write();
394  h_photonScEtaWidth_->Write();
395 
396  // Composite or Other Histograms
397  h_photonInAnyGap_->Write();
398  h_nPassingPho_->Write();
399  h_nPassEM_->Write();
400  h_nPho_->Write();
401 
402  // Write the root file (really writes the TTree)
403  rootFile_->Write();
404  rootFile_->Close();
405 }
406 
407 //define this as a plug-in
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
T const * product() const
Definition: Handle.h:70
PhotonIDSimpleAnalyzer(const edm::ParameterSet &)
void analyze(const edm::Event &, const edm::EventSetup &) override
#define M_PI
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
HLT enums.
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500