ElectronIdentification/src/ElectronMVAEstimator.cc

0001 #include "RecoEgamma/ElectronIdentification/interface/ElectronMVAEstimator.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0005 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0006 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0007
0008 ElectronMVAEstimator::ElectronMVAEstimator() : cfg_{} {}
0009
0010 ElectronMVAEstimator::ElectronMVAEstimator(const std::string& fileName) : cfg_{} {
0011   // Taken from Daniele (his mail from the 30/11)
0012   //  tmvaReader.BookMVA("BDTSimpleCat","../Training/weights_Root527b_3Depth_DanVarConvRej_2PtBins_10Pt_800TPrune5_Min100Events_NoBjets_half/TMVA_BDTSimpleCat.weights.xm");
0013   // training of the 7/12 with Nvtx added
0014   gbr_.push_back(createGBRForest(fileName));
0015 }
0016
0017 ElectronMVAEstimator::ElectronMVAEstimator(const Configuration& cfg) : cfg_(cfg) {
0018   for (const auto& weightsfile : cfg_.vweightsfiles) {
0019     gbr_.push_back(createGBRForest(weightsfile));
0020   }
0021 }
0022
0023 double ElectronMVAEstimator::mva(const reco::GsfElectron& myElectron, int nvertices) const {
0024   float vars[18];
0025
0026   vars[0] = myElectron.fbrem();
0027   vars[1] = std::abs(myElectron.deltaEtaSuperClusterTrackAtVtx());
0028   vars[2] = std::abs(myElectron.deltaPhiSuperClusterTrackAtVtx());
0029   vars[3] = myElectron.sigmaIetaIeta();
0030   vars[4] = myElectron.hcalOverEcal();
0031   vars[5] = myElectron.eSuperClusterOverP();
0032   vars[6] = (myElectron.e5x5()) != 0. ? 1. - (myElectron.e1x5() / myElectron.e5x5()) : -1.;
0033   vars[7] = myElectron.eEleClusterOverPout();
0034   vars[8] = std::abs(myElectron.deltaEtaEleClusterTrackAtCalo());
0035
0036   bool validKF = false;
0037
0038   reco::TrackRef myTrackRef = myElectron.closestCtfTrackRef();
0039   validKF = (myTrackRef.isNonnull() && myTrackRef.isAvailable());
0040
0041   vars[9] = (validKF) ? myTrackRef->normalizedChi2() : 0;
0042   vars[10] = (validKF) ? myTrackRef->hitPattern().trackerLayersWithMeasurement() : -1.;
0043   vars[11] = myElectron.gsfTrack()->missingInnerHits();
0044   vars[12] = std::abs(myElectron.convDist());
0045   vars[13] = std::abs(myElectron.convDcot());
0046   vars[14] = nvertices;
0047   vars[15] = myElectron.eta();
0048   vars[16] = myElectron.pt();
0049   vars[17] = myElectron.ecalDrivenSeed();
0050
0051   bindVariables(vars);
0052
0053   //0 pt &lt; 10 &amp;&amp; abs(eta)&lt;=1.485
0054   //1 pt &gt;= 10 &amp;&amp; abs(eta)&lt;=1.485
0055   //2 pt &lt; 10 &amp;&amp; abs(eta)&gt; 1.485
0056   //3 pt &gt;= 10 &amp;&amp;  abs(eta)&gt; 1.485
0057
0058   const unsigned index = (unsigned)(myElectron.pt() >= 10) + 2 * (unsigned)(std::abs(myElectron.eta()) > 1.485);
0059
0060   double result = gbr_[index]->GetAdaBoostClassifier(vars);
0061   //
0062   //  std::cout << "fbrem" << vars[0] << std::endl;
0063   //  std::cout << "detain"<< vars[1] << std::endl;
0064   //  std::cout << "dphiin"<< vars[2] << std::endl;
0065   //  std::cout << "sieie"<< vars[3] << std::endl;
0066   //  std::cout << "hoe"<< vars[4] << std::endl;
0067   //  std::cout << "eop"<< vars[5] << std::endl;
0068   //  std::cout << "e1x5e5x5"<< vars[6] << std::endl;
0069   //  std::cout << "eleopout"<< vars[7] << std::endl;
0070   //  std::cout << "detaeleout"<< vars[8] << std::endl;
0071   //  std::cout << "kfchi2"<< vars[9] << std::endl;
0072   //  std::cout << "kfhits"<< vars[10] << std::endl;
0073   //  std::cout << "mishits"<<vars[11] << std::endl;
0074   //  std::cout << "dist"<< vars[12] << std::endl;
0075   //  std::cout << "dcot"<< vars[13] << std::endl;
0076   //  std::cout << "nvtx"<< vars[14] << std::endl;
0077   //  std::cout << "eta"<< vars[15] << std::endl;
0078   //  std::cout << "pt"<< vars[16] << std::endl;
0079   //  std::cout << "ecalseed"<< vars[17] << std::endl;
0080   //
0081   //  std::cout << " MVA " << result << std::endl;
0082   return result;
0083 }
0084
0085 void ElectronMVAEstimator::bindVariables(float vars[18]) const {
0086   if (vars[0] < -1.)
0087     vars[1] = -1.;
0088
0089   if (vars[1] > 0.06)
0090     vars[1] = 0.06;
0091
0092   if (vars[2] > 0.6)
0093     vars[2] = 0.6;
0094
0095   if (vars[5] > 20.)
0096     vars[5] = 20.;
0097
0098   if (vars[7] > 20.)
0099     vars[7] = 20;
0100
0101   if (vars[8] > 0.2)
0102     vars[8] = 0.2;
0103
0104   if (vars[9] < 0.)
0105     vars[9] = 0.;
0106
0107   if (vars[9] > 15.)
0108     vars[9] = 15.;
0109
0110   if (vars[6] < -1.)
0111     vars[6] = -1;
0112
0113   if (vars[6] > 2.)
0114     vars[6] = 2.;
0115
0116   if (vars[12] > 15.)
0117     vars[12] = 15.;
0118
0119   if (vars[13] > 3.)
0120     vars[13] = 3.;
0121 }