File indexing completed on 2021-04-01 02:23:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <memory>
0018
0019
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
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
0089 edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0090
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
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
0164
0165
0166
0167
0168
0169
0170
0171
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
0178 vertexToken_(consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"))),
0179
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
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
0257
0258 }
0259
0260
0261
0262
0263
0264
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
0275
0276
0277
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 {
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
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
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
0383 myMVAVar_see,
0384 myMVAVar_spp,
0385 myMVAVar_etawidth,
0386 myMVAVar_phiwidth,
0387 myMVAVar_e1x5e5x5,
0388 myMVAVar_R9,
0389
0390 myMVAVar_HoE,
0391 myMVAVar_EoP,
0392 myMVAVar_IoEmIoP,
0393 myMVAVar_eleEoPout,
0394 myMVAVar_PreShowerOverRaw,
0395
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, 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
0418 myMVAVar_see,
0419 myMVAVar_spp,
0420 myMVAVar_etawidth,
0421 myMVAVar_phiwidth,
0422 myMVAVar_e1x5e5x5,
0423 myMVAVar_R9,
0424
0425 myMVAVar_HoE,
0426 myMVAVar_EoP,
0427 myMVAVar_IoEmIoP,
0428 myMVAVar_eleEoPout,
0429 myMVAVar_PreShowerOverRaw,
0430
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 }
0445 }
0446 }
0447 }
0448
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
0463 myMVAVar_gsfchi2 = ele.gsfTrack()->normalizedChi2();
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();
0471 const auto& vCov = myEcalCluster.localCovariances(*(ele.superCluster()->seed()));
0472 if (edm::isFinite(vCov[2]))
0473 myMVAVar_spp = sqrt(vCov[2]);
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());
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
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
0503 myMVAVar_ip3d = -999.0;
0504
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
0513 myMVAVar_ip3d = ip3d;
0514
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
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
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
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
0611
0612
0613
0614
0615
0616
0617
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
0710 DEFINE_FWK_MODULE(ElectronTestAnalyzer);