Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:21:17 

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   ParameterSet conf_;
0060 
0061   unsigned int ev;
0062   // ----------member data ---------------------------
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 // constants, enums and typedefs
0084 //
0085 
0086 //
0087 // static data member definitions
0088 //
0089 
0090 //
0091 // constructors and destructor
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   // Efficiency plots
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   // Ecal map
0135   etaphi_nonreco = fs->make<TH2F>("etaphi_nonreco", " ", 50, -2.5, 2.5, 50, -3.2, 3.2);
0136 
0137   // Energies
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   // do anything here that needs to be done at desctruction time
0159   // (e.g. close files, deallocate resources etc.)
0160 }
0161 
0162 //
0163 // member functions
0164 //
0165 
0166 // ------------ method called to for each event  ------------
0167 void EgGEDPhotonAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0168   // Candidate info
0169   Handle<reco::PFCandidateCollection> collection;
0170   InputTag label("particleFlow");  // <- Special electron coll.
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   // Validation from generator events
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             // all for the moment
0263           }
0264         }  // End PFCandidates Electron Selection
0265       }    // End Loop PFCandidates
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       //Add loop on the GED electrons
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       }  //End Loop Generator Particles
0376 
0377     }  //End IF Generator Particles
0378 
0379   }  //End Loop Generator Particles
0380 }
0381 // ------------ method called once each job just before starting event loop  ------------
0382 void EgGEDPhotonAnalyzer::beginJob() { ev = 0; }
0383 
0384 // ------------ method called once each job just after ending the event loop  ------------
0385 void EgGEDPhotonAnalyzer::endJob() { cout << " endJob:: #events " << ev << endl; }
0386 //define this as a plug-in
0387 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!