ElectronIdentification/src/SoftElectronMVAEstimator.cc

0001 #include "RecoEgamma/ElectronIdentification/interface/SoftElectronMVAEstimator.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 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 SoftElectronMVAEstimator::SoftElectronMVAEstimator(const Configuration& cfg) : cfg_(cfg) {
0010   //Check number of weight files given
0011   if (ExpectedNBins != cfg_.vweightsfiles.size()) {
0012     edm::LogError("Soft Electron MVA Error")
0013         << "Expected Number of bins = " << ExpectedNBins
0014         << " does not equal to weightsfiles.size() = " << cfg_.vweightsfiles.size() << std::endl;
0015   }
0016
0017   for (auto& weightsfile : cfg_.vweightsfiles) {
0018     // Taken from Daniele (his mail from the 30/11)
0019     // training of the 7/12 with Nvtx added
0020     gbr_.push_back(createGBRForest(weightsfile));
0021   }
0022 }
0023
0024 SoftElectronMVAEstimator::~SoftElectronMVAEstimator() {}
0025
0026 double SoftElectronMVAEstimator::mva(const reco::GsfElectron& myElectron, const reco::VertexCollection& pvc) const {
0027   float vars[25];
0028
0029   vars[0] = myElectron.fbrem();                // fbrem
0030   vars[1] = myElectron.eSuperClusterOverP();   //EtotOvePin
0031   vars[2] = myElectron.eEleClusterOverPout();  //eleEoPout
0032
0033   float etot = myElectron.eSuperClusterOverP() * myElectron.trackMomentumAtVtx().R();
0034   float eEcal = myElectron.eEleClusterOverPout() * myElectron.trackMomentumAtEleClus().R();
0035   float dP = myElectron.trackMomentumAtVtx().R() - myElectron.trackMomentumAtEleClus().R();
0036   vars[3] = (etot - eEcal) / dP;                         //EBremOverDeltaP
0037   vars[4] = std::log(myElectron.sigmaEtaEta());          //logSigmaEtaEta
0038   vars[5] = myElectron.deltaEtaEleClusterTrackAtCalo();  //DeltaEtaTrackEcalSeed
0039   vars[6] = myElectron.hcalOverEcalBc();                 //HoE
0040
0041   bool validKF = false;
0042   reco::TrackRef myTrackRef = myElectron.closestCtfTrackRef();
0043   validKF = (myTrackRef.isAvailable() && myTrackRef.isNonnull());
0044   vars[7] = myElectron.gsfTrack()->normalizedChi2();                                    //gsfchi2
0045   vars[8] = (validKF) ? myTrackRef->normalizedChi2() : 0;                               //kfchi2
0046   vars[9] = (validKF) ? myTrackRef->hitPattern().trackerLayersWithMeasurement() : -1.;  //kfhits
0047
0048   vars[10] = myElectron.gsfTrack().get()->ptModeError() / myElectron.gsfTrack().get()->ptMode();     //SigmaPtOverPt
0049   vars[11] = myElectron.deltaEtaSuperClusterTrackAtVtx();                                            //deta
0050   vars[12] = myElectron.deltaPhiSuperClusterTrackAtVtx();                                            //dphi
0051   vars[13] = myElectron.deltaEtaSeedClusterTrackAtCalo();                                            //detacalo
0052   vars[14] = myElectron.sigmaIetaIeta();                                                             //see
0053   vars[15] = myElectron.sigmaIphiIphi();                                                             //spp
0054   vars[16] = myElectron.r9();                                                                        //R9
0055   vars[17] = myElectron.superCluster()->etaWidth();                                                  //etawidth
0056   vars[18] = myElectron.superCluster()->phiWidth();                                                  //phiwidth
0057   vars[19] = (myElectron.e5x5()) != 0. ? 1. - (myElectron.e1x5() / myElectron.e5x5()) : -1.;         //OneMinusE1x5E5x5
0058   vars[20] = (1.0 / myElectron.ecalEnergy()) - (1.0 / myElectron.p());                               // IoEmIoP
0059   vars[21] = myElectron.superCluster()->preshowerEnergy() / myElectron.superCluster()->rawEnergy();  //PreShowerOverRaw
0060   vars[22] = pvc.size();                                                                             // nPV
0061   vars[23] = myElectron.pt();                                                                        //pt
0062   vars[24] = myElectron.eta();                                                                       //eta
0063
0064   /*
0065   std::cout<<"fbrem "<<fbrem<<std::endl;
0066   std::cout<<"EtotOvePin "<<EtotOvePin<<std::endl;
0067   std::cout<<"eleEoPout "<<eleEoPout<<std::endl;
0068   std::cout<<"EBremOverDeltaP "<<EBremOverDeltaP<<std::endl;
0069   std::cout<<"logSigmaEtaEta "<<logSigmaEtaEta<<std::endl;
0070   std::cout<<"DeltaEtaTrackEcalSeed "<<DeltaEtaTrackEcalSeed<<std::endl;
0071   std::cout<<"HoE "<<HoE<<std::endl;
0072   std::cout<<"kfchi2 "<<kfchi2<<std::endl;
0073   std::cout<<"kfhits "<<kfhits<<std::endl;
0074   std::cout<<"gsfchi2 "<<gsfchi2<<std::endl;
0075   std::cout<<"SigmaPtOverPt "<<SigmaPtOverPt<<std::endl;
0076   std::cout<<"deta "<<deta<<std::endl;
0077   std::cout<<"dphi "<<dphi<<std::endl;
0078   std::cout<<"detacalo "<<detacalo<<std::endl;
0079   std::cout<<"see "<<see<<std::endl;
0080   std::cout<< "spp "             <<          spp<< std::endl;
0081   std::cout<< "R9 "             <<         R9<< std::endl;
0082   std::cout<< "IoEmIoP "        <<         IoEmIoP<< std::endl;
0083   std::cout<<"etawidth "<<etawidth<<std::endl;
0084   std::cout<<"phiwidth "<<phiwidth<<std::endl;
0085   std::cout<<"OneMinusE1x5E5x5 "<<OneMinusE1x5E5x5<<std::endl;
0086   std::cout<<"PreShowerOverRaw "<<PreShowerOverRaw<<std::endl;
0087 */
0088   bindVariables(vars);
0089
0090   double result = gbr_[0]->GetClassifier(vars);
0091
0092   return result;
0093 }
0094
0095 void SoftElectronMVAEstimator::bindVariables(float vars[25]) const {
0096   if (vars[0] < -1.)  //fbrem
0097     vars[0] = -1.;
0098
0099   vars[11] = std::abs(vars[11]);  // deta
0100   if (vars[11] > 0.06)
0101     vars[11] = 0.06;
0102
0103   vars[12] = std::abs(vars[12]);
0104   if (vars[12] > 0.6)
0105     vars[12] = 0.6;
0106
0107   //if(EoP > 20.)
0108   //  EoP = 20.;
0109
0110   if (vars[2] > 20.)  //eleEoPout
0111     vars[2] = 20.;
0112
0113   vars[13] = std::abs(vars[13]);  //detacalo
0114   if (vars[13] > 0.2)
0115     vars[13] = 0.2;
0116
0117   if (vars[19] < -1.)  //OneMinusE1x5E5x5
0118     vars[19] = -1;
0119
0120   if (vars[19] > 2.)  //OneMinusE1x5E5x5
0121     vars[19] = 2.;
0122
0123   if (vars[7] > 200.)  //gsfchi2
0124     vars[7] = 200;
0125
0126   if (vars[8] > 10.)  //kfchi2
0127     vars[8] = 10.;
0128 }