File indexing completed on 2023-03-17 11:21:17
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 ParameterSet conf_;
0060
0061 unsigned int ev;
0062
0063
0064 TH1F *h_etamc_ele, *h_etaec_ele, *h_etaeg_ele, *h_etaged_ele;
0065 TH1F *h_ptmc_ele, *h_ptec_ele, *h_pteg_ele, *h_ptged_ele;
0066
0067 TH2F *etaphi_nonreco;
0068
0069 TH1F *eg_EEcalEtrue_1, *eg_EEcalEtrue_2, *eg_EEcalEtrue_3, *eg_EEcalEtrue_4, *eg_EEcalEtrue_5;
0070 TH1F *pf_EEcalEtrue_1, *pf_EEcalEtrue_2, *pf_EEcalEtrue_3, *pf_EEcalEtrue_4, *pf_EEcalEtrue_5;
0071 TH1F *ged_EEcalEtrue_1, *ged_EEcalEtrue_2, *ged_EEcalEtrue_3, *ged_EEcalEtrue_4, *ged_EEcalEtrue_5;
0072
0073 TH1F *eg_hoe, *pf_hoe, *ged_hoe;
0074
0075 TH1F *eg_sigmaetaeta_eb, *pf_sigmaetaeta_eb, *ged_sigmaetaeta_eb;
0076 TH1F *eg_sigmaetaeta_ee, *pf_sigmaetaeta_ee, *ged_sigmaetaeta_ee;
0077
0078 TH1F *eg_r9_eb, *pf_r9_eb, *ged_r9_eb;
0079 TH1F *eg_r9_ee, *pf_r9_ee, *ged_r9_ee;
0080 };
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093 EgGEDPhotonAnalyzer::EgGEDPhotonAnalyzer(const edm::ParameterSet &iConfig)
0094 : conf_(iConfig)
0095
0096 {
0097 usesResource(TFileService::kSharedResource);
0098
0099 edm::Service<TFileService> fs;
0100
0101
0102
0103 h_etamc_ele = fs->make<TH1F>("h_etamc_ele", " ", 50, -2.5, 2.5);
0104 h_etaec_ele = fs->make<TH1F>("h_etaec_ele", " ", 50, -2.5, 2.5);
0105 h_etaeg_ele = fs->make<TH1F>("h_etaeg_ele", " ", 50, -2.5, 2.5);
0106 h_etaged_ele = fs->make<TH1F>("h_etaged_ele", " ", 50, -2.5, 2.5);
0107
0108 h_ptmc_ele = fs->make<TH1F>("h_ptmc_ele", " ", 50, 0, 100.);
0109 h_ptec_ele = fs->make<TH1F>("h_ptec_ele", " ", 50, 0, 100.);
0110 h_pteg_ele = fs->make<TH1F>("h_pteg_ele", " ", 50, 0, 100.);
0111
0112 h_ptged_ele = fs->make<TH1F>("h_ptged_ele", " ", 50, 0, 100.);
0113
0114 pf_hoe = fs->make<TH1F>("pf_hoe", " ", 100, 0, 1.);
0115 eg_hoe = fs->make<TH1F>("eg_hoe", " ", 100, 0, 1.);
0116 ged_hoe = fs->make<TH1F>("ged_hoe", " ", 100, 0, 1.);
0117
0118 pf_sigmaetaeta_eb = fs->make<TH1F>("pf_sigmaetaeta_eb", " ", 100, 0., 0.03);
0119 eg_sigmaetaeta_eb = fs->make<TH1F>("eg_sigmaetaeta_eb", " ", 100, 0., 0.03);
0120 ged_sigmaetaeta_eb = fs->make<TH1F>("ged_sigmaetaeta_eb", " ", 100, 0., 0.03);
0121
0122 pf_sigmaetaeta_ee = fs->make<TH1F>("pf_sigmaetaeta_ee", " ", 100, 0., 0.1);
0123 eg_sigmaetaeta_ee = fs->make<TH1F>("eg_sigmaetaeta_ee", " ", 100, 0., 0.1);
0124 ged_sigmaetaeta_ee = fs->make<TH1F>("ged_sigmaetaeta_ee", " ", 100, 0., 0.1);
0125
0126 pf_r9_eb = fs->make<TH1F>("pf_r9_eb", " ", 100, 0., 1.0);
0127 eg_r9_eb = fs->make<TH1F>("eg_r9_eb", " ", 100, 0., 1.0);
0128 ged_r9_eb = fs->make<TH1F>("ged_r9_eb", " ", 100, 0., 1.0);
0129
0130 pf_r9_ee = fs->make<TH1F>("pf_r9_ee", " ", 100, 0., 1.0);
0131 eg_r9_ee = fs->make<TH1F>("eg_r9_ee", " ", 100, 0., 1.0);
0132 ged_r9_ee = fs->make<TH1F>("ged_r9_ee", " ", 100, 0., 1.0);
0133
0134
0135 etaphi_nonreco = fs->make<TH2F>("etaphi_nonreco", " ", 50, -2.5, 2.5, 50, -3.2, 3.2);
0136
0137
0138 eg_EEcalEtrue_1 = fs->make<TH1F>("eg_EEcalEtrue_1", " ", 50, 0., 2.);
0139 eg_EEcalEtrue_2 = fs->make<TH1F>("eg_EEcalEtrue_2", " ", 50, 0., 2.);
0140 eg_EEcalEtrue_3 = fs->make<TH1F>("eg_EEcalEtrue_3", " ", 50, 0., 2.);
0141 eg_EEcalEtrue_4 = fs->make<TH1F>("eg_EEcalEtrue_4", " ", 50, 0., 2.);
0142 eg_EEcalEtrue_5 = fs->make<TH1F>("eg_EEcalEtrue_5", " ", 50, 0., 2.);
0143
0144 pf_EEcalEtrue_1 = fs->make<TH1F>("pf_EEcalEtrue_1", " ", 50, 0., 2.);
0145 pf_EEcalEtrue_2 = fs->make<TH1F>("pf_EEcalEtrue_2", " ", 50, 0., 2.);
0146 pf_EEcalEtrue_3 = fs->make<TH1F>("pf_EEcalEtrue_3", " ", 50, 0., 2.);
0147 pf_EEcalEtrue_4 = fs->make<TH1F>("pf_EEcalEtrue_4", " ", 50, 0., 2.);
0148 pf_EEcalEtrue_5 = fs->make<TH1F>("pf_EEcalEtrue_5", " ", 50, 0., 2.);
0149
0150 ged_EEcalEtrue_1 = fs->make<TH1F>("ged_EEcalEtrue_1", " ", 50, 0., 2.);
0151 ged_EEcalEtrue_2 = fs->make<TH1F>("ged_EEcalEtrue_2", " ", 50, 0., 2.);
0152 ged_EEcalEtrue_3 = fs->make<TH1F>("ged_EEcalEtrue_3", " ", 50, 0., 2.);
0153 ged_EEcalEtrue_4 = fs->make<TH1F>("ged_EEcalEtrue_4", " ", 50, 0., 2.);
0154 ged_EEcalEtrue_5 = fs->make<TH1F>("ged_EEcalEtrue_5", " ", 50, 0., 2.);
0155 }
0156
0157 EgGEDPhotonAnalyzer::~EgGEDPhotonAnalyzer() {
0158
0159
0160 }
0161
0162
0163
0164
0165
0166
0167 void EgGEDPhotonAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0168
0169 Handle<reco::PFCandidateCollection> collection;
0170 InputTag label("particleFlow");
0171 iEvent.getByLabel(label, collection);
0172 std::vector<reco::PFCandidate> candidates = (*collection.product());
0173
0174 InputTag recoPhotonLabel(string("photons"));
0175 Handle<PhotonCollection> theRecoPhotonCollection;
0176 iEvent.getByLabel(recoPhotonLabel, theRecoPhotonCollection);
0177 const PhotonCollection theRecoPh = *(theRecoPhotonCollection.product());
0178
0179 InputTag MCTruthCollection(string("VtxSmeared"));
0180 edm::Handle<edm::HepMCProduct> pMCTruth;
0181 iEvent.getByLabel(MCTruthCollection, pMCTruth);
0182 const HepMC::GenEvent *genEvent = pMCTruth->GetEvent();
0183
0184 InputTag gedPhotonLabel(string("gedPhotons"));
0185 Handle<PhotonCollection> theGedPhotonCollection;
0186 iEvent.getByLabel(gedPhotonLabel, theGedPhotonCollection);
0187 const PhotonCollection theGedPh = *(theGedPhotonCollection.product());
0188
0189 bool debug = true;
0190
0191 ev++;
0192
0193 if (debug)
0194 cout << "************************* New Event:: " << ev << " *************************" << endl;
0195
0196
0197
0198 for (HepMC::GenEvent::particle_const_iterator cP = genEvent->particles_begin(); cP != genEvent->particles_end();
0199 cP++) {
0200 float etamc = (*cP)->momentum().eta();
0201 float phimc = (*cP)->momentum().phi();
0202 float ptmc = (*cP)->momentum().perp();
0203 float Emc = (*cP)->momentum().e();
0204
0205 if (abs((*cP)->pdg_id()) == 22 && (*cP)->status() == 1 && (*cP)->momentum().perp() > 8. &&
0206 fabs((*cP)->momentum().eta()) < 2.5) {
0207 h_etamc_ele->Fill(etamc);
0208 h_ptmc_ele->Fill(ptmc);
0209
0210 float MindR = 1000.;
0211
0212 if (debug)
0213 cout << " MC particle: pt " << ptmc << " eta,phi " << etamc << ", " << phimc << endl;
0214
0215 std::vector<reco::PFCandidate>::iterator it;
0216 for (it = candidates.begin(); it != candidates.end(); ++it) {
0217 reco::PFCandidate::ParticleType type = (*it).particleId();
0218
0219 if (type == reco::PFCandidate::gamma) {
0220 reco::PhotonRef phref = (*it).photonRef();
0221 float eta = (*it).eta();
0222 float phi = (*it).phi();
0223
0224 float deta = etamc - eta;
0225 float dphi = normalizedPhi(phimc - phi);
0226 float dR = sqrt(deta * deta + dphi * dphi);
0227
0228 if (dR < 0.1) {
0229 MindR = dR;
0230
0231 if (phref.isNonnull()) {
0232 float SCPF = phref->energy();
0233 float EpfEtrue = SCPF / Emc;
0234 if (fabs(etamc) < 0.5) {
0235 pf_EEcalEtrue_1->Fill(EpfEtrue);
0236 }
0237 if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0238 pf_EEcalEtrue_2->Fill(EpfEtrue);
0239 }
0240 if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0241 pf_EEcalEtrue_3->Fill(EpfEtrue);
0242 }
0243 if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0244 pf_EEcalEtrue_4->Fill(EpfEtrue);
0245 }
0246 if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0247 pf_EEcalEtrue_5->Fill(EpfEtrue);
0248 }
0249
0250 pf_hoe->Fill(phref->hadronicOverEm());
0251 if (fabs(etamc) < 1.479) {
0252 pf_r9_eb->Fill(phref->r9());
0253 pf_sigmaetaeta_eb->Fill(phref->sigmaEtaEta());
0254 } else if (fabs(etamc) > 1.479 && fabs(etamc) < 2.5) {
0255 pf_r9_ee->Fill(phref->r9());
0256 pf_sigmaetaeta_ee->Fill(phref->sigmaEtaEta());
0257 }
0258 }
0259
0260 if (debug)
0261 cout << " PF photon matched: pt " << (*it).pt() << " eta,phi " << eta << ", " << phi << endl;
0262
0263 }
0264 }
0265 }
0266
0267 if (MindR < 0.1) {
0268 h_etaec_ele->Fill(etamc);
0269 h_ptec_ele->Fill(ptmc);
0270
0271 } else {
0272 etaphi_nonreco->Fill(etamc, phimc);
0273 }
0274
0275 float MindREG = 1000;
0276 for (uint j = 0; j < theRecoPh.size(); j++) {
0277 float etareco = theRecoPh[j].eta();
0278 float phireco = theRecoPh[j].phi();
0279
0280 float deta = etamc - etareco;
0281 float dphi = normalizedPhi(phimc - phireco);
0282 float dR = sqrt(deta * deta + dphi * dphi);
0283
0284 float SCEnergy = (theRecoPh[j]).energy();
0285 float ErecoEtrue = SCEnergy / Emc;
0286
0287 if (dR < 0.1) {
0288 if (debug)
0289 cout << " EG ele matched: pt " << theRecoPh[j].pt() << " eta,phi " << etareco << ", " << phireco << endl;
0290
0291 if (fabs(etamc) < 0.5) {
0292 eg_EEcalEtrue_1->Fill(ErecoEtrue);
0293 }
0294 if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0295 eg_EEcalEtrue_2->Fill(ErecoEtrue);
0296 }
0297 if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0298 eg_EEcalEtrue_3->Fill(ErecoEtrue);
0299 }
0300 if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0301 eg_EEcalEtrue_4->Fill(ErecoEtrue);
0302 }
0303 if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0304 eg_EEcalEtrue_5->Fill(ErecoEtrue);
0305 }
0306
0307 eg_hoe->Fill(theRecoPh[j].hadronicOverEm());
0308 if (fabs(etamc) < 1.479) {
0309 eg_r9_eb->Fill(theRecoPh[j].r9());
0310 eg_sigmaetaeta_eb->Fill(theRecoPh[j].sigmaEtaEta());
0311 } else if (fabs(etamc) > 1.479 && fabs(etamc) < 2.5) {
0312 eg_r9_ee->Fill(theRecoPh[j].r9());
0313 eg_sigmaetaeta_ee->Fill(theRecoPh[j].sigmaEtaEta());
0314 }
0315
0316 MindREG = dR;
0317 }
0318 }
0319 if (MindREG < 0.1) {
0320 h_etaeg_ele->Fill(etamc);
0321 h_pteg_ele->Fill(ptmc);
0322 }
0323
0324
0325
0326 float MindRGedEg = 1000;
0327
0328 for (uint j = 0; j < theGedPh.size(); j++) {
0329 reco::GsfTrackRef egGsfTrackRef = (theGedPh[j]).gsfTrack();
0330 float etareco = theGedPh[j].eta();
0331 float phireco = theGedPh[j].phi();
0332
0333 float deta = etamc - etareco;
0334 float dphi = normalizedPhi(phimc - phireco);
0335 float dR = sqrt(deta * deta + dphi * dphi);
0336
0337 float SCEnergy = (theGedPh[j]).energy();
0338 float ErecoEtrue = SCEnergy / Emc;
0339
0340 if (dR < 0.1) {
0341 if (debug)
0342 cout << " GED ele matched: pt " << theGedPh[j].pt() << " eta,phi " << etareco << ", " << phireco << endl;
0343
0344 if (fabs(etamc) < 0.5) {
0345 ged_EEcalEtrue_1->Fill(ErecoEtrue);
0346 }
0347 if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0348 ged_EEcalEtrue_2->Fill(ErecoEtrue);
0349 }
0350 if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0351 ged_EEcalEtrue_3->Fill(ErecoEtrue);
0352 }
0353 if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0354 ged_EEcalEtrue_4->Fill(ErecoEtrue);
0355 }
0356 if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0357 ged_EEcalEtrue_5->Fill(ErecoEtrue);
0358 }
0359
0360 ged_hoe->Fill(theGedPh[j].hadronicOverEm());
0361 if (fabs(etamc) < 1.479) {
0362 ged_r9_eb->Fill(theGedPh[j].r9());
0363 ged_sigmaetaeta_eb->Fill(theGedPh[j].sigmaEtaEta());
0364 } else if (fabs(etamc) > 1.479 && fabs(etamc) < 2.5) {
0365 ged_r9_ee->Fill(theGedPh[j].r9());
0366 ged_sigmaetaeta_ee->Fill(theGedPh[j].sigmaEtaEta());
0367 }
0368
0369 MindRGedEg = dR;
0370 }
0371 }
0372 if (MindRGedEg < 0.1) {
0373 h_etaged_ele->Fill(etamc);
0374 h_ptged_ele->Fill(ptmc);
0375 }
0376
0377 }
0378
0379 }
0380 }
0381
0382 void EgGEDPhotonAnalyzer::beginJob() { ev = 0; }
0383
0384
0385 void EgGEDPhotonAnalyzer::endJob() { cout << " endJob:: #events " << ev << endl; }
0386
0387 DEFINE_FWK_MODULE(EgGEDPhotonAnalyzer);