File indexing completed on 2024-09-07 04:37:52
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <memory>
0018
0019
0020 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ServiceRegistry/interface/Service.h"
0025 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0026 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0027 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0028 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0029 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0030 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0031 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0033 #include "DataFormats/Math/interface/normalizedPhi.h"
0034
0035 #include <vector>
0036 #include <TROOT.h>
0037 #include <TFile.h>
0038 #include <TTree.h>
0039 #include <TH1F.h>
0040 #include <TH2F.h>
0041 #include <TLorentzVector.h>
0042
0043
0044
0045
0046 using namespace edm;
0047 using namespace reco;
0048 using namespace std;
0049 class EgGEDPhotonAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0050 public:
0051 explicit EgGEDPhotonAnalyzer(const edm::ParameterSet &);
0052 ~EgGEDPhotonAnalyzer() override;
0053
0054 private:
0055 void beginJob() override;
0056 void analyze(const edm::Event &, const edm::EventSetup &) override;
0057 void endJob() override;
0058
0059 unsigned int ev;
0060
0061
0062 TH1F *h_etamc_ele, *h_etaec_ele, *h_etaeg_ele, *h_etaged_ele;
0063 TH1F *h_ptmc_ele, *h_ptec_ele, *h_pteg_ele, *h_ptged_ele;
0064
0065 TH2F *etaphi_nonreco;
0066
0067 TH1F *eg_EEcalEtrue_1, *eg_EEcalEtrue_2, *eg_EEcalEtrue_3, *eg_EEcalEtrue_4, *eg_EEcalEtrue_5;
0068 TH1F *pf_EEcalEtrue_1, *pf_EEcalEtrue_2, *pf_EEcalEtrue_3, *pf_EEcalEtrue_4, *pf_EEcalEtrue_5;
0069 TH1F *ged_EEcalEtrue_1, *ged_EEcalEtrue_2, *ged_EEcalEtrue_3, *ged_EEcalEtrue_4, *ged_EEcalEtrue_5;
0070
0071 TH1F *eg_hoe, *pf_hoe, *ged_hoe;
0072
0073 TH1F *eg_sigmaetaeta_eb, *pf_sigmaetaeta_eb, *ged_sigmaetaeta_eb;
0074 TH1F *eg_sigmaetaeta_ee, *pf_sigmaetaeta_ee, *ged_sigmaetaeta_ee;
0075
0076 TH1F *eg_r9_eb, *pf_r9_eb, *ged_r9_eb;
0077 TH1F *eg_r9_ee, *pf_r9_ee, *ged_r9_ee;
0078
0079 const edm::EDGetTokenT<reco::PFCandidateCollection> pfToken_;
0080 const edm::EDGetTokenT<reco::PhotonCollection> phoToken_;
0081 const edm::EDGetTokenT<edm::HepMCProduct> mcToken_;
0082 const edm::EDGetTokenT<reco::PhotonCollection> gedPhoToken_;
0083 };
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096 EgGEDPhotonAnalyzer::EgGEDPhotonAnalyzer(const edm::ParameterSet &iConfig)
0097 : pfToken_(consumes<reco::PFCandidateCollection>(edm::InputTag("particleFlow"))),
0098 phoToken_(consumes<reco::PhotonCollection>(edm::InputTag("photons"))),
0099 mcToken_(consumes<edm::HepMCProduct>(edm::InputTag("VtxSmeared"))),
0100 gedPhoToken_(consumes<reco::PhotonCollection>(edm::InputTag("gedPhotons"))) {
0101 usesResource(TFileService::kSharedResource);
0102
0103 edm::Service<TFileService> fs;
0104
0105
0106
0107 h_etamc_ele = fs->make<TH1F>("h_etamc_ele", " ", 50, -2.5, 2.5);
0108 h_etaec_ele = fs->make<TH1F>("h_etaec_ele", " ", 50, -2.5, 2.5);
0109 h_etaeg_ele = fs->make<TH1F>("h_etaeg_ele", " ", 50, -2.5, 2.5);
0110 h_etaged_ele = fs->make<TH1F>("h_etaged_ele", " ", 50, -2.5, 2.5);
0111
0112 h_ptmc_ele = fs->make<TH1F>("h_ptmc_ele", " ", 50, 0, 100.);
0113 h_ptec_ele = fs->make<TH1F>("h_ptec_ele", " ", 50, 0, 100.);
0114 h_pteg_ele = fs->make<TH1F>("h_pteg_ele", " ", 50, 0, 100.);
0115
0116 h_ptged_ele = fs->make<TH1F>("h_ptged_ele", " ", 50, 0, 100.);
0117
0118 pf_hoe = fs->make<TH1F>("pf_hoe", " ", 100, 0, 1.);
0119 eg_hoe = fs->make<TH1F>("eg_hoe", " ", 100, 0, 1.);
0120 ged_hoe = fs->make<TH1F>("ged_hoe", " ", 100, 0, 1.);
0121
0122 pf_sigmaetaeta_eb = fs->make<TH1F>("pf_sigmaetaeta_eb", " ", 100, 0., 0.03);
0123 eg_sigmaetaeta_eb = fs->make<TH1F>("eg_sigmaetaeta_eb", " ", 100, 0., 0.03);
0124 ged_sigmaetaeta_eb = fs->make<TH1F>("ged_sigmaetaeta_eb", " ", 100, 0., 0.03);
0125
0126 pf_sigmaetaeta_ee = fs->make<TH1F>("pf_sigmaetaeta_ee", " ", 100, 0., 0.1);
0127 eg_sigmaetaeta_ee = fs->make<TH1F>("eg_sigmaetaeta_ee", " ", 100, 0., 0.1);
0128 ged_sigmaetaeta_ee = fs->make<TH1F>("ged_sigmaetaeta_ee", " ", 100, 0., 0.1);
0129
0130 pf_r9_eb = fs->make<TH1F>("pf_r9_eb", " ", 100, 0., 1.0);
0131 eg_r9_eb = fs->make<TH1F>("eg_r9_eb", " ", 100, 0., 1.0);
0132 ged_r9_eb = fs->make<TH1F>("ged_r9_eb", " ", 100, 0., 1.0);
0133
0134 pf_r9_ee = fs->make<TH1F>("pf_r9_ee", " ", 100, 0., 1.0);
0135 eg_r9_ee = fs->make<TH1F>("eg_r9_ee", " ", 100, 0., 1.0);
0136 ged_r9_ee = fs->make<TH1F>("ged_r9_ee", " ", 100, 0., 1.0);
0137
0138
0139 etaphi_nonreco = fs->make<TH2F>("etaphi_nonreco", " ", 50, -2.5, 2.5, 50, -3.2, 3.2);
0140
0141
0142 eg_EEcalEtrue_1 = fs->make<TH1F>("eg_EEcalEtrue_1", " ", 50, 0., 2.);
0143 eg_EEcalEtrue_2 = fs->make<TH1F>("eg_EEcalEtrue_2", " ", 50, 0., 2.);
0144 eg_EEcalEtrue_3 = fs->make<TH1F>("eg_EEcalEtrue_3", " ", 50, 0., 2.);
0145 eg_EEcalEtrue_4 = fs->make<TH1F>("eg_EEcalEtrue_4", " ", 50, 0., 2.);
0146 eg_EEcalEtrue_5 = fs->make<TH1F>("eg_EEcalEtrue_5", " ", 50, 0., 2.);
0147
0148 pf_EEcalEtrue_1 = fs->make<TH1F>("pf_EEcalEtrue_1", " ", 50, 0., 2.);
0149 pf_EEcalEtrue_2 = fs->make<TH1F>("pf_EEcalEtrue_2", " ", 50, 0., 2.);
0150 pf_EEcalEtrue_3 = fs->make<TH1F>("pf_EEcalEtrue_3", " ", 50, 0., 2.);
0151 pf_EEcalEtrue_4 = fs->make<TH1F>("pf_EEcalEtrue_4", " ", 50, 0., 2.);
0152 pf_EEcalEtrue_5 = fs->make<TH1F>("pf_EEcalEtrue_5", " ", 50, 0., 2.);
0153
0154 ged_EEcalEtrue_1 = fs->make<TH1F>("ged_EEcalEtrue_1", " ", 50, 0., 2.);
0155 ged_EEcalEtrue_2 = fs->make<TH1F>("ged_EEcalEtrue_2", " ", 50, 0., 2.);
0156 ged_EEcalEtrue_3 = fs->make<TH1F>("ged_EEcalEtrue_3", " ", 50, 0., 2.);
0157 ged_EEcalEtrue_4 = fs->make<TH1F>("ged_EEcalEtrue_4", " ", 50, 0., 2.);
0158 ged_EEcalEtrue_5 = fs->make<TH1F>("ged_EEcalEtrue_5", " ", 50, 0., 2.);
0159 }
0160
0161 EgGEDPhotonAnalyzer::~EgGEDPhotonAnalyzer() {
0162
0163
0164 }
0165
0166
0167
0168
0169
0170
0171 void EgGEDPhotonAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0172
0173 const auto &candidates = iEvent.get(pfToken_);
0174 const auto &theRecoPh = iEvent.get(phoToken_);
0175 const auto &genEvent = iEvent.get(mcToken_).GetEvent();
0176 const auto &theGedPh = iEvent.get(gedPhoToken_);
0177
0178 bool debug = true;
0179
0180 ev++;
0181
0182 if (debug)
0183 cout << "************************* New Event:: " << ev << " *************************" << endl;
0184
0185
0186
0187 for (HepMC::GenEvent::particle_const_iterator cP = genEvent->particles_begin(); cP != genEvent->particles_end();
0188 cP++) {
0189 float etamc = (*cP)->momentum().eta();
0190 float phimc = (*cP)->momentum().phi();
0191 float ptmc = (*cP)->momentum().perp();
0192 float Emc = (*cP)->momentum().e();
0193
0194 if (abs((*cP)->pdg_id()) == 22 && (*cP)->status() == 1 && (*cP)->momentum().perp() > 8. &&
0195 fabs((*cP)->momentum().eta()) < 2.5) {
0196 h_etamc_ele->Fill(etamc);
0197 h_ptmc_ele->Fill(ptmc);
0198
0199 float MindR = 1000.;
0200
0201 if (debug)
0202 cout << " MC particle: pt " << ptmc << " eta,phi " << etamc << ", " << phimc << endl;
0203
0204 for (const auto &cand : candidates) {
0205 reco::PFCandidate::ParticleType type = cand.particleId();
0206
0207 if (type == reco::PFCandidate::gamma) {
0208 reco::PhotonRef phref = cand.photonRef();
0209 float eta = cand.eta();
0210 float phi = cand.phi();
0211
0212 float deta = etamc - eta;
0213 float dphi = normalizedPhi(phimc - phi);
0214 float dR = sqrt(deta * deta + dphi * dphi);
0215
0216 if (dR < 0.1) {
0217 MindR = dR;
0218
0219 if (phref.isNonnull()) {
0220 float SCPF = phref->energy();
0221 float EpfEtrue = SCPF / Emc;
0222 if (fabs(etamc) < 0.5) {
0223 pf_EEcalEtrue_1->Fill(EpfEtrue);
0224 }
0225 if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0226 pf_EEcalEtrue_2->Fill(EpfEtrue);
0227 }
0228 if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0229 pf_EEcalEtrue_3->Fill(EpfEtrue);
0230 }
0231 if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0232 pf_EEcalEtrue_4->Fill(EpfEtrue);
0233 }
0234 if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0235 pf_EEcalEtrue_5->Fill(EpfEtrue);
0236 }
0237
0238 pf_hoe->Fill(phref->hadronicOverEm());
0239 if (fabs(etamc) < 1.479) {
0240 pf_r9_eb->Fill(phref->r9());
0241 pf_sigmaetaeta_eb->Fill(phref->sigmaEtaEta());
0242 } else if (fabs(etamc) > 1.479 && fabs(etamc) < 2.5) {
0243 pf_r9_ee->Fill(phref->r9());
0244 pf_sigmaetaeta_ee->Fill(phref->sigmaEtaEta());
0245 }
0246 }
0247
0248 if (debug)
0249 cout << " PF photon matched: pt " << cand.pt() << " eta,phi " << eta << ", " << phi << endl;
0250
0251 }
0252 }
0253 }
0254
0255 if (MindR < 0.1) {
0256 h_etaec_ele->Fill(etamc);
0257 h_ptec_ele->Fill(ptmc);
0258
0259 } else {
0260 etaphi_nonreco->Fill(etamc, phimc);
0261 }
0262
0263 float MindREG = 1000;
0264 for (uint j = 0; j < theRecoPh.size(); j++) {
0265 float etareco = theRecoPh[j].eta();
0266 float phireco = theRecoPh[j].phi();
0267
0268 float deta = etamc - etareco;
0269 float dphi = normalizedPhi(phimc - phireco);
0270 float dR = sqrt(deta * deta + dphi * dphi);
0271
0272 float SCEnergy = (theRecoPh[j]).energy();
0273 float ErecoEtrue = SCEnergy / Emc;
0274
0275 if (dR < 0.1) {
0276 if (debug)
0277 cout << " EG ele matched: pt " << theRecoPh[j].pt() << " eta,phi " << etareco << ", " << phireco << endl;
0278
0279 if (fabs(etamc) < 0.5) {
0280 eg_EEcalEtrue_1->Fill(ErecoEtrue);
0281 }
0282 if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0283 eg_EEcalEtrue_2->Fill(ErecoEtrue);
0284 }
0285 if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0286 eg_EEcalEtrue_3->Fill(ErecoEtrue);
0287 }
0288 if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0289 eg_EEcalEtrue_4->Fill(ErecoEtrue);
0290 }
0291 if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0292 eg_EEcalEtrue_5->Fill(ErecoEtrue);
0293 }
0294
0295 eg_hoe->Fill(theRecoPh[j].hadronicOverEm());
0296 if (fabs(etamc) < 1.479) {
0297 eg_r9_eb->Fill(theRecoPh[j].r9());
0298 eg_sigmaetaeta_eb->Fill(theRecoPh[j].sigmaEtaEta());
0299 } else if (fabs(etamc) > 1.479 && fabs(etamc) < 2.5) {
0300 eg_r9_ee->Fill(theRecoPh[j].r9());
0301 eg_sigmaetaeta_ee->Fill(theRecoPh[j].sigmaEtaEta());
0302 }
0303
0304 MindREG = dR;
0305 }
0306 }
0307 if (MindREG < 0.1) {
0308 h_etaeg_ele->Fill(etamc);
0309 h_pteg_ele->Fill(ptmc);
0310 }
0311
0312
0313
0314 float MindRGedEg = 1000;
0315
0316 for (uint j = 0; j < theGedPh.size(); j++) {
0317 reco::GsfTrackRef egGsfTrackRef = (theGedPh[j]).gsfTrack();
0318 float etareco = theGedPh[j].eta();
0319 float phireco = theGedPh[j].phi();
0320
0321 float deta = etamc - etareco;
0322 float dphi = normalizedPhi(phimc - phireco);
0323 float dR = sqrt(deta * deta + dphi * dphi);
0324
0325 float SCEnergy = (theGedPh[j]).energy();
0326 float ErecoEtrue = SCEnergy / Emc;
0327
0328 if (dR < 0.1) {
0329 if (debug)
0330 cout << " GED ele matched: pt " << theGedPh[j].pt() << " eta,phi " << etareco << ", " << phireco << endl;
0331
0332 if (fabs(etamc) < 0.5) {
0333 ged_EEcalEtrue_1->Fill(ErecoEtrue);
0334 }
0335 if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0336 ged_EEcalEtrue_2->Fill(ErecoEtrue);
0337 }
0338 if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0339 ged_EEcalEtrue_3->Fill(ErecoEtrue);
0340 }
0341 if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0342 ged_EEcalEtrue_4->Fill(ErecoEtrue);
0343 }
0344 if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0345 ged_EEcalEtrue_5->Fill(ErecoEtrue);
0346 }
0347
0348 ged_hoe->Fill(theGedPh[j].hadronicOverEm());
0349 if (fabs(etamc) < 1.479) {
0350 ged_r9_eb->Fill(theGedPh[j].r9());
0351 ged_sigmaetaeta_eb->Fill(theGedPh[j].sigmaEtaEta());
0352 } else if (fabs(etamc) > 1.479 && fabs(etamc) < 2.5) {
0353 ged_r9_ee->Fill(theGedPh[j].r9());
0354 ged_sigmaetaeta_ee->Fill(theGedPh[j].sigmaEtaEta());
0355 }
0356
0357 MindRGedEg = dR;
0358 }
0359 }
0360 if (MindRGedEg < 0.1) {
0361 h_etaged_ele->Fill(etamc);
0362 h_ptged_ele->Fill(ptmc);
0363 }
0364
0365 }
0366
0367 }
0368 }
0369
0370 void EgGEDPhotonAnalyzer::beginJob() { ev = 0; }
0371
0372
0373 void EgGEDPhotonAnalyzer::endJob() { cout << " endJob:: #events " << ev << endl; }
0374
0375 DEFINE_FWK_MODULE(EgGEDPhotonAnalyzer);