Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:52 

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgGEDPhotonAnalyzer
0004 // Class:      EgGEDPhotonAnalyzer
0005 //
0006 /**\class EgGEDPhotonAnalyzer
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Daniele Benedetti
0015 
0016 // system include files
0017 #include <memory>
0018 
0019 // user include files
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 // class decleration
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   // ----------member data ---------------------------
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 // constants, enums and typedefs
0087 //
0088 
0089 //
0090 // static data member definitions
0091 //
0092 
0093 //
0094 // constructors and destructor
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   // Efficiency plots
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   // Ecal map
0139   etaphi_nonreco = fs->make<TH2F>("etaphi_nonreco", " ", 50, -2.5, 2.5, 50, -3.2, 3.2);
0140 
0141   // Energies
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   // do anything here that needs to be done at desctruction time
0163   // (e.g. close files, deallocate resources etc.)
0164 }
0165 
0166 //
0167 // member functions
0168 //
0169 
0170 // ------------ method called to for each event  ------------
0171 void EgGEDPhotonAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0172   // Candidate info
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   // Validation from generator events
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             // all for the moment
0251           }
0252         }  // End PFCandidates Electron Selection
0253       }  // End Loop PFCandidates
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       //Add loop on the GED electrons
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       }  //End Loop Generator Particles
0364 
0365     }  //End IF Generator Particles
0366 
0367   }  //End Loop Generator Particles
0368 }
0369 // ------------ method called once each job just before starting event loop  ------------
0370 void EgGEDPhotonAnalyzer::beginJob() { ev = 0; }
0371 
0372 // ------------ method called once each job just after ending the event loop  ------------
0373 void EgGEDPhotonAnalyzer::endJob() { cout << " endJob:: #events " << ev << endl; }
0374 //define this as a plug-in
0375 DEFINE_FWK_MODULE(EgGEDPhotonAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!