Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-22 02:24:07

0001 // -*- C++ -*-
0002 //
0003 // Package:    GsfGEDElectronAnalyzer
0004 // Class:      GsfGEDElectronAnalyzer
0005 //
0006 /**\class GsfGEDElectronAnalyzer
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/GsfElectron.h"
0031 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0032 #include "DataFormats/METReco/interface/PFMET.h"
0033 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0034 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0035 #include "DataFormats/Math/interface/normalizedPhi.h"
0036 
0037 #include <vector>
0038 #include <TROOT.h>
0039 #include <TFile.h>
0040 #include <TTree.h>
0041 #include <TH1F.h>
0042 #include <TH2F.h>
0043 #include <TLorentzVector.h>
0044 //
0045 // class decleration
0046 //
0047 
0048 using namespace edm;
0049 using namespace reco;
0050 using namespace std;
0051 class GsfGEDElectronAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0052 public:
0053   explicit GsfGEDElectronAnalyzer(const edm::ParameterSet &);
0054   ~GsfGEDElectronAnalyzer() override;
0055 
0056 private:
0057   void beginJob() override;
0058   void analyze(const edm::Event &, const edm::EventSetup &) override;
0059   void endJob() override;
0060 
0061   ParameterSet conf_;
0062 
0063   unsigned int ev;
0064   // ----------member data ---------------------------
0065 
0066   TH1F *h_etamc_ele, *h_etaec_ele, *h_etaecTrDr_ele, *h_etaecEcDr_ele, *h_etaeg_ele, *h_etaegTrDr_ele, *h_etaegEcDr_ele,
0067       *h_etaged_ele, *h_etagedTrDr_ele, *h_etagedEcDr_ele;
0068 
0069   TH1F *h_ptmc_ele, *h_ptec_ele, *h_ptecEcDr_ele, *h_ptecTrDr_ele, *h_pteg_ele, *h_ptegTrDr_ele, *h_ptegEcDr_ele,
0070       *h_ptged_ele, *h_ptgedTrDr_ele, *h_ptgedEcDr_ele;
0071 
0072   TH1F *h_mva_ele;
0073   TH2F *etaphi_nonreco;
0074 
0075   TH1F *eg_EEcalEtrue_1, *eg_EEcalEtrue_2, *eg_EEcalEtrue_3, *eg_EEcalEtrue_4, *eg_EEcalEtrue_5;
0076   TH1F *pf_EEcalEtrue_1, *pf_EEcalEtrue_2, *pf_EEcalEtrue_3, *pf_EEcalEtrue_4, *pf_EEcalEtrue_5;
0077   TH1F *ged_EEcalEtrue_1, *ged_EEcalEtrue_2, *ged_EEcalEtrue_3, *ged_EEcalEtrue_4, *ged_EEcalEtrue_5;
0078 
0079   TH1F *eg_e1x5_all, *eg_e1x5_eb, *eg_e1x5_ee, *ged_e1x5_all, *ged_e1x5_eb, *ged_e1x5_ee, *pf_e1x5_all, *pf_e1x5_eb,
0080       *pf_e1x5_ee;
0081 
0082   TH1F *eg_sihih_all, *eg_sihih_eb, *eg_sihih_ee, *ged_sihih_all, *ged_sihih_eb, *ged_sihih_ee, *pf_sihih_all,
0083       *pf_sihih_eb, *pf_sihih_ee;
0084 
0085   TH1F *eg_r9_all, *eg_r9_eb, *eg_r9_ee, *ged_r9_all, *ged_r9_eb, *ged_r9_ee, *pf_r9_all, *pf_r9_eb, *pf_r9_ee;
0086 
0087   TH1F *eg_scshh_all, *eg_scshh_eb, *eg_scshh_ee, *ged_scshh_all, *ged_scshh_eb, *ged_scshh_ee, *pf_scshh_all,
0088       *pf_scshh_eb, *pf_scshh_ee;
0089 
0090   TH1F *eg_scsff_all, *eg_scsff_eb, *eg_scsff_ee, *ged_scsff_all, *ged_scsff_eb, *ged_scsff_ee, *pf_scsff_all,
0091       *pf_scsff_eb, *pf_scsff_ee;
0092 
0093   TH1F *pf_MET_reco, *pf_MET_rereco;
0094 
0095   const edm::EDGetTokenT<reco::PFCandidateCollection> pfToken_;
0096   const edm::EDGetTokenT<reco::GsfElectronCollection> gsfEleToken_;
0097   const edm::EDGetTokenT<reco::GenParticleCollection> mcToken_;
0098   const edm::EDGetTokenT<reco::GsfElectronCollection> gedEleToken_;
0099   const edm::EDGetTokenT<std::vector<reco::PFMET>> metToken_;
0100   const edm::EDGetTokenT<std::vector<reco::PFMET>> reMetToken_;
0101 };
0102 
0103 //
0104 // constants, enums and typedefs
0105 //
0106 
0107 //
0108 // static data member definitions
0109 //
0110 
0111 //
0112 // constructors and destructor
0113 //
0114 GsfGEDElectronAnalyzer::GsfGEDElectronAnalyzer(const edm::ParameterSet &iConfig)
0115     : pfToken_(consumes<reco::PFCandidateCollection>(edm::InputTag("particleFlow::reRECO"))),
0116       gsfEleToken_(consumes<reco::GsfElectronCollection>(edm::InputTag("gsfElectrons"))),
0117       mcToken_(consumes<reco::GenParticleCollection>(edm::InputTag("genParticles"))),
0118       gedEleToken_(consumes<reco::GsfElectronCollection>(edm::InputTag("gedGsfElectrons::reRECO"))),
0119       metToken_(consumes<std::vector<reco::PFMET>>(edm::InputTag("pfMet::RECO"))),
0120       reMetToken_(consumes<std::vector<reco::PFMET>>(edm::InputTag("pfMet::reRECO"))) {
0121   usesResource(TFileService::kSharedResource);
0122 
0123   edm::Service<TFileService> fs;
0124 
0125   // Efficiency plots
0126 
0127   h_etamc_ele = fs->make<TH1F>("h_etamc_ele", " ", 50, -2.5, 2.5);
0128   h_etaec_ele = fs->make<TH1F>("h_etaec_ele", " ", 50, -2.5, 2.5);
0129   h_etaecTrDr_ele = fs->make<TH1F>("h_etaecTrDr_ele", " ", 50, -2.5, 2.5);
0130   h_etaecEcDr_ele = fs->make<TH1F>("h_etaecEcDr_ele", " ", 50, -2.5, 2.5);
0131   h_etaeg_ele = fs->make<TH1F>("h_etaeg_ele", " ", 50, -2.5, 2.5);
0132   h_etaegTrDr_ele = fs->make<TH1F>("h_etaegTrDr_ele", " ", 50, -2.5, 2.5);
0133   h_etaegEcDr_ele = fs->make<TH1F>("h_etaegEcDr_ele", " ", 50, -2.5, 2.5);
0134   h_etaged_ele = fs->make<TH1F>("h_etaged_ele", " ", 50, -2.5, 2.5);
0135   h_etagedTrDr_ele = fs->make<TH1F>("h_etagedTrDr_ele", " ", 50, -2.5, 2.5);
0136   h_etagedEcDr_ele = fs->make<TH1F>("h_etagedEcDr_ele", " ", 50, -2.5, 2.5);
0137 
0138   h_ptmc_ele = fs->make<TH1F>("h_ptmc_ele", " ", 50, 0, 100.);
0139   h_ptec_ele = fs->make<TH1F>("h_ptec_ele", " ", 50, 0, 100.);
0140   h_ptecTrDr_ele = fs->make<TH1F>("h_ptecTrDr_ele", " ", 50, 0, 100.);
0141   h_ptecEcDr_ele = fs->make<TH1F>("h_ptecEcDr_ele", " ", 50, 0, 100.);
0142   h_pteg_ele = fs->make<TH1F>("h_pteg_ele", " ", 50, 0, 100.);
0143   h_ptegTrDr_ele = fs->make<TH1F>("h_ptegTrDr_ele", " ", 50, 0, 100.);
0144   h_ptegEcDr_ele = fs->make<TH1F>("h_ptegEcDr_ele", " ", 50, 0, 100.);
0145   h_ptged_ele = fs->make<TH1F>("h_ptged_ele", " ", 50, 0, 100.);
0146   h_ptgedTrDr_ele = fs->make<TH1F>("h_ptgedTrDr_ele", " ", 50, 0, 100.);
0147   h_ptgedEcDr_ele = fs->make<TH1F>("h_ptgedEcDr_ele", " ", 50, 0, 100.);
0148 
0149   h_mva_ele = fs->make<TH1F>("h_mva_ele", " ", 50, -1.1, 1.1);
0150 
0151   // Ecal map
0152   etaphi_nonreco = fs->make<TH2F>("etaphi_nonreco", " ", 50, -2.5, 2.5, 50, -3.2, 3.2);
0153 
0154   // Energies
0155   eg_EEcalEtrue_1 = fs->make<TH1F>("eg_EEcalEtrue_1", "  ", 50, 0., 2.);
0156   eg_EEcalEtrue_2 = fs->make<TH1F>("eg_EEcalEtrue_2", "  ", 50, 0., 2.);
0157   eg_EEcalEtrue_3 = fs->make<TH1F>("eg_EEcalEtrue_3", "  ", 50, 0., 2.);
0158   eg_EEcalEtrue_4 = fs->make<TH1F>("eg_EEcalEtrue_4", "  ", 50, 0., 2.);
0159   eg_EEcalEtrue_5 = fs->make<TH1F>("eg_EEcalEtrue_5", "  ", 50, 0., 2.);
0160 
0161   pf_EEcalEtrue_1 = fs->make<TH1F>("pf_EEcalEtrue_1", "  ", 50, 0., 2.);
0162   pf_EEcalEtrue_2 = fs->make<TH1F>("pf_EEcalEtrue_2", "  ", 50, 0., 2.);
0163   pf_EEcalEtrue_3 = fs->make<TH1F>("pf_EEcalEtrue_3", "  ", 50, 0., 2.);
0164   pf_EEcalEtrue_4 = fs->make<TH1F>("pf_EEcalEtrue_4", "  ", 50, 0., 2.);
0165   pf_EEcalEtrue_5 = fs->make<TH1F>("pf_EEcalEtrue_5", "  ", 50, 0., 2.);
0166 
0167   ged_EEcalEtrue_1 = fs->make<TH1F>("ged_EEcalEtrue_1", "  ", 50, 0., 2.);
0168   ged_EEcalEtrue_2 = fs->make<TH1F>("ged_EEcalEtrue_2", "  ", 50, 0., 2.);
0169   ged_EEcalEtrue_3 = fs->make<TH1F>("ged_EEcalEtrue_3", "  ", 50, 0., 2.);
0170   ged_EEcalEtrue_4 = fs->make<TH1F>("ged_EEcalEtrue_4", "  ", 50, 0., 2.);
0171   ged_EEcalEtrue_5 = fs->make<TH1F>("ged_EEcalEtrue_5", "  ", 50, 0., 2.);
0172 
0173   //shower shapes
0174   eg_e1x5_all = fs->make<TH1F>("eg_e1x5_all", "  ", 60, 0., 300.);
0175   eg_e1x5_eb = fs->make<TH1F>("eg_e1x5_eb", "  ", 60, 0., 300.);
0176   eg_e1x5_ee = fs->make<TH1F>("eg_e1x5_ee", "  ", 60, 0., 300.);
0177 
0178   ged_e1x5_all = fs->make<TH1F>("ged_e1x5_all", "  ", 60, 0., 300.);
0179   ged_e1x5_eb = fs->make<TH1F>("ged_e1x5_eb", "  ", 60, 0., 300.);
0180   ged_e1x5_ee = fs->make<TH1F>("ged_e1x5_ee", "  ", 60, 0., 300.);
0181 
0182   pf_e1x5_all = fs->make<TH1F>("pf_e1x5_all", "  ", 60, 0., 300.);
0183   pf_e1x5_eb = fs->make<TH1F>("pf_e1x5_eb", "  ", 60, 0., 300.);
0184   pf_e1x5_ee = fs->make<TH1F>("pf_e1x5_ee", "  ", 60, 0., 300.);
0185 
0186   eg_sihih_all = fs->make<TH1F>("eg_sihih_all", "  ", 100, 0.0, 0.05);
0187   eg_sihih_eb = fs->make<TH1F>("eg_sihih_eb", "  ", 100, 0.0, 0.05);
0188   eg_sihih_ee = fs->make<TH1F>("eg_sihih_ee", "  ", 100, 0.0, 0.05);
0189 
0190   ged_sihih_all = fs->make<TH1F>("ged_sihih_all", "  ", 100, 0.0, 0.05);
0191   ged_sihih_eb = fs->make<TH1F>("ged_sihih_eb", "  ", 100, 0.0, 0.05);
0192   ged_sihih_ee = fs->make<TH1F>("ged_sihih_ee", "  ", 100, 0.0, 0.05);
0193 
0194   pf_sihih_all = fs->make<TH1F>("pf_sihih_all", "  ", 100, 0.0, 0.05);
0195   pf_sihih_eb = fs->make<TH1F>("pf_sihih_eb", "  ", 100, 0.0, 0.05);
0196   pf_sihih_ee = fs->make<TH1F>("pf_sihih_ee", "  ", 100, 0.0, 0.05);
0197 
0198   eg_r9_all = fs->make<TH1F>("eg_r9_all", "  ", 200, 0.0, 2.0);
0199   eg_r9_eb = fs->make<TH1F>("eg_r9_eb", "  ", 200, 0.0, 2.0);
0200   eg_r9_ee = fs->make<TH1F>("eg_r9_ee", "  ", 200, 0.0, 2.0);
0201 
0202   ged_r9_all = fs->make<TH1F>("ged_r9_all", "  ", 200, 0.0, 2.0);
0203   ged_r9_eb = fs->make<TH1F>("ged_r9_eb", "  ", 200, 0.0, 2.0);
0204   ged_r9_ee = fs->make<TH1F>("ged_r9_ee", "  ", 200, 0.0, 2.0);
0205 
0206   pf_r9_all = fs->make<TH1F>("pf_r9_all", "  ", 200, 0.0, 2.0);
0207   pf_r9_eb = fs->make<TH1F>("pf_r9_eb", "  ", 200, 0.0, 2.0);
0208   pf_r9_ee = fs->make<TH1F>("pf_r9_ee", "  ", 200, 0.0, 2.0);
0209 
0210   eg_scshh_all = fs->make<TH1F>("eg_scshh_all", "  ", 100, 0.0, 0.06);
0211   eg_scshh_eb = fs->make<TH1F>("eg_scshh_eb", "  ", 100, 0.0, 0.06);
0212   eg_scshh_ee = fs->make<TH1F>("eg_scshh_ee", "  ", 100, 0.0, 0.06);
0213 
0214   ged_scshh_all = fs->make<TH1F>("ged_scshh_all", "  ", 100, 0.0, 0.06);
0215   ged_scshh_eb = fs->make<TH1F>("ged_scshh_eb", "  ", 100, 0.0, 0.06);
0216   ged_scshh_ee = fs->make<TH1F>("ged_scshh_ee", "  ", 100, 0.0, 0.06);
0217 
0218   pf_scshh_all = fs->make<TH1F>("pf_scshh_all", "  ", 100, 0.0, 0.06);
0219   pf_scshh_eb = fs->make<TH1F>("pf_scshh_eb", "  ", 100, 0.0, 0.06);
0220   pf_scshh_ee = fs->make<TH1F>("pf_scshh_ee", "  ", 100, 0.0, 0.06);
0221 
0222   eg_scsff_all = fs->make<TH1F>("eg_scsff_all", "  ", 150, 0.0, 0.15);
0223   eg_scsff_eb = fs->make<TH1F>("eg_scsff_eb", "  ", 150, 0.0, 0.15);
0224   eg_scsff_ee = fs->make<TH1F>("eg_scsff_ee", "  ", 150, 0.0, 0.15);
0225 
0226   ged_scsff_all = fs->make<TH1F>("ged_scsff_all", "  ", 150, 0.0, 0.15);
0227   ged_scsff_eb = fs->make<TH1F>("ged_scsff_eb", "  ", 150, 0.0, 0.15);
0228   ged_scsff_ee = fs->make<TH1F>("ged_scsff_ee", "  ", 150, 0.0, 0.15);
0229 
0230   pf_scsff_all = fs->make<TH1F>("pf_scsff_all", "  ", 150, 0.0, 0.15);
0231   pf_scsff_eb = fs->make<TH1F>("pf_scsff_eb", "  ", 150, 0.0, 0.15);
0232   pf_scsff_ee = fs->make<TH1F>("pf_scsff_ee", "  ", 150, 0.0, 0.15);
0233 
0234   // MET
0235   pf_MET_reco = fs->make<TH1F>("pf_MET_reco", "  ", 10, 0, 1000);
0236   pf_MET_rereco = fs->make<TH1F>("pf_MET_rereco", "  ", 10, 0, 1000);
0237 }
0238 
0239 GsfGEDElectronAnalyzer::~GsfGEDElectronAnalyzer() {
0240   // do anything here that needs to be done at desctruction time
0241   // (e.g. close files, deallocate resources etc.)
0242 }
0243 
0244 //
0245 // member functions
0246 //
0247 
0248 // ------------ method called to for each event  ------------
0249 void GsfGEDElectronAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0250   // Candidate info
0251   const auto &candidates = iEvent.get(pfToken_);
0252   const auto &theGsfEle = iEvent.get(gsfEleToken_);
0253   const auto &pMCTruth = iEvent.get(mcToken_);
0254   const auto &theGedEle = iEvent.get(gedEleToken_);
0255   const auto &recoMet = iEvent.get(metToken_);
0256   const auto &rerecoMet = iEvent.get(reMetToken_);
0257 
0258   pf_MET_reco->Fill(recoMet[0].et());
0259   pf_MET_rereco->Fill(rerecoMet[0].et());
0260 
0261   bool debug = true;
0262 
0263   ev++;
0264 
0265   if (debug)
0266     cout << "************************* New Event:: " << ev << " *************************" << endl;
0267 
0268   // Validation from generator events
0269 
0270   for (const auto &cP : pMCTruth) {
0271     float etamc = cP.eta();
0272     float phimc = cP.phi();
0273     float ptmc = cP.pt();
0274     float Emc = cP.energy();
0275 
0276     if (abs(cP.pdgId()) == 11 && cP.status() == 1 && cP.pt() > 2. && fabs(cP.eta()) < 2.5) {
0277       h_etamc_ele->Fill(etamc);
0278       h_ptmc_ele->Fill(ptmc);
0279 
0280       float MindR = 1000.;
0281       float mvaCutEle = -1.;
0282       bool isPfTrDr = false;
0283       bool isPfEcDr = false;
0284 
0285       if (debug)
0286         cout << " MC particle:  pt " << ptmc << " eta,phi " << etamc << ", " << phimc << endl;
0287 
0288       for (const auto &cand : candidates) {
0289         reco::PFCandidate::ParticleType type = cand.particleId();
0290 
0291         if (type == reco::PFCandidate::e) {
0292           float eta = cand.eta();
0293           float phi = cand.phi();
0294           float pfmva = cand.mva_e_pi();
0295 
0296           reco::GsfTrackRef refGsf = cand.gsfTrackRef();
0297           //ElectronSeedRef seedRef= refGsf->extra()->seedRef().castTo<ElectronSeedRef>();
0298 
0299           float deta = etamc - eta;
0300           float dphi = normalizedPhi(phimc - phi);
0301           float dR = sqrt(deta * deta + dphi * dphi);
0302 
0303           if (dR < 0.05) {
0304             MindR = dR;
0305             mvaCutEle = cand.mva_e_pi();
0306             /*
0307         if(seedRef->isEcalDriven())
0308           isPfEcDr = true;
0309         if(seedRef->isTrackerDriven())
0310           isPfTrDr = true;
0311         */
0312 
0313             if (debug)
0314               cout << " PF ele matched:  pt " << cand.pt() << " (" << cand.ecalEnergy() / std::cosh(eta) << ") "
0315                    << " eta,phi " << eta << ", " << phi << " pfmva " << pfmva << endl;
0316             // all for the moment
0317           }
0318         }  // End PFCandidates Electron Selection
0319       }    // End Loop PFCandidates
0320 
0321       if (MindR < 0.05) {
0322         h_mva_ele->Fill(mvaCutEle);
0323         h_etaec_ele->Fill(etamc);
0324         h_ptec_ele->Fill(ptmc);
0325 
0326         if (isPfEcDr) {
0327           h_etaecEcDr_ele->Fill(etamc);
0328           h_ptecEcDr_ele->Fill(ptmc);
0329         }
0330         if (isPfTrDr) {
0331           h_etaecTrDr_ele->Fill(etamc);
0332           h_ptecTrDr_ele->Fill(ptmc);
0333         }
0334       } else {
0335         etaphi_nonreco->Fill(etamc, phimc);
0336       }
0337 
0338       float MindREG = 1000;
0339       bool isEgTrDr = false;
0340       bool isEgEcDr = false;
0341 
0342       for (const auto &gsfEle : theGsfEle) {
0343         reco::GsfTrackRef egGsfTrackRef = gsfEle.gsfTrack();
0344         float etareco = gsfEle.eta();
0345         float phireco = gsfEle.phi();
0346 
0347         float pfmva = gsfEle.mva_e_pi();
0348 
0349         reco::GsfTrackRef refGsf = gsfEle.gsfTrack();
0350         //ElectronSeedRef seedRef= refGsf->extra()->seedRef().castTo<ElectronSeedRef>();
0351 
0352         float deta = etamc - etareco;
0353         float dphi = normalizedPhi(phimc - phireco);
0354         float dR = sqrt(deta * deta + dphi * dphi);
0355 
0356         float SCEnergy = gsfEle.superCluster()->energy();
0357         float ErecoEtrue = SCEnergy / Emc;
0358 
0359         if (dR < 0.05) {
0360           if (debug)
0361             cout << " EG ele matched: pt " << gsfEle.pt() << " (" << SCEnergy / std::cosh(etareco) << ") "
0362                  << " eta,phi " << etareco << ", " << phireco << " pfmva " << pfmva << endl;
0363 
0364           if (fabs(etamc) < 0.5) {
0365             eg_EEcalEtrue_1->Fill(ErecoEtrue);
0366           }
0367           if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0368             eg_EEcalEtrue_2->Fill(ErecoEtrue);
0369           }
0370           if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0371             eg_EEcalEtrue_3->Fill(ErecoEtrue);
0372           }
0373           if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0374             eg_EEcalEtrue_4->Fill(ErecoEtrue);
0375           }
0376           if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0377             eg_EEcalEtrue_5->Fill(ErecoEtrue);
0378           }
0379 
0380           if ((gsfEle.parentSuperCluster()).isNonnull()) {
0381             float SCPF = gsfEle.parentSuperCluster()->rawEnergy();
0382             float EpfEtrue = SCPF / Emc;
0383             if (fabs(etamc) < 0.5) {
0384               pf_EEcalEtrue_1->Fill(EpfEtrue);
0385             }
0386             if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0387               pf_EEcalEtrue_2->Fill(EpfEtrue);
0388             }
0389             if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0390               pf_EEcalEtrue_3->Fill(EpfEtrue);
0391             }
0392             if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0393               pf_EEcalEtrue_4->Fill(EpfEtrue);
0394             }
0395             if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0396               pf_EEcalEtrue_5->Fill(EpfEtrue);
0397             }
0398             const reco::GsfElectron::ShowerShape &pfshapes = gsfEle.showerShape();
0399 
0400             pf_e1x5_all->Fill(pfshapes.e1x5);
0401             pf_sihih_all->Fill(pfshapes.sigmaIetaIeta);
0402             pf_r9_all->Fill(pfshapes.r9);
0403             pf_scshh_all->Fill(gsfEle.parentSuperCluster()->etaWidth());
0404             pf_scsff_all->Fill(gsfEle.parentSuperCluster()->phiWidth());
0405             if (std::abs(etareco) < 1.479) {
0406               pf_e1x5_eb->Fill(pfshapes.e1x5);
0407               pf_sihih_eb->Fill(pfshapes.sigmaIetaIeta);
0408               pf_r9_eb->Fill(pfshapes.r9);
0409               pf_scshh_eb->Fill(gsfEle.parentSuperCluster()->etaWidth());
0410               pf_scsff_eb->Fill(gsfEle.parentSuperCluster()->phiWidth());
0411             }
0412             if (std::abs(etareco) >= 1.479) {
0413               pf_e1x5_ee->Fill(pfshapes.e1x5);
0414               pf_sihih_ee->Fill(pfshapes.sigmaIetaIeta);
0415               pf_r9_ee->Fill(pfshapes.r9);
0416               pf_scshh_ee->Fill(gsfEle.parentSuperCluster()->etaWidth());
0417               pf_scsff_ee->Fill(gsfEle.parentSuperCluster()->phiWidth());
0418             }
0419           }
0420 
0421           eg_e1x5_all->Fill(gsfEle.e1x5());
0422           eg_sihih_all->Fill(gsfEle.sigmaIetaIeta());
0423           eg_r9_all->Fill(gsfEle.r9());
0424           eg_scshh_all->Fill(gsfEle.superCluster()->etaWidth());
0425           eg_scsff_all->Fill(gsfEle.superCluster()->phiWidth());
0426           if (std::abs(etareco) < 1.479) {
0427             eg_e1x5_eb->Fill(gsfEle.e1x5());
0428             eg_sihih_eb->Fill(gsfEle.sigmaIetaIeta());
0429             eg_r9_eb->Fill(gsfEle.r9());
0430             eg_scshh_eb->Fill(gsfEle.superCluster()->etaWidth());
0431             eg_scsff_eb->Fill(gsfEle.superCluster()->phiWidth());
0432           }
0433           if (std::abs(etareco) >= 1.479) {
0434             eg_e1x5_ee->Fill(gsfEle.e1x5());
0435             eg_sihih_ee->Fill(gsfEle.sigmaIetaIeta());
0436             eg_r9_ee->Fill(gsfEle.r9());
0437             eg_scshh_ee->Fill(gsfEle.superCluster()->etaWidth());
0438             eg_scsff_ee->Fill(gsfEle.superCluster()->phiWidth());
0439           }
0440 
0441           MindREG = dR;
0442           /*
0443       if(seedRef->isEcalDriven())
0444         isEgEcDr = true;
0445       if(seedRef->isTrackerDriven())
0446         isEgTrDr = true;
0447       */
0448         }
0449       }
0450       if (MindREG < 0.05) {
0451         h_etaeg_ele->Fill(etamc);
0452         h_pteg_ele->Fill(ptmc);
0453 
0454         // This is to understand the contribution of each seed type.
0455         if (isEgEcDr) {
0456           h_etaegEcDr_ele->Fill(etamc);
0457           h_ptegEcDr_ele->Fill(ptmc);
0458         }
0459         if (isEgTrDr) {
0460           h_etaegTrDr_ele->Fill(etamc);
0461           h_ptegTrDr_ele->Fill(ptmc);
0462         }
0463       }
0464 
0465       //Add loop on the GED electrons
0466 
0467       float MindRGedEg = 1000;
0468       bool isGedEgTrDr = false;
0469       bool isGedEgEcDr = false;
0470 
0471       for (const auto &gedEle : theGedEle) {
0472         reco::GsfTrackRef egGsfTrackRef = gedEle.gsfTrack();
0473         float etareco = gedEle.eta();
0474         float phireco = gedEle.phi();
0475         float pfmva = gedEle.mva_e_pi();
0476 
0477         reco::GsfTrackRef refGsf = gedEle.gsfTrack();
0478         //ElectronSeedRef seedRef= refGsf->extra()->seedRef().castTo<ElectronSeedRef>();
0479 
0480         float deta = etamc - etareco;
0481         float dphi = normalizedPhi(phimc - phireco);
0482         float dR = sqrt(deta * deta + dphi * dphi);
0483 
0484         float SCEnergy = gedEle.superCluster()->energy();
0485         float ErecoEtrue = SCEnergy / Emc;
0486 
0487         if (dR < 0.05) {
0488           const reco::PFCandidate *matchPF = NULL;
0489           if (debug)
0490             cout << " GED ele matched: pt " << gedEle.pt() << " (" << SCEnergy / std::cosh(etareco) << ") "
0491                  << " eta,phi " << etareco << ", " << phireco << " pfmva " << pfmva << endl;
0492 
0493           for (unsigned k = 0; k < candidates.size(); ++k) {
0494             if (std::abs(candidates[k].pdgId()) == 11) {
0495               reco::GsfTrackRef gsfEleGsfTrackRef = candidates[k].gsfTrackRef();
0496               if (gsfEleGsfTrackRef->ptMode() == egGsfTrackRef->ptMode()) {
0497                 matchPF = &candidates[k];
0498               }
0499             }
0500           }
0501 
0502           if (debug && matchPF)
0503             std::cout << "GED-PF match!" << std::endl;
0504 
0505           if (ErecoEtrue < 0.6)
0506             std::cout << "++bad Ereco/Etrue: " << ErecoEtrue << std::endl;
0507 
0508           if (fabs(etamc) < 0.5) {
0509             ged_EEcalEtrue_1->Fill(ErecoEtrue);
0510           }
0511           if (fabs(etamc) >= 0.5 && fabs(etamc) < 1.0) {
0512             ged_EEcalEtrue_2->Fill(ErecoEtrue);
0513           }
0514           if (fabs(etamc) >= 1.0 && fabs(etamc) < 1.5) {
0515             ged_EEcalEtrue_3->Fill(ErecoEtrue);
0516           }
0517           if (fabs(etamc) >= 1.5 && fabs(etamc) < 2.0) {
0518             ged_EEcalEtrue_4->Fill(ErecoEtrue);
0519           }
0520           if (fabs(etamc) >= 2.0 && fabs(etamc) < 2.5) {
0521             ged_EEcalEtrue_5->Fill(ErecoEtrue);
0522             /*
0523         if( debug && matchPF && ErecoEtrue > 1.30 ) {
0524           cout << " -PF ele matched: pt " 
0525            << matchPF->pt()
0526          << " eta,phi " <<  matchPF->eta() 
0527          << ", " << matchPF->phi() << " pfmva " 
0528          <<  matchPF->mva_e_pi() << endl;
0529           std::cout << "GED Clusters: " << std::endl;
0530           for( const auto& c : gedEle.superCluster()->clusters() ) {
0531         std::cout << " E/eta/phi : " << c->energy() 
0532               << '/' << c->eta() 
0533               << '/' << c->phi() << std::endl;      
0534           }
0535           std::cout << "-PF Clusters: " << std::endl;
0536           for( const auto& c : matchPF->superClusterRef()->clusters() ) {
0537         std::cout << " E/eta/phi : " << c->energy() 
0538               << '/' << c->eta() 
0539               << '/' << c->phi() << std::endl;
0540           }
0541           
0542         }
0543         */
0544           }
0545 
0546           ged_e1x5_all->Fill(gedEle.e1x5());
0547           ged_sihih_all->Fill(gedEle.sigmaIetaIeta());
0548           ged_r9_all->Fill(gedEle.r9());
0549           ged_scshh_all->Fill(gedEle.superCluster()->etaWidth());
0550           ged_scsff_all->Fill(gedEle.superCluster()->phiWidth());
0551           if (std::abs(etareco) < 1.479) {
0552             ged_e1x5_eb->Fill(gedEle.e1x5());
0553             ged_sihih_eb->Fill(gedEle.sigmaIetaIeta());
0554             ged_r9_eb->Fill(gedEle.r9());
0555             ged_scshh_eb->Fill(gedEle.superCluster()->etaWidth());
0556             ged_scsff_eb->Fill(gedEle.superCluster()->phiWidth());
0557           }
0558           if (std::abs(etareco) >= 1.479) {
0559             ged_e1x5_ee->Fill(gedEle.e1x5());
0560             ged_sihih_ee->Fill(gedEle.sigmaIetaIeta());
0561             ged_r9_ee->Fill(gedEle.r9());
0562             ged_scshh_ee->Fill(gedEle.superCluster()->etaWidth());
0563             ged_scsff_eb->Fill(gedEle.superCluster()->phiWidth());
0564           }
0565 
0566           MindRGedEg = dR;
0567           /*
0568       if(seedRef->isEcalDriven())
0569         isGedEgEcDr = true;
0570       if(seedRef->isTrackerDriven())
0571         isGedEgTrDr = true;
0572       */
0573         }
0574       }
0575       if (MindRGedEg < 0.05) {
0576         h_etaged_ele->Fill(etamc);
0577         h_ptged_ele->Fill(ptmc);
0578 
0579         // This is to understand the contribution of each seed type.
0580         if (isGedEgEcDr) {
0581           h_etagedEcDr_ele->Fill(etamc);
0582           h_ptgedEcDr_ele->Fill(ptmc);
0583         }
0584         if (isGedEgTrDr) {
0585           h_etagedTrDr_ele->Fill(etamc);
0586           h_ptgedTrDr_ele->Fill(ptmc);
0587         }
0588       }  //End Loop Generator Particles
0589 
0590     }  //End IF Generator Particles
0591 
0592   }  //End Loop Generator Particles
0593 }
0594 // ------------ method called once each job just before starting event loop  ------------
0595 void GsfGEDElectronAnalyzer::beginJob() { ev = 0; }
0596 
0597 // ------------ method called once each job just after ending the event loop  ------------
0598 void GsfGEDElectronAnalyzer::endJob() { cout << " endJob:: #events " << ev << endl; }
0599 //define this as a plug-in
0600 DEFINE_FWK_MODULE(GsfGEDElectronAnalyzer);