Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-04-01 02:23:10 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ElectronTestAnalyzer
0004 // Class:      ElectronTestAnalyzer
0005 //
0006 /**\class ElectronTestAnalyzer
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/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/EDAnalyzer.h"
0022 
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 
0026 #include "DataFormats/MuonReco/interface/Muon.h"
0027 
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0031 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0032 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0033 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0034 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0035 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0036 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0037 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0038 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0039 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0040 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0041 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0042 #include "EgammaAnalysis/ElectronTools/interface/EGammaMvaEleEstimator.h"
0043 #include "DataFormats/VertexReco/interface/Vertex.h"
0044 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0045 #include "TrackingTools/IPTools/interface/IPTools.h"
0046 #include "FWCore/Utilities/interface/isFinite.h"
0047 #include "DataFormats/Math/interface/normalizedPhi.h"
0048 
0049 #include <cmath>
0050 #include <vector>
0051 #include <TROOT.h>
0052 #include <TFile.h>
0053 #include <TTree.h>
0054 #include <TH1F.h>
0055 #include <TH2F.h>
0056 #include <TLorentzVector.h>
0057 #include "TMVA/Factory.h"
0058 #include "TMVA/Tools.h"
0059 #include "TMVA/Reader.h"
0060 
0061 //
0062 // class decleration
0063 //
0064 
0065 class ElectronTestAnalyzer : public edm::EDAnalyzer {
0066 public:
0067   explicit ElectronTestAnalyzer(const edm::ParameterSet&);
0068   ~ElectronTestAnalyzer();
0069 
0070 private:
0071   virtual void beginJob() override;
0072   virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
0073   virtual void endJob() override;
0074   virtual void myBindVariables();
0075   virtual void myVar(const reco::GsfElectron& ele,
0076                      const reco::Vertex& vertex,
0077                      const TransientTrackBuilder& transientTrackBuilder,
0078                      EcalClusterLazyTools const& myEcalCluster,
0079                      bool printDebug = kFALSE);
0080   virtual void evaluate_mvas(const edm::Event& iEvent, const edm::EventSetup& iSetup);
0081 
0082   bool trainTrigPresel(const reco::GsfElectron& ele);
0083 
0084   edm::ParameterSet conf_;
0085 
0086   edm::EDGetTokenT<reco::GsfElectronCollection> gsfEleToken_;
0087   edm::EDGetTokenT<reco::GenParticleCollection> genToken_;
0088   //edm::EDGetTokenT<edm::HepMCProduct>  mcTruthToken_;
0089   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0090   //edm::EDGetTokenT<reco::PFCandidateCollection> pfCandToken_;
0091   edm::EDGetTokenT<double> eventrhoToken_;
0092   edm::EDGetTokenT<reco::MuonCollection> muonToken_;
0093   edm::EDGetTokenT<EcalRecHitCollection> reducedEBRecHitCollectionToken_;
0094   edm::EDGetTokenT<EcalRecHitCollection> reducedEERecHitCollectionToken_;
0095   const EcalClusterLazyTools::ESGetTokens ecalClusterESGetTokens_;
0096 
0097   EGammaMvaEleEstimator* myMVATrigV0;
0098   EGammaMvaEleEstimator* myMVATrigNoIPV0;
0099   EGammaMvaEleEstimator* myMVANonTrigV0;
0100 
0101   TMVA::Reader* myTMVAReader;
0102   Float_t myMVAVar_fbrem;
0103   Float_t myMVAVar_kfchi2;
0104   Float_t myMVAVar_kfhits;
0105   Float_t myMVAVar_gsfchi2;
0106 
0107   Float_t myMVAVar_deta;
0108   Float_t myMVAVar_dphi;
0109   Float_t myMVAVar_detacalo;
0110   Float_t myMVAVar_dphicalo;
0111 
0112   Float_t myMVAVar_see;
0113   Float_t myMVAVar_spp;
0114   Float_t myMVAVar_etawidth;
0115   Float_t myMVAVar_phiwidth;
0116   Float_t myMVAVar_e1x5e5x5;
0117   Float_t myMVAVar_R9;
0118   Float_t myMVAVar_nbrems;
0119 
0120   Float_t myMVAVar_HoE;
0121   Float_t myMVAVar_EoP;
0122   Float_t myMVAVar_IoEmIoP;
0123   Float_t myMVAVar_eleEoPout;
0124   Float_t myMVAVar_PreShowerOverRaw;
0125   Float_t myMVAVar_EoPout;
0126 
0127   Float_t myMVAVar_d0;
0128   Float_t myMVAVar_ip3d;
0129 
0130   Float_t myMVAVar_eta;
0131   Float_t myMVAVar_pt;
0132 
0133   double _Rho;
0134   unsigned int ev;
0135   // ----------member data ---------------------------
0136 
0137   TH1F *h_mva_nonTrig, *h_mva_trig, *h_mva_trigNonIp;
0138   TH1F* h_fbrem;
0139   TH1F* h_kfchi2;
0140   TH1F* h_kfhits;
0141   TH1F* h_gsfchi2;
0142   TH1F* h_deta;
0143   TH1F* h_dphi;
0144   TH1F* h_detacalo;
0145   TH1F* h_dphicalo;
0146   TH1F* h_see;
0147   TH1F* h_spp;
0148   TH1F* h_etawidth;
0149   TH1F* h_phiwidth;
0150   TH1F* h_e1x5e5x5;
0151   TH1F* h_nbrems;
0152   TH1F* h_R9;
0153   TH1F* h_HoE;
0154   TH1F* h_EoP;
0155   TH1F* h_IoEmIoP;
0156   TH1F* h_eleEoPout;
0157   TH1F* h_EoPout;
0158   TH1F* h_PreShowerOverRaw;
0159   TH1F* h_pt;
0160 };
0161 
0162 //
0163 // constants, enums and typedefs
0164 //
0165 
0166 //
0167 // static data member definitions
0168 //
0169 
0170 //
0171 // constructors and destructor
0172 //
0173 ElectronTestAnalyzer::ElectronTestAnalyzer(const edm::ParameterSet& iConfig)
0174     : conf_(iConfig),
0175       gsfEleToken_(consumes<reco::GsfElectronCollection>(edm::InputTag("gsfElectrons"))),
0176       genToken_(consumes<reco::GenParticleCollection>(edm::InputTag("genParticles"))),
0177       //mcTruthToken_(consumes<edm::HepMCProduct>(edm::InputTag("generator"))),
0178       vertexToken_(consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"))),
0179       //pfCandToken_(consumes<reco::PFCandidateCollection>(edm::InputTag("particleFlow"))),
0180       eventrhoToken_(consumes<double>(edm::InputTag("kt6PFJets", "rho"))),
0181       muonToken_(consumes<reco::MuonCollection>(edm::InputTag("muons"))),
0182       reducedEBRecHitCollectionToken_(consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEB"))),
0183       reducedEERecHitCollectionToken_(consumes<EcalRecHitCollection>(edm::InputTag("reducedEcalRecHitsEE"))),
0184       ecalClusterESGetTokens_{consumesCollector()} {
0185   Bool_t manualCat = true;
0186 
0187   ev = 0;
0188 
0189   // NOTE: it is better if you copy the MVA weight files locally
0190 
0191   std::vector<std::string> catWTrigV0;
0192   catWTrigV0.push_back("../data/Electrons_BDTG_TrigV0_Cat1.weights.xml");
0193   catWTrigV0.push_back("../data/Electrons_BDTG_TrigV0_Cat2.weights.xml");
0194   catWTrigV0.push_back("../data/Electrons_BDTG_TrigV0_Cat3.weights.xml");
0195   catWTrigV0.push_back("../data/Electrons_BDTG_TrigV0_Cat4.weights.xml");
0196   catWTrigV0.push_back("../data/Electrons_BDTG_TrigV0_Cat5.weights.xml");
0197   catWTrigV0.push_back("../data/Electrons_BDTG_TrigV0_Cat6.weights.xml");
0198 
0199   myMVATrigV0 = new EGammaMvaEleEstimator();
0200   myMVATrigV0->initialize("BDT", EGammaMvaEleEstimator::kTrig, manualCat, catWTrigV0);
0201 
0202   std::vector<std::string> catWTrigNoIPV0;
0203   catWTrigNoIPV0.push_back("../data/Electrons_BDTG_TrigNoIPV0_2012_Cat1.weights.xml");
0204   catWTrigNoIPV0.push_back("../data/Electrons_BDTG_TrigNoIPV0_2012_Cat2.weights.xml");
0205   catWTrigNoIPV0.push_back("../data/Electrons_BDTG_TrigNoIPV0_2012_Cat3.weights.xml");
0206   catWTrigNoIPV0.push_back("../data/Electrons_BDTG_TrigNoIPV0_2012_Cat4.weights.xml");
0207   catWTrigNoIPV0.push_back("../data/Electrons_BDTG_TrigNoIPV0_2012_Cat5.weights.xml");
0208   catWTrigNoIPV0.push_back("../data/Electrons_BDTG_TrigNoIPV0_2012_Cat6.weights.xml");
0209 
0210   myMVATrigNoIPV0 = new EGammaMvaEleEstimator();
0211   myMVATrigNoIPV0->initialize("BDT", EGammaMvaEleEstimator::kTrigNoIP, manualCat, catWTrigNoIPV0);
0212 
0213   std::vector<std::string> catWNonTrigV0;
0214   catWNonTrigV0.push_back("../data/Electrons_BDTG_NonTrigV0_Cat1.weights.xml");
0215   catWNonTrigV0.push_back("../data/Electrons_BDTG_NonTrigV0_Cat2.weights.xml");
0216   catWNonTrigV0.push_back("../data/Electrons_BDTG_NonTrigV0_Cat3.weights.xml");
0217   catWNonTrigV0.push_back("../data/Electrons_BDTG_NonTrigV0_Cat4.weights.xml");
0218   catWNonTrigV0.push_back("../data/Electrons_BDTG_NonTrigV0_Cat5.weights.xml");
0219   catWNonTrigV0.push_back("../data/Electrons_BDTG_NonTrigV0_Cat6.weights.xml");
0220 
0221   myMVANonTrigV0 = new EGammaMvaEleEstimator();
0222   myMVANonTrigV0->initialize("BDT", EGammaMvaEleEstimator::kNonTrig, manualCat, catWNonTrigV0);
0223 
0224   edm::Service<TFileService> fs;
0225 
0226   h_mva_nonTrig = fs->make<TH1F>("h_mva_nonTrig", " ", 100, -1.1, 1.1);
0227   h_mva_trig = fs->make<TH1F>("h_mva_trig", " ", 100, -1.1, 1.1);
0228   h_mva_trigNonIp = fs->make<TH1F>("h_mva_trigNonIp", " ", 100, -1.1, 1.1);
0229 
0230   h_fbrem = fs->make<TH1F>("h_fbrem", " ", 100, -1., 1.);
0231   h_kfchi2 = fs->make<TH1F>("h_kfchi2", " ", 100, 0, 15);
0232   h_kfhits = fs->make<TH1F>("h_kfhits", " ", 25, 0, 25);
0233   h_gsfchi2 = fs->make<TH1F>("h_gsfchi2", " ", 100, 0., 50);
0234 
0235   h_deta = fs->make<TH1F>("h_deta", " ", 100, 0., 0.06);
0236   h_dphi = fs->make<TH1F>("h_dphi", " ", 100, 0., 0.3);
0237   h_detacalo = fs->make<TH1F>("h_detacalo", " ", 100, 0., 0.05);
0238   h_dphicalo = fs->make<TH1F>("h_dphicalo", " ", 100, 0., 0.2);
0239   h_see = fs->make<TH1F>("h_see", " ", 100, 0., 0.06);
0240   h_spp = fs->make<TH1F>("h_spp", " ", 100, 0., 0.09);
0241   h_etawidth = fs->make<TH1F>("h_etawidth", " ", 100, 0., 0.1);
0242   h_phiwidth = fs->make<TH1F>("h_phiwidth", " ", 100, 0., 0.2);
0243   h_e1x5e5x5 = fs->make<TH1F>("h_e1x5e5x5", " ", 100, -0.1, 1.1);
0244   h_R9 = fs->make<TH1F>("h_R9", " ", 100, 0., 2.);
0245   h_nbrems = fs->make<TH1F>("h_nbrems", " ", 100, 0., 10);
0246   h_HoE = fs->make<TH1F>("h_HoE", " ", 100, 0., 0.5);
0247   h_EoP = fs->make<TH1F>("h_EoP", " ", 100, 0., 5.);
0248   h_IoEmIoP = fs->make<TH1F>("h_IoEmIoP", " ", 100, -0.15, 0.15);
0249   h_eleEoPout = fs->make<TH1F>("h_eleEoPout", " ", 100, 0., 5.);
0250   h_EoPout = fs->make<TH1F>("h_EoPout", " ", 100, 0., 5.);
0251   h_PreShowerOverRaw = fs->make<TH1F>("h_PreShowerOverRaw", " ", 100, 0., 0.3);
0252   h_pt = fs->make<TH1F>("h_pt", " ", 100, 0., 50);
0253 }
0254 
0255 ElectronTestAnalyzer::~ElectronTestAnalyzer() {
0256   // do anything here that needs to be done at desctruction time
0257   // (e.g. close files, deallocate resources etc.)
0258 }
0259 
0260 //
0261 // member functions
0262 //
0263 
0264 // ------------ method called to for each event  ------------
0265 void ElectronTestAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0266   ElectronTestAnalyzer::evaluate_mvas(iEvent, iSetup);
0267 
0268   edm::Handle<reco::GsfElectronCollection> theEGammaCollection;
0269   iEvent.getByToken(gsfEleToken_, theEGammaCollection);
0270   const reco::GsfElectronCollection theEGamma = *(theEGammaCollection.product());
0271 
0272   edm::Handle<reco::GenParticleCollection> genParticles;
0273   iEvent.getByToken(genToken_, genParticles);
0274   //InputTag  mcTruthToken_(string("generator"));
0275   //edm::Handle<edm::HepMCProduct> pMCTruth;
0276   //iEvent.getByToken(mcTruthToken_,pMCTruth);
0277   //const HepMC::GenEvent* genEvent = pMCTruth->GetEvent();
0278 
0279   edm::Handle<reco::VertexCollection> thePrimaryVertexColl;
0280   iEvent.getByToken(vertexToken_, thePrimaryVertexColl);
0281 
0282   _Rho = 0;
0283   edm::Handle<double> rhoPtr;
0284   iEvent.getByToken(eventrhoToken_, rhoPtr);
0285   _Rho = *rhoPtr;
0286 
0287   reco::Vertex dummy;
0288   const reco::Vertex* pv = &dummy;
0289   if (thePrimaryVertexColl->size() != 0) {
0290     pv = &*thePrimaryVertexColl->begin();
0291   } else {  // create a dummy PV
0292     reco::Vertex::Error e;
0293     e(0, 0) = 0.0015 * 0.0015;
0294     e(1, 1) = 0.0015 * 0.0015;
0295     e(2, 2) = 15. * 15.;
0296     reco::Vertex::Point p(0, 0, 0);
0297     dummy = reco::Vertex(p, e, 0, 0, 0);
0298   }
0299   EcalClusterLazyTools lazyTools(
0300       iEvent, ecalClusterESGetTokens_.get(iSetup), reducedEBRecHitCollectionToken_, reducedEERecHitCollectionToken_);
0301 
0302   edm::ESHandle<TransientTrackBuilder> builder;
0303   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
0304   TransientTrackBuilder thebuilder = *(builder.product());
0305 
0306   bool debug = true;
0307   bool debugMVAclass = false;
0308   bool debugMyVar = false;
0309 
0310   ev++;
0311 
0312   // Validation from generator events
0313 
0314   for (size_t i(0); i < genParticles->size(); i++) {
0315     const reco::GenParticle& genPtcl = (*genParticles)[i];
0316     float etamc = genPtcl.eta();
0317     float phimc = genPtcl.phi();
0318     float ptmc = genPtcl.pt();
0319 
0320     if (abs(genPtcl.pdgId()) == 11 && genPtcl.status() == 1 && ptmc > 10. && std::abs(etamc) < 2.5) {
0321       for (unsigned int j = 0; j < theEGamma.size(); j++) {
0322         float etareco = theEGamma[j].eta();
0323         float phireco = theEGamma[j].phi();
0324         float deta = etamc - etareco;
0325         float dphi = normalizedPhi(phimc - phireco);
0326         float dR = sqrt(deta * deta + dphi * dphi);
0327 
0328         if (dR < 0.1) {
0329           myVar((theEGamma[j]), *pv, thebuilder, lazyTools, debugMyVar);
0330           myBindVariables();
0331 
0332           h_fbrem->Fill(myMVAVar_fbrem);
0333           h_kfchi2->Fill(myMVAVar_kfchi2);
0334           h_kfhits->Fill(myMVAVar_kfhits);
0335           h_gsfchi2->Fill(myMVAVar_gsfchi2);
0336 
0337           h_deta->Fill(myMVAVar_deta);
0338           h_dphi->Fill(myMVAVar_dphi);
0339           h_detacalo->Fill(myMVAVar_detacalo);
0340           h_dphicalo->Fill(myMVAVar_dphicalo);
0341 
0342           h_see->Fill(myMVAVar_see);
0343           h_spp->Fill(myMVAVar_spp);
0344           h_etawidth->Fill(myMVAVar_etawidth);
0345           h_phiwidth->Fill(myMVAVar_phiwidth);
0346           h_e1x5e5x5->Fill(myMVAVar_e1x5e5x5);
0347           h_R9->Fill(myMVAVar_R9);
0348           h_nbrems->Fill(myMVAVar_nbrems);
0349 
0350           h_HoE->Fill(myMVAVar_HoE);
0351           h_EoP->Fill(myMVAVar_EoP);
0352           h_IoEmIoP->Fill(myMVAVar_IoEmIoP);
0353           h_eleEoPout->Fill(myMVAVar_eleEoPout);
0354           h_EoPout->Fill(myMVAVar_EoPout);
0355           h_PreShowerOverRaw->Fill(myMVAVar_PreShowerOverRaw);
0356           h_pt->Fill(myMVAVar_pt);
0357 
0358           if (debug)
0359             std::cout << "************************* New Good Event:: " << ev << " *************************"
0360                       << std::endl;
0361 
0362           // ********************* Triggering electrons
0363 
0364           bool elePresel = trainTrigPresel(theEGamma[j]);
0365 
0366           double mvaTrigMthd1 = -11.;
0367           double mvaTrigMthd2 = -11.;
0368 
0369           double mvaTrigNonIp = -11.;
0370 
0371           double mvaNonTrigMthd1 = -11;
0372           double mvaNonTrigMthd2 = -11;
0373 
0374           mvaNonTrigMthd1 = myMVANonTrigV0->mvaValue((theEGamma[j]), *pv, thebuilder, lazyTools, debugMVAclass);
0375           mvaNonTrigMthd2 = myMVANonTrigV0->mvaValue(myMVAVar_fbrem,
0376                                                      myMVAVar_kfchi2,
0377                                                      myMVAVar_kfhits,
0378                                                      myMVAVar_gsfchi2,
0379                                                      myMVAVar_deta,
0380                                                      myMVAVar_dphi,
0381                                                      myMVAVar_detacalo,
0382                                                      // myMVAVar_dphicalo,
0383                                                      myMVAVar_see,
0384                                                      myMVAVar_spp,
0385                                                      myMVAVar_etawidth,
0386                                                      myMVAVar_phiwidth,
0387                                                      myMVAVar_e1x5e5x5,
0388                                                      myMVAVar_R9,
0389                                                      //myMVAVar_nbrems,
0390                                                      myMVAVar_HoE,
0391                                                      myMVAVar_EoP,
0392                                                      myMVAVar_IoEmIoP,
0393                                                      myMVAVar_eleEoPout,
0394                                                      myMVAVar_PreShowerOverRaw,
0395                                                      // myMVAVar_EoPout,
0396                                                      myMVAVar_eta,
0397                                                      myMVAVar_pt,
0398                                                      debugMyVar);
0399           h_mva_nonTrig->Fill(mvaNonTrigMthd1);
0400 
0401           if (debug)
0402             std::cout << "Non-Triggering:: MyMVA Method-1 " << mvaNonTrigMthd1 << " MyMVA Method-2 " << mvaNonTrigMthd2
0403                       << std::endl;
0404 
0405           if (elePresel) {
0406             mvaTrigNonIp =
0407                 myMVATrigNoIPV0->mvaValue((theEGamma[j]), *pv, _Rho, /*thebuilder,*/ lazyTools, debugMVAclass);
0408 
0409             mvaTrigMthd1 = myMVATrigV0->mvaValue((theEGamma[j]), *pv, thebuilder, lazyTools, debugMVAclass);
0410             mvaTrigMthd2 = myMVATrigV0->mvaValue(myMVAVar_fbrem,
0411                                                  myMVAVar_kfchi2,
0412                                                  myMVAVar_kfhits,
0413                                                  myMVAVar_gsfchi2,
0414                                                  myMVAVar_deta,
0415                                                  myMVAVar_dphi,
0416                                                  myMVAVar_detacalo,
0417                                                  // myMVAVar_dphicalo,
0418                                                  myMVAVar_see,
0419                                                  myMVAVar_spp,
0420                                                  myMVAVar_etawidth,
0421                                                  myMVAVar_phiwidth,
0422                                                  myMVAVar_e1x5e5x5,
0423                                                  myMVAVar_R9,
0424                                                  //myMVAVar_nbrems,
0425                                                  myMVAVar_HoE,
0426                                                  myMVAVar_EoP,
0427                                                  myMVAVar_IoEmIoP,
0428                                                  myMVAVar_eleEoPout,
0429                                                  myMVAVar_PreShowerOverRaw,
0430                                                  // myMVAVar_EoPout,
0431                                                  myMVAVar_d0,
0432                                                  myMVAVar_ip3d,
0433                                                  myMVAVar_eta,
0434                                                  myMVAVar_pt,
0435                                                  debugMyVar);
0436             h_mva_trig->Fill(mvaTrigMthd1);
0437             h_mva_trigNonIp->Fill(mvaTrigNonIp);
0438           }
0439 
0440           if (debug)
0441             std::cout << "Triggering:: ElePreselection " << elePresel << " MyMVA Method-1 " << mvaTrigMthd1
0442                       << " MyMVA Method-2 " << mvaTrigMthd2 << std::endl;
0443         }
0444       }  // End Loop on RECO electrons
0445     }    // End if MC electrons selection
0446   }      //End Loop Generator Particles
0447 }
0448 // ------------ method called once each job just before starting event loop  ------------
0449 void ElectronTestAnalyzer::myVar(const reco::GsfElectron& ele,
0450                                  const reco::Vertex& vertex,
0451                                  const TransientTrackBuilder& transientTrackBuilder,
0452                                  EcalClusterLazyTools const& myEcalCluster,
0453                                  bool printDebug) {
0454   bool validKF = false;
0455   reco::TrackRef myTrackRef = ele.closestCtfTrackRef();
0456   validKF = (myTrackRef.isAvailable());
0457   validKF = (myTrackRef.isNonnull());
0458 
0459   myMVAVar_fbrem = ele.fbrem();
0460   myMVAVar_kfchi2 = (validKF) ? myTrackRef->normalizedChi2() : 0;
0461   myMVAVar_kfhits = (validKF) ? myTrackRef->hitPattern().trackerLayersWithMeasurement() : -1.;
0462   //  myMVAVar_kfhits          =  (validKF) ? myTrackRef->numberOfValidHits() : -1. ;   // for analysist save also this
0463   myMVAVar_gsfchi2 = ele.gsfTrack()->normalizedChi2();  // to be checked
0464 
0465   myMVAVar_deta = ele.deltaEtaSuperClusterTrackAtVtx();
0466   myMVAVar_dphi = ele.deltaPhiSuperClusterTrackAtVtx();
0467   myMVAVar_detacalo = ele.deltaEtaSeedClusterTrackAtCalo();
0468   myMVAVar_dphicalo = ele.deltaPhiSeedClusterTrackAtCalo();
0469 
0470   myMVAVar_see = ele.sigmaIetaIeta();  //EleSigmaIEtaIEta
0471   const auto& vCov = myEcalCluster.localCovariances(*(ele.superCluster()->seed()));
0472   if (edm::isFinite(vCov[2]))
0473     myMVAVar_spp = sqrt(vCov[2]);  //EleSigmaIPhiIPhi
0474   else
0475     myMVAVar_spp = 0.;
0476   myMVAVar_etawidth = ele.superCluster()->etaWidth();
0477   myMVAVar_phiwidth = ele.superCluster()->phiWidth();
0478   myMVAVar_e1x5e5x5 = (ele.e5x5()) != 0. ? 1. - (ele.e1x5() / ele.e5x5()) : -1.;
0479   myMVAVar_R9 = myEcalCluster.e3x3(*(ele.superCluster()->seed())) / ele.superCluster()->rawEnergy();
0480   myMVAVar_nbrems = std::abs(ele.numberOfBrems());
0481 
0482   myMVAVar_HoE = ele.hadronicOverEm();
0483   myMVAVar_EoP = ele.eSuperClusterOverP();
0484   myMVAVar_IoEmIoP =
0485       (1.0 / ele.ecalEnergy()) - (1.0 / ele.p());  // in the future to be changed with ele.gsfTrack()->p()
0486   myMVAVar_eleEoPout = ele.eEleClusterOverPout();
0487   myMVAVar_EoPout = ele.eSeedClusterOverPout();
0488   myMVAVar_PreShowerOverRaw = ele.superCluster()->preshowerEnergy() / ele.superCluster()->rawEnergy();
0489 
0490   myMVAVar_eta = ele.superCluster()->eta();
0491   myMVAVar_pt = ele.pt();
0492 
0493   //d0
0494   if (ele.gsfTrack().isNonnull()) {
0495     myMVAVar_d0 = (-1.0) * ele.gsfTrack()->dxy(vertex.position());
0496   } else if (ele.closestCtfTrackRef().isNonnull()) {
0497     myMVAVar_d0 = (-1.0) * ele.closestCtfTrackRef()->dxy(vertex.position());
0498   } else {
0499     myMVAVar_d0 = -9999.0;
0500   }
0501 
0502   //default values for IP3D
0503   myMVAVar_ip3d = -999.0;
0504   // myMVAVar_ip3dSig = 0.0;
0505   if (ele.gsfTrack().isNonnull()) {
0506     const double gsfsign = ((-ele.gsfTrack()->dxy(vertex.position())) >= 0) ? 1. : -1.;
0507 
0508     const reco::TransientTrack& tt = transientTrackBuilder.build(ele.gsfTrack());
0509     const std::pair<bool, Measurement1D>& ip3dpv = IPTools::absoluteImpactParameter3D(tt, vertex);
0510     if (ip3dpv.first) {
0511       double ip3d = gsfsign * ip3dpv.second.value();
0512       //double ip3derr = ip3dpv.second.error();
0513       myMVAVar_ip3d = ip3d;
0514       // myMVAVar_ip3dSig = ip3d/ip3derr;
0515     }
0516   }
0517 
0518   if (printDebug) {
0519     std::cout << " My Local Variables " << std::endl;
0520     std::cout << " fbrem " << myMVAVar_fbrem << " kfchi2 " << myMVAVar_kfchi2 << " mykfhits " << myMVAVar_kfhits
0521               << " gsfchi2 " << myMVAVar_gsfchi2 << " deta " << myMVAVar_deta << " dphi " << myMVAVar_dphi
0522               << " detacalo " << myMVAVar_detacalo << " dphicalo " << myMVAVar_dphicalo << " see " << myMVAVar_see
0523               << " spp " << myMVAVar_spp << " etawidth " << myMVAVar_etawidth << " phiwidth " << myMVAVar_phiwidth
0524               << " e1x5e5x5 " << myMVAVar_e1x5e5x5 << " R9 " << myMVAVar_R9 << " mynbrems " << myMVAVar_nbrems
0525               << " HoE " << myMVAVar_HoE << " EoP " << myMVAVar_EoP << " IoEmIoP " << myMVAVar_IoEmIoP << " eleEoPout "
0526               << myMVAVar_eleEoPout << " EoPout " << myMVAVar_EoPout << " PreShowerOverRaw "
0527               << myMVAVar_PreShowerOverRaw << " d0 " << myMVAVar_d0 << " ip3d " << myMVAVar_ip3d << " eta "
0528               << myMVAVar_eta << " pt " << myMVAVar_pt << std::endl;
0529   }
0530   return;
0531 }
0532 void ElectronTestAnalyzer::myBindVariables() {
0533   // this binding is needed for variables that sometime diverge.
0534 
0535   if (myMVAVar_fbrem < -1.)
0536     myMVAVar_fbrem = -1.;
0537 
0538   myMVAVar_deta = std::abs(myMVAVar_deta);
0539   if (myMVAVar_deta > 0.06)
0540     myMVAVar_deta = 0.06;
0541 
0542   myMVAVar_dphi = std::abs(myMVAVar_dphi);
0543   if (myMVAVar_dphi > 0.6)
0544     myMVAVar_dphi = 0.6;
0545 
0546   if (myMVAVar_EoPout > 20.)
0547     myMVAVar_EoPout = 20.;
0548 
0549   if (myMVAVar_EoP > 20.)
0550     myMVAVar_EoP = 20.;
0551 
0552   if (myMVAVar_eleEoPout > 20.)
0553     myMVAVar_eleEoPout = 20.;
0554 
0555   myMVAVar_detacalo = std::abs(myMVAVar_detacalo);
0556   if (myMVAVar_detacalo > 0.2)
0557     myMVAVar_detacalo = 0.2;
0558 
0559   myMVAVar_dphicalo = std::abs(myMVAVar_dphicalo);
0560   if (myMVAVar_dphicalo > 0.4)
0561     myMVAVar_dphicalo = 0.4;
0562 
0563   if (myMVAVar_e1x5e5x5 < -1.)
0564     myMVAVar_e1x5e5x5 = -1;
0565 
0566   if (myMVAVar_e1x5e5x5 > 2.)
0567     myMVAVar_e1x5e5x5 = 2.;
0568 
0569   if (myMVAVar_R9 > 5)
0570     myMVAVar_R9 = 5;
0571 
0572   if (myMVAVar_gsfchi2 > 200.)
0573     myMVAVar_gsfchi2 = 200;
0574 
0575   if (myMVAVar_kfchi2 > 10.)
0576     myMVAVar_kfchi2 = 10.;
0577 
0578   // Needed for a bug in CMSSW_420, fixed in more recent CMSSW versions
0579   if (edm::isNotFinite(myMVAVar_spp))
0580     myMVAVar_spp = 0.;
0581 
0582   return;
0583 }
0584 bool ElectronTestAnalyzer::trainTrigPresel(const reco::GsfElectron& ele) {
0585   bool myTrigPresel = false;
0586   if (std::abs(ele.superCluster()->eta()) < 1.479) {
0587     if (ele.sigmaIetaIeta() < 0.014 && ele.hadronicOverEm() < 0.15 && ele.dr03TkSumPt() / ele.pt() < 0.2 &&
0588         ele.dr03EcalRecHitSumEt() / ele.pt() < 0.2 && ele.dr03HcalTowerSumEt() / ele.pt() < 0.2 &&
0589         ele.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) == 0)
0590       myTrigPresel = true;
0591   } else {
0592     if (ele.sigmaIetaIeta() < 0.035 && ele.hadronicOverEm() < 0.10 && ele.dr03TkSumPt() / ele.pt() < 0.2 &&
0593         ele.dr03EcalRecHitSumEt() / ele.pt() < 0.2 && ele.dr03HcalTowerSumEt() / ele.pt() < 0.2 &&
0594         ele.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) == 0)
0595       myTrigPresel = true;
0596   }
0597 
0598   return myTrigPresel;
0599 }
0600 void ElectronTestAnalyzer::beginJob() { ev = 0; }
0601 
0602 // ------------ method called once each job just after ending the event loop  ------------
0603 void ElectronTestAnalyzer::endJob() { std::cout << " endJob:: #events " << ev << std::endl; }
0604 
0605 void ElectronTestAnalyzer::evaluate_mvas(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0606   edm::Handle<reco::VertexCollection> hVertex;
0607   iEvent.getByToken(vertexToken_, hVertex);
0608   const reco::VertexCollection* pvCol = hVertex.product();
0609 
0610   // FIXME: unused variable giving compilation warnings/errors
0611   //   Handle<double> hRho;
0612   //   iEvent.getByToken(eventrhoToken_,hRho);
0613   //   double Rho = *hRho;
0614   //
0615   //   Handle<reco::PFCandidateCollection> hPfCandProduct;
0616   //    iEvent.getByToken(pfCandToken_, hPfCandProduct);
0617   //   const reco::PFCandidateCollection &inPfCands = *(hPfCandProduct.product());
0618 
0619   edm::Handle<reco::GsfElectronCollection> theEGammaCollection;
0620   iEvent.getByToken(gsfEleToken_, theEGammaCollection);
0621   const reco::GsfElectronCollection inElectrons = *(theEGammaCollection.product());
0622 
0623   edm::Handle<reco::MuonCollection> hMuonProduct;
0624   iEvent.getByToken(muonToken_, hMuonProduct);
0625   const reco::MuonCollection inMuons = *(hMuonProduct.product());
0626 
0627   reco::MuonCollection IdentifiedMuons;
0628   reco::GsfElectronCollection IdentifiedElectrons;
0629 
0630   for (reco::GsfElectronCollection::const_iterator iE = inElectrons.begin(); iE != inElectrons.end(); ++iE) {
0631     double electronTrackZ = 0;
0632     if (iE->gsfTrack().isNonnull()) {
0633       electronTrackZ = iE->gsfTrack()->dz(pvCol->at(0).position());
0634     } else if (iE->closestCtfTrackRef().isNonnull()) {
0635       electronTrackZ = iE->closestCtfTrackRef()->dz(pvCol->at(0).position());
0636     }
0637     if (std::abs(electronTrackZ) > 0.2)
0638       continue;
0639 
0640     if (std::abs(iE->superCluster()->eta()) < 1.479) {
0641       if (iE->pt() > 20) {
0642         if (iE->sigmaIetaIeta() > 0.01)
0643           continue;
0644         if (std::abs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.007)
0645           continue;
0646         if (std::abs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8)
0647           continue;
0648         if (iE->hadronicOverEm() > 0.15)
0649           continue;
0650       } else {
0651         if (iE->sigmaIetaIeta() > 0.012)
0652           continue;
0653         if (std::abs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.007)
0654           continue;
0655         if (std::abs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8)
0656           continue;
0657         if (iE->hadronicOverEm() > 0.15)
0658           continue;
0659       }
0660     } else {
0661       if (iE->pt() > 20) {
0662         if (iE->sigmaIetaIeta() > 0.03)
0663           continue;
0664         if (std::abs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.010)
0665           continue;
0666         if (std::abs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8)
0667           continue;
0668       } else {
0669         if (iE->sigmaIetaIeta() > 0.032)
0670           continue;
0671         if (std::abs(iE->deltaEtaSuperClusterTrackAtVtx()) > 0.010)
0672           continue;
0673         if (std::abs(iE->deltaPhiSuperClusterTrackAtVtx()) > 0.8)
0674           continue;
0675       }
0676     }
0677     IdentifiedElectrons.push_back(*iE);
0678   }
0679 
0680   for (reco::MuonCollection::const_iterator iM = inMuons.begin(); iM != inMuons.end(); ++iM) {
0681     if (!(iM->innerTrack().isNonnull())) {
0682       continue;
0683     }
0684 
0685     if (!(iM->isGlobalMuon() || iM->isTrackerMuon()))
0686       continue;
0687     if (iM->innerTrack()->numberOfValidHits() < 11)
0688       continue;
0689 
0690     IdentifiedMuons.push_back(*iM);
0691   }
0692 
0693   EcalClusterLazyTools lazyTools(
0694       iEvent, ecalClusterESGetTokens_.get(iSetup), reducedEBRecHitCollectionToken_, reducedEERecHitCollectionToken_);
0695 
0696   edm::ESHandle<TransientTrackBuilder> builder;
0697   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
0698   TransientTrackBuilder thebuilder = *(builder.product());
0699 
0700   for (reco::GsfElectronCollection::const_iterator iE = inElectrons.begin(); iE != inElectrons.end(); ++iE) {
0701     reco::GsfElectron ele = *iE;
0702 
0703     double idmva = myMVATrigV0->mvaValue(ele, pvCol->at(0), thebuilder, lazyTools);
0704 
0705     std::cout << "idmva = " << idmva << std::endl;
0706   }
0707 }
0708 
0709 //define this as a plug-in
0710 DEFINE_FWK_MODULE(ElectronTestAnalyzer);
[ 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!