CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoTauTag/RecoTau/src/AntiElectronIDMVA.cc

Go to the documentation of this file.
00001 #include <TFile.h>
00002 #include <TMath.h>
00003 #include "RecoTauTag/RecoTau//interface/AntiElectronIDMVA.h"
00004 #include "RecoTauTag/RecoTau/interface/TMVAZipReader.h"
00005 
00006 AntiElectronIDMVA::AntiElectronIDMVA():
00007   isInitialized_(kFALSE),
00008   methodName_("BDT")
00009 {
00010   for(UInt_t i=0; i<6; ++i) {
00011     fTMVAReader_[i] = 0;
00012   }
00013 
00014   verbosity_ = 1;
00015 }
00016 
00017 
00018 AntiElectronIDMVA::~AntiElectronIDMVA()
00019 {
00020   for(UInt_t i=0; i<6; ++i) {
00021     if (fTMVAReader_[i]) delete fTMVAReader_[i];
00022   }
00023 }
00024 
00025 void AntiElectronIDMVA::Initialize(std::string methodName,
00026                                    std::string oneProng0Pi0_BL,
00027                                    std::string oneProng1pi0wGSF_BL,
00028                                    std::string oneProng1pi0woGSF_BL,
00029                                    std::string oneProng0Pi0_EC,
00030                                    std::string oneProng1pi0wGSF_EC,
00031                                    std::string oneProng1pi0woGSF_EC
00032                                    ){
00033   TauSignalPFGammaCands_ = 0.; 
00034   TauVisMass_ = 0.; 
00035   GammadEta_ = 0.; 
00036   GammadPhi_ = 0.; 
00037   GammadPt_ = 0.;
00038   TauLeadPFChargedHadrMva_ = 0.;
00039   TauLeadPFChargedHadrHoP_ = 0.;
00040   TauLeadPFChargedHadrEoP_ = 0.;
00041   TauEmFraction_ = 0.;
00042 
00043   for(UInt_t i=0; i<6; ++i) {
00044     if (fTMVAReader_[i]) delete fTMVAReader_[i];
00045   }
00046 
00047   isInitialized_ = kTRUE;
00048   methodName_    = methodName;
00049 
00050   //TMVA::Tools::Instance();
00051 
00052   TMVA::Reader *readerX0BL = new TMVA::Reader( "!Color:Silent:Error" );
00053   readerX0BL->AddVariable("HoP",       &TauLeadPFChargedHadrHoP_);
00054   readerX0BL->AddVariable("EoP",       &TauLeadPFChargedHadrEoP_);
00055   //readerX0BL->AddVariable("emFraction",&TauEmFraction_);
00056   readerX0BL->SetVerbose(verbosity_);
00057   reco::details::loadTMVAWeights(readerX0BL,  methodName_, oneProng0Pi0_BL );
00058 
00059   TMVA::Reader *reader11BL = new TMVA::Reader( "!Color:Silent:Error" );
00060   reader11BL->AddVariable("mva",               &TauLeadPFChargedHadrMva_);
00061   reader11BL->AddVariable("visMass",           &TauVisMass_);
00062   reader11BL->AddVariable("etaMom2*TMath::Sqrt(gammaFrac)*pt", &GammadEta_);
00063   reader11BL->AddVariable("phiMom2*TMath::Sqrt(gammaFrac)*pt", &GammadPhi_);
00064   reader11BL->AddVariable("gammaFrac",         &GammadPt_);
00065   reader11BL->SetVerbose(verbosity_);
00066   reco::details::loadTMVAWeights(reader11BL,  methodName_, oneProng1pi0wGSF_BL );
00067 
00068   TMVA::Reader *reader01BL = new TMVA::Reader( "!Color:Silent:Error" );
00069   reader01BL->AddVariable("visMass",           &TauVisMass_);
00070   reader01BL->AddVariable("etaMom2*TMath::Sqrt(gammaFrac)*pt", &GammadEta_);
00071   reader01BL->AddVariable("phiMom2*TMath::Sqrt(gammaFrac)*pt", &GammadPhi_);
00072   reader01BL->AddVariable("gammaFrac",         &GammadPt_);
00073   reader01BL->SetVerbose(verbosity_);
00074   reco::details::loadTMVAWeights(reader01BL,  methodName_, oneProng1pi0woGSF_BL );
00075 
00077 
00078   TMVA::Reader *readerX0EC = new TMVA::Reader( "!Color:Silent:Error" );
00079   readerX0EC->AddVariable("HoP",       &TauLeadPFChargedHadrHoP_);
00080   readerX0EC->AddVariable("EoP",       &TauLeadPFChargedHadrEoP_);
00081   //readerX0EC->AddVariable("emFraction",&TauEmFraction_);
00082   readerX0EC->SetVerbose(verbosity_);
00083   reco::details::loadTMVAWeights(readerX0EC,  methodName_, oneProng0Pi0_EC );
00084 
00085   TMVA::Reader *reader11EC = new TMVA::Reader( "!Color:Silent:Error" );
00086   reader11EC->AddVariable("mva",               &TauLeadPFChargedHadrMva_);
00087   reader11EC->AddVariable("visMass",           &TauVisMass_);
00088   reader11EC->AddVariable("etaMom2*TMath::Sqrt(gammaFrac)*pt", &GammadEta_);
00089   reader11EC->AddVariable("phiMom2*TMath::Sqrt(gammaFrac)*pt", &GammadPhi_);
00090   reader11EC->AddVariable("gammaFrac",         &GammadPt_);
00091   reader11EC->SetVerbose(verbosity_);
00092   reco::details::loadTMVAWeights(reader11EC,  methodName_, oneProng1pi0wGSF_EC );
00093 
00094   TMVA::Reader *reader01EC = new TMVA::Reader( "!Color:Silent:Error" );
00095   reader01EC->AddVariable("visMass",           &TauVisMass_);
00096   reader01EC->AddVariable("etaMom2*TMath::Sqrt(gammaFrac)*pt", &GammadEta_);
00097   reader01EC->AddVariable("phiMom2*TMath::Sqrt(gammaFrac)*pt", &GammadPhi_);
00098   reader01EC->AddVariable("gammaFrac",         &GammadPt_);
00099   reader01EC->SetVerbose(verbosity_);
00100   reco::details::loadTMVAWeights(reader01EC,  methodName_, oneProng1pi0woGSF_EC );
00101 
00102   fTMVAReader_[0] = readerX0BL;
00103   fTMVAReader_[1] = reader11BL;
00104   fTMVAReader_[2] = reader01BL;
00105   fTMVAReader_[3] = readerX0EC;
00106   fTMVAReader_[4] = reader11EC;
00107   fTMVAReader_[5] = reader01EC;
00108 }
00109 
00110 
00111 double AntiElectronIDMVA::MVAValue(Float_t TauEta, Float_t TauPt,
00112                                    Float_t TauSignalPFChargedCands, Float_t TauSignalPFGammaCands,
00113                                    Float_t TauLeadPFChargedHadrMva,
00114                                    Float_t TauLeadPFChargedHadrHoP, Float_t TauLeadPFChargedHadrEoP,
00115                                    Float_t TauHasGsf, Float_t TauVisMass,  Float_t TauEmFraction,
00116                                    std::vector<Float_t>* GammasdEta, std::vector<Float_t>* GammasdPhi, std::vector<Float_t>* GammasPt
00117                                    ){
00118 
00119   if (!isInitialized_) {
00120     std::cout << "Error: AntiElectronMVA with method 1 not properly initialized.\n";
00121     return -999;
00122   }
00123 
00124   double mva;
00125 
00126   TauVisMass_              = TauVisMass;
00127   TauLeadPFChargedHadrMva_ = TMath::Max(TauLeadPFChargedHadrMva,float(-1.0));
00128   TauLeadPFChargedHadrHoP_ = TauLeadPFChargedHadrHoP;
00129   TauLeadPFChargedHadrEoP_ = TauLeadPFChargedHadrEoP;
00130   TauEmFraction_           = TMath::Max(TauEmFraction,float(0.0));
00131 
00132   float sumPt  = 0;
00133   float dEta   = 0;
00134   float dEta2  = 0;
00135   float dPhi   = 0;
00136   float dPhi2  = 0;
00137   float sumPt2 = 0;
00138 
00139   for(unsigned int k = 0 ; k < GammasPt->size() ; k++){
00140     float pt_k  = (*GammasPt)[k];
00141     float phi_k = (*GammasdPhi)[k];
00142     if ((*GammasdPhi)[k] > TMath::Pi()) phi_k = (*GammasdPhi)[k] - 2*TMath::Pi();
00143     else if((*GammasdPhi)[k] < -TMath::Pi()) phi_k = (*GammasdPhi)[k] + 2*TMath::Pi();
00144     float eta_k = (*GammasdEta)[k];
00145     sumPt  +=  pt_k;
00146     sumPt2 += (pt_k*pt_k);
00147     dEta   += (pt_k*eta_k);
00148     dEta2  += (pt_k*eta_k*eta_k);
00149     dPhi   += (pt_k*phi_k);
00150     dPhi2  += (pt_k*phi_k*phi_k);
00151   }
00152 
00153   GammadPt_ = sumPt/TauPt;
00154 
00155   if(sumPt>0){
00156     dEta  /= sumPt;
00157     dPhi  /= sumPt;
00158     dEta2 /= sumPt;
00159     dPhi2 /= sumPt;
00160 
00161   }
00162 
00163   //GammadEta_ = dEta;
00164   //GammadPhi_ = dPhi;
00165 
00166   GammadEta_ = TMath::Sqrt(dEta2)*TMath::Sqrt(GammadPt_)*TauPt;
00167   GammadPhi_ = TMath::Sqrt(dPhi2)*TMath::Sqrt(GammadPt_)*TauPt;
00168 
00169 
00170   if( TauSignalPFChargedCands==3 )
00171     mva = 1.0;
00172   else if( TauSignalPFChargedCands==1 && TauSignalPFGammaCands==0){
00173     if(TMath::Abs(TauEta)<1.5)
00174       mva = fTMVAReader_[0]->EvaluateMVA( methodName_ );
00175     else
00176       mva = fTMVAReader_[3]->EvaluateMVA( methodName_ );
00177   }
00178   else if( TauSignalPFChargedCands==1 && TauSignalPFGammaCands>0 && TauHasGsf>0.5){
00179     if(TMath::Abs(TauEta)<1.5)
00180       mva = fTMVAReader_[1]->EvaluateMVA( methodName_ );
00181     else
00182       mva = fTMVAReader_[4]->EvaluateMVA( methodName_ );
00183   }
00184   else if( TauSignalPFChargedCands==1 && TauSignalPFGammaCands>0 && TauHasGsf<0.5){
00185     if(TMath::Abs(TauEta)<1.5)
00186       mva = fTMVAReader_[2]->EvaluateMVA( methodName_ );
00187     else
00188       mva = fTMVAReader_[5]->EvaluateMVA( methodName_ );
00189   }
00190   else{
00191     mva = -99;
00192   }
00193 
00194   return mva;
00195 
00196 }
00197 
00198 double AntiElectronIDMVA::MVAValue(Float_t TauEta, Float_t TauPt,
00199                                    Float_t TauSignalPFChargedCands, Float_t TauSignalPFGammaCands,
00200                                    Float_t TauLeadPFChargedHadrMva,
00201                                    Float_t TauLeadPFChargedHadrHoP, Float_t TauLeadPFChargedHadrEoP,
00202                                    Float_t TauHasGsf, Float_t TauVisMass,  Float_t TauEmFraction,
00203                                    Float_t GammaEtaMom, Float_t GammaPhiMom, Float_t GammaEnFrac
00204                                    ){
00205 
00206   if (!isInitialized_) {
00207     std::cout << "Error: AntiElectronMVA with method 2 not properly initialized.\n";
00208     return -999;
00209   }
00210 
00211 
00212   double mva;
00213 
00214   TauVisMass_              = TauVisMass;
00215   TauLeadPFChargedHadrMva_ = TMath::Max(TauLeadPFChargedHadrMva,float(-1.0));
00216   TauLeadPFChargedHadrHoP_ = TauLeadPFChargedHadrHoP;
00217   TauLeadPFChargedHadrEoP_ = TauLeadPFChargedHadrEoP;
00218   TauEmFraction_           = TMath::Max(TauEmFraction,float(0.0));
00219   GammadPt_                = GammaEnFrac;
00220   GammadEta_               = GammaEtaMom;
00221   GammadPhi_               = GammaPhiMom;
00222 
00223   if( TauSignalPFChargedCands==3 )
00224     mva = 1.0;
00225   else if( TauSignalPFChargedCands==1 && TauSignalPFGammaCands==0){
00226     if(TMath::Abs(TauEta)<1.5)
00227       mva = fTMVAReader_[0]->EvaluateMVA( methodName_ );
00228     else
00229       mva = fTMVAReader_[3]->EvaluateMVA( methodName_ );
00230   }
00231   else if( TauSignalPFChargedCands==1 && TauSignalPFGammaCands>0 && TauHasGsf>0.5){
00232     if(TMath::Abs(TauEta)<1.5)
00233       mva = fTMVAReader_[1]->EvaluateMVA( methodName_ );
00234     else
00235       mva = fTMVAReader_[4]->EvaluateMVA( methodName_ );
00236   }
00237   else if( TauSignalPFChargedCands==1 && TauSignalPFGammaCands>0 && TauHasGsf<0.5){
00238     if(TMath::Abs(TauEta)<1.5)
00239       mva = fTMVAReader_[2]->EvaluateMVA( methodName_ );
00240     else
00241       mva = fTMVAReader_[5]->EvaluateMVA( methodName_ );
00242   }
00243   else{
00244     mva = -99.;
00245   }
00246 
00247   return mva;
00248 
00249 }
00250 
00251 
00252 
00253 double AntiElectronIDMVA::MVAValue(const reco::PFTauRef& thePFTauRef){
00254 
00255   if (!isInitialized_) {
00256     std::cout << "Error: AntiElectronMVA with method 3 not properly initialized.\n";
00257     return -999;
00258   }
00259 
00260   double mva;
00261 
00262   TauVisMass_              = (*thePFTauRef).mass();
00263   TauLeadPFChargedHadrMva_ = TMath::Max((*thePFTauRef).electronPreIDOutput(),float(-1.0));
00264   TauLeadPFChargedHadrHoP_ = ((*thePFTauRef).leadPFChargedHadrCand())->hcalEnergy()/(*thePFTauRef).leadPFChargedHadrCand()->p();
00265   TauLeadPFChargedHadrEoP_ = ((*thePFTauRef).leadPFChargedHadrCand())->ecalEnergy()/(*thePFTauRef).leadPFChargedHadrCand()->p();
00266   TauEmFraction_           = TMath::Max((*thePFTauRef).emFraction(),float(0.0));
00267 
00268   std::vector<float> GammasdEta;
00269   std::vector<float> GammasdPhi;
00270   std::vector<float> GammasPt;
00271 
00272   for(unsigned int k = 0 ; k < ((*thePFTauRef).signalPFGammaCands()).size() ; k++){
00273     reco::PFCandidateRef gamma = ((*thePFTauRef).signalPFGammaCands()).at(k);
00274     if( ((*thePFTauRef).leadPFChargedHadrCand()).isNonnull() ){
00275       GammasdEta.push_back( gamma->eta() - (*thePFTauRef).leadPFChargedHadrCand()->eta() );
00276       GammasdPhi.push_back( gamma->phi() - (*thePFTauRef).leadPFChargedHadrCand()->phi() );
00277     }
00278     else{
00279       GammasdEta.push_back( gamma->eta() - (*thePFTauRef).eta() );
00280       GammasdPhi.push_back( gamma->phi() - (*thePFTauRef).phi() );
00281     }
00282     GammasPt.push_back(  gamma->pt() );
00283   }
00284 
00285   float sumPt  = 0;
00286   float dEta   = 0;
00287   float dEta2  = 0;
00288   float dPhi   = 0;
00289   float dPhi2  = 0;
00290   float sumPt2 = 0;
00291 
00292   for(unsigned int k = 0 ; k < GammasPt.size() ; k++){
00293     float pt_k  = GammasPt[k];
00294     float phi_k = GammasdPhi[k];
00295     if (GammasdPhi[k] > TMath::Pi()) phi_k = GammasdPhi[k] - 2*TMath::Pi();
00296     else if(GammasdPhi[k] < -TMath::Pi()) phi_k = GammasdPhi[k] + 2*TMath::Pi();
00297     float eta_k = GammasdEta[k];
00298     sumPt  +=  pt_k;
00299     sumPt2 += (pt_k*pt_k);
00300     dEta   += (pt_k*eta_k);
00301     dEta2  += (pt_k*eta_k*eta_k);
00302     dPhi   += (pt_k*phi_k);
00303     dPhi2  += (pt_k*phi_k*phi_k);
00304   }
00305 
00306   GammadPt_ = sumPt/(*thePFTauRef).pt();
00307 
00308   if(sumPt>0){
00309     dEta  /= sumPt;
00310     dPhi  /= sumPt;
00311     dEta2 /= sumPt;
00312     dPhi2 /= sumPt;
00313   }
00314 
00315   //GammadEta_ = dEta;
00316   //GammadPhi_ = dPhi;
00317 
00318   GammadEta_ = TMath::Sqrt(dEta2)*TMath::Sqrt(GammadPt_)*(*thePFTauRef).pt();
00319   GammadPhi_ = TMath::Sqrt(dPhi2)*TMath::Sqrt(GammadPt_)*(*thePFTauRef).pt();
00320 
00321   if( ((*thePFTauRef).signalPFChargedHadrCands()).size() == 3)
00322     mva = 1.0;
00323   else if( ((*thePFTauRef).signalPFChargedHadrCands()).size()==1 && ((*thePFTauRef).signalPFGammaCands()).size()==0){
00324     if(TMath::Abs((*thePFTauRef).eta())<1.5)
00325       mva = fTMVAReader_[0]->EvaluateMVA( methodName_ );
00326     else
00327       mva = fTMVAReader_[3]->EvaluateMVA( methodName_ );
00328   }
00329   else if( ((*thePFTauRef).signalPFChargedHadrCands()).size()==1 && ((*thePFTauRef).signalPFGammaCands()).size()>0 && (((*thePFTauRef).leadPFChargedHadrCand())->gsfTrackRef()).isNonnull()){
00330     if(TMath::Abs((*thePFTauRef).eta())<1.5)
00331       mva = fTMVAReader_[1]->EvaluateMVA( methodName_ );
00332     else
00333       mva = fTMVAReader_[4]->EvaluateMVA( methodName_ );
00334   }
00335   else if( ((*thePFTauRef).signalPFChargedHadrCands()).size()==1 && ((*thePFTauRef).signalPFGammaCands()).size()>0 && !(((*thePFTauRef).leadPFChargedHadrCand())->gsfTrackRef()).isNonnull()){
00336     if(TMath::Abs((*thePFTauRef).eta())<1.5)
00337       mva = fTMVAReader_[2]->EvaluateMVA( methodName_ );
00338     else
00339       mva = fTMVAReader_[5]->EvaluateMVA( methodName_ );
00340   }
00341   else{
00342     mva = -99;
00343   }
00344 
00345   return mva;
00346 
00347 }