ElectronTools/test/ElectronTestAnalyzer.cc

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