File indexing completed on 2024-04-06 12:25:09
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
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
0045
0046 #include "TH1F.h"
0047 #include "TH2F.h"
0048 #include "TFile.h"
0049 #include "TTree.h"
0050
0051
0052 #include <memory>
0053
0054
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
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_;
0083 double minPhotonEt_;
0084 double minPhotonAbsEta_;
0085 double maxPhotonAbsEta_;
0086 double minPhotonR9_;
0087 double maxPhotonHoverE_;
0088 bool createPhotonTTree_;
0089
0090
0091
0092
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
0114 TFile* rootFile_;
0115
0116
0117
0118
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
0133 TH1F* h_photonEt_;
0134 TH1F* h_photonEta_;
0135 TH1F* h_photonPhi_;
0136 TH1F* h_hadoverem_;
0137
0138
0139 TH1F* h_photonScEt_;
0140 TH1F* h_photonScEta_;
0141 TH1F* h_photonScPhi_;
0142 TH1F* h_photonScEtaWidth_;
0143
0144
0145 TH1F* h_photonInAnyGap_;
0146 TH1F* h_nPassingPho_;
0147 TH1F* h_nPassEM_;
0148 TH1F* h_nPho_;
0149
0150
0151 TTree* tree_PhotonAll_;
0152 };
0153
0154 using namespace std;
0155
0156
0157
0158
0159 PhotonIDSimpleAnalyzer::PhotonIDSimpleAnalyzer(const edm::ParameterSet& ps) {
0160
0161
0162
0163 outputFile_ = ps.getParameter<std::string>("outputFile");
0164
0165
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
0173
0174 createPhotonTTree_ = ps.getParameter<bool>("createPhotonTTree");
0175
0176
0177 rootFile_ = TFile::Open(outputFile_.c_str(), "RECREATE");
0178 }
0179
0180
0181
0182
0183 PhotonIDSimpleAnalyzer::~PhotonIDSimpleAnalyzer() {
0184
0185
0186
0187 delete rootFile_;
0188 }
0189
0190
0191
0192
0193 void PhotonIDSimpleAnalyzer::beginJob() {
0194
0195 rootFile_->cd();
0196
0197
0198
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
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
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
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
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
0243
0244 void PhotonIDSimpleAnalyzer::analyze(const edm::Event& evt, const edm::EventSetup& es) {
0245 using namespace std;
0246 using namespace edm;
0247
0248
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
0257
0258
0259
0260
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
0270
0271
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
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
0287
0288
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
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
0309
0310
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
0317
0318
0319
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
0340
0341 tree_PhotonAll_->Fill();
0342 }
0343
0344
0345
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
0357 }
0358 idxpho++;
0359 }
0360 h_nPho_->Fill(photonCounter);
0361 }
0362
0363
0364
0365
0366 void PhotonIDSimpleAnalyzer::endJob() {
0367
0368 rootFile_->cd();
0369
0370
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
0385 h_photonEt_->Write();
0386 h_photonEta_->Write();
0387 h_photonPhi_->Write();
0388 h_hadoverem_->Write();
0389
0390
0391 h_photonScEt_->Write();
0392 h_photonScEta_->Write();
0393 h_photonScPhi_->Write();
0394 h_photonScEtaWidth_->Write();
0395
0396
0397 h_photonInAnyGap_->Write();
0398 h_nPassingPho_->Write();
0399 h_nPassEM_->Write();
0400 h_nPho_->Write();
0401
0402
0403 rootFile_->Write();
0404 rootFile_->Close();
0405 }
0406
0407
0408 DEFINE_FWK_MODULE(PhotonIDSimpleAnalyzer);