Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.13 2012/05/29 21:20:12 mingyang Exp $

#include "MitPhysics/Utils/interface/PhotonTools.h"
#include "MitPhysics/Utils/interface/MVATools.h"
#include "MitPhysics/Utils/interface/ElectronTools.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitAna/DataTree/interface/PFCandidateCol.h"
#include "MitAna/DataTree/interface/StableData.h"
#include <TMath.h> 
#include <TFile.h>
#include <TRandom3.h>
#include "TMVA/Tools.h"//MVA
#include "TMVA/Reader.h"//MVA


ClassImp(mithep::MVATools)

using namespace mithep;


//--------------------------------------------------------------------------------------------------
MVATools::MVATools():
  fReaderEndcap(0),
  fReaderBarrel(0),
  
  //MVA Variables v4
  HoE(0),
  covIEtaIEta(0),
  tIso1abs(0),
  tIso3abs(0),
  tIso2abs(0),
  R9(0),
  
  absIsoEcal(0),
  absIsoHcal(0),
  RelEMax(0),
  RelETop(0),
  RelEBottom(0),
  RelELeft(0),
  RelERight(0),
  RelE2x5Max(0),
  RelE2x5Top(0),
  RelE2x5Bottom(0),
  RelE2x5Left(0),
  RelE2x5Right(0),
  RelE5x5(0),
  
  EtaWidth(0),
  PhiWidth(0),
  CoviEtaiPhi(0),
  CoviPhiiPhi(0),
  
  NVertexes(0),
  RelPreshowerEnergy(0),

  //MVA v2 and v1
  RelIsoEcal(0),
  RelIsoHcal(0),
  tIso1(0),
  tIso3(0),
  tIso2(0),
  ScEta(0.)
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void MVATools::InitializeMVA(int VariableType, TString EndcapWeights,TString BarrelWeights) {
  
  if (fReaderEndcap) delete fReaderEndcap;  
  if (fReaderBarrel) delete fReaderBarrel; 
  
  fReaderEndcap = new TMVA::Reader( "!Color:!Silent:Error" );       
  fReaderBarrel = new TMVA::Reader( "!Color:!Silent:Error" );   
  
  TMVA::Reader *readers[2];
  readers[0]  = fReaderEndcap;
  readers[1]  = fReaderBarrel;  
  
  for (UInt_t i=0; i<2; ++i) {

    if(VariableType==0||VariableType==1||VariableType==2){
      readers[i]->AddVariable( "HoE", &HoE );
      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
      readers[i]->AddVariable( "tIso1", &tIso1 );
      readers[i]->AddVariable( "tIso3", &tIso3 );
      readers[i]->AddVariable( "tIso2", &tIso2 );
      readers[i]->AddVariable( "R9", &R9 );
    }
    
    if(VariableType==3||VariableType==4){
      readers[i]->AddVariable( "HoE", &HoE );
      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
      readers[i]->AddVariable( "tIso1abs", &tIso1abs );
      readers[i]->AddVariable( "tIso3abs", &tIso3abs );
      readers[i]->AddVariable( "tIso2abs", &tIso2abs );
      readers[i]->AddVariable( "R9", &R9 );
    }
    
    if(VariableType==1||VariableType==2){
      readers[i]->AddVariable( "RelIsoEcal", &RelIsoEcal );
      readers[i]->AddVariable( "RelIsoHcal", &RelIsoHcal );
      readers[i]->AddVariable( "RelEMax", &RelEMax );
      readers[i]->AddVariable( "RelETop", &RelETop );
      readers[i]->AddVariable( "RelEBottom", &RelEBottom );
      readers[i]->AddVariable( "RelELeft", &RelELeft );
      readers[i]->AddVariable( "RelERight", &RelERight );
      readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
      readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
      readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
      readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
      readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
      readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
    }
    
    if(VariableType==3||VariableType==4){
      readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
      readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
      readers[i]->AddVariable( "RelEMax", &RelEMax );
      readers[i]->AddVariable( "RelETop", &RelEBottom);
      readers[i]->AddVariable( "RelEBottom", &RelEBottom );
      readers[i]->AddVariable( "RelELeft", &RelELeft );
      readers[i]->AddVariable( "RelERight", &RelERight );
      readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
      readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
      readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
      readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
      readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
      readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
    }
    
    if(VariableType==2||VariableType==3||VariableType==4){
      readers[i]->AddVariable( "EtaWidth", &EtaWidth );
      readers[i]->AddVariable( "PhiWidth", &PhiWidth );
      readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
      readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
      if(VariableType==4){
        readers[i]->AddVariable( "NVertexes", &NVertexes );
      }
      if(i==0){
        readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
      }
    }

    if(VariableType==6){
      readers[i]->AddVariable( "HoE", &HoE );
      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
      readers[i]->AddVariable( "tIso1abs", &tIso1abs );
      readers[i]->AddVariable( "tIso3abs", &tIso3abs );
      readers[i]->AddVariable( "tIso2abs", &tIso2abs );
      readers[i]->AddVariable( "R9", &R9 );
      readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
      readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
      readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
      readers[i]->AddVariable( "EtaWidth", &EtaWidth );
      readers[i]->AddVariable( "PhiWidth", &PhiWidth );
      readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
      readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
      readers[i]->AddVariable( "NVertexes", &NVertexes );
      if(i==0){
      readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
      }
    }
    
    if(VariableType==7){
      readers[i]->AddVariable( "HoE", &HoE );
      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
      readers[i]->AddVariable( "tIso1abs", &tIso1abs );
      readers[i]->AddVariable( "tIso3abs", &tIso3abs );
      readers[i]->AddVariable( "tIso2abs", &tIso2abs );
      readers[i]->AddVariable( "R9", &R9 );
      readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
      readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
      readers[i]->AddVariable( "NVertexes", &NVertexes );
      readers[i]->AddVariable( "ScEta", &ScEta );
    }    

    if(VariableType==10){
      readers[i]->AddVariable( "HoE", &HoE );
      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
      readers[i]->AddVariable( "tIso1abs", &tIso1abs );
      readers[i]->AddVariable( "tIso3abs", &tIso3abs );
      readers[i]->AddVariable( "tIso2abs", &tIso2abs );
      readers[i]->AddVariable( "R9", &R9 );
      readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
      readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
      readers[i]->AddVariable( "NVertexes", &NVertexes );
      readers[i]->AddVariable( "ScEta", &ScEta );
      readers[i]->AddVariable( "EtaWidth", &EtaWidth );
      readers[i]->AddVariable( "PhiWidth", &PhiWidth );      
    }    

    if(VariableType==1201){
      readers[i]->AddVariable( "myphoton_pfchargedisogood03", &myphoton_pfchargedisogood03);
      readers[i]->AddVariable( "myphoton_pfchargedisobad03", &myphoton_pfchargedisobad03); 
      readers[i]->AddVariable( "myphoton_pfphotoniso03", &myphoton_pfphotoniso03 );
      readers[i]->AddVariable( "myphoton_sieie", &myphoton_sieie ); 
      readers[i]->AddVariable( "myphoton_sieip", &myphoton_sieip );
      readers[i]->AddVariable( "myphoton_etawidth", &myphoton_etawidth ); 
      readers[i]->AddVariable( "myphoton_phiwidth", &myphoton_phiwidth );
      readers[i]->AddVariable( "myphoton_r9", &myphoton_r9 ); 
      readers[i]->AddVariable( "myphoton_s4ratio", &myphoton_s4ratio );
      readers[i]->AddVariable( "myphoton_SCeta", &myphoton_SCeta ); 
      readers[i]->AddVariable( "event_rho", &event_rho );
      if(i==0){
        readers[i]->AddVariable( "myphoton_ESEffSigmaRR", &myphoton_ESEffSigmaRR);
      } 
    }
    
  }
  
  fReaderEndcap->BookMVA("BDT method",EndcapWeights);
  fReaderBarrel->BookMVA("BDT method",BarrelWeights);
  
  assert(fReaderEndcap);
  assert(fReaderBarrel);
  
}

Bool_t MVATools::PassMVASelection(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,Float_t bdtCutBarrel, Float_t bdtCutEndcap, const ElectronCol* els, Bool_t applyElectronVeto) {
  
  //initilize the bool value
  PassMVA=kFALSE;
  
  Float_t photon_bdt =  MVATools::GetMVAbdtValue_2011(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
  
  if (isbarrel) {
    if(bdt>bdtCutBarrel){
      PassMVA=kTRUE;    
    }
  }
  else {
    if(bdt>bdtCutEndcap){
      PassMVA=kTRUE;    
    }
  } 
  return PassMVA;
}

//---------------------------------------------------------------------------------
Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
  
  // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
  float cic4_allcuts_temp_sublead[] = { 
    3.8,         2.2,         1.77,        1.29,
    11.7,        3.4,         3.9,         1.84,
    3.5,         2.2,         2.3,         1.45,
    0.0106,      0.0097,      0.028,       0.027,
    0.082,       0.062,       0.065,       0.048,
    0.94,        0.36,        0.94,        0.32,
    1.,          0.062,       0.97,        0.97,
    1.5,         1.5,         1.5,         1.5 };  // the last line is PixelmatchVeto and un-used
  
  //initilize the bool value
  PassElecVetoInt=0;
  
  dRTrack = PhotonTools::ElectronVetoCiC(p, els);
  
  ScEta_MVA=p->SCluster()->Eta();
  
  isbarrel = (fabs(ScEta_MVA)<1.4442);

  R9 = p->R9();
  //R9 = p->E33()/p->SCluster()->RawEnergy();
  
  // check which category it is ...
  _tCat = 1;
  if ( !isbarrel ) _tCat = 3;
  if ( R9 < 0.94 ) _tCat++;
  
  //Electron Veto
  if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
    PassElecVetoInt=1;
  }
  
  return  PassElecVetoInt;
  
}

//--------------------------------------------------------------------------------------------------

Float_t MVATools::GetMVAbdtValue_2012_globe(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho, const PFCandidateCol *fPFCands,const ElectronCol* els,Bool_t applyElectronVeto) {
  
  //get the variables used to compute MVA variables
  ecalIso3 = p->EcalRecHitIsoDr03();
  ecalIso4 = p->EcalRecHitIsoDr04();
  hcalIso4 = p->HcalTowerSumEtDr04();
  
  wVtxInd = 0;
  
  trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
    
  // track iso worst vtx
  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
  
  combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
  combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
  
  RawEnergy = p->SCluster()->RawEnergy();
  
  ScEta = p->SCluster()->Eta();
  
  //mva varialbes v1 and v2
  tIso1 = (combIso1) *50./p->Et();
  tIso3 = (trackIso1)*50./p->Et();
  tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
  RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
  RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();

  //compute mva variables for v3
  HoE = p->HadOverEm();
  covIEtaIEta = p->CoviEtaiEta();
  tIso1abs = combIso1;
  tIso3abs = trackIso1;
  tIso2abs = combIso2;
  R9 = p->R9();

  absIsoEcal=(ecalIso3-0.17*_tRho);
  absIsoHcal=(hcalIso4-0.17*_tRho);
  RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
  RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
  RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
  RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
  RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
  RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
  RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
  RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
  RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
  RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
  RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
  
  EtaWidth=p->EtaWidth();
  PhiWidth=p->PhiWidth();
  CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi(); 
  CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();

  RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
  NVertexes=vtxCol->GetEntries();
  
  //spectator variables
  Pt_MVA=p->Pt();
  ScEta_MVA=p->SCluster()->Eta();

  //

  isbarrel = (fabs(ScEta_MVA)<1.4442);
  
  //variable 1201
  myphoton_pfchargedisogood03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
  myphoton_pfchargedisobad03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol); 
  myphoton_pfphotoniso03=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
  myphoton_sieie=covIEtaIEta; 
  myphoton_sieip=CoviEtaiPhi;
  myphoton_etawidth=EtaWidth; 
  myphoton_phiwidth=PhiWidth;
  myphoton_r9=R9; 
  myphoton_s4ratio=p->SCluster()->Seed()->E2x2()/p->SCluster()->Seed()->E5x5();
  myphoton_SCeta=ScEta_MVA; 
  event_rho= _tRho;

  myphoton_ESEffSigmaRR=-99;

  if(!isbarrel){
    myphoton_ESEffSigmaRR=float(sqrt((p->SCluster()->PsEffWidthSigmaXX())*(p->SCluster()->PsEffWidthSigmaXX())+(p->SCluster()->PsEffWidthSigmaYY())*(p->SCluster()->PsEffWidthSigmaYY())));
  }
  
  if (isbarrel) {
    reader = fReaderBarrel;
  }
  else {
    reader = fReaderEndcap;
  }
  
  assert(reader); 

  bdt = reader->EvaluateMVA("BDT method");

  /* printf("HoE: %f\n",HoE);
  printf("covIEtaIEta: %f\n",covIEtaIEta);
  printf("tIso1abs: %f\n",tIso1abs);
  printf("tIso3abs: %f\n",tIso3abs);
  printf("tIso2abs: %f\n",tIso2abs);
  
  printf("absIsoEcal: %f\n",absIsoEcal);
  printf("absIsoHcal: %f\n",absIsoHcal);
  printf("RelEMax: %f\n",RelEMax);
  printf("RelETop: %f\n",RelETop);
  printf("RelEBottom: %f\n",RelEBottom);
  printf("RelELeft: %f\n",RelELeft);
  printf("RelERight: %f\n",RelERight);
  printf("RelE2x5Max: %f\n",RelE2x5Max);
  printf("RelE2x5Top: %f\n",RelE2x5Top);
  printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
  printf("RelE2x5Left: %f\n",RelE2x5Left);
  printf("RelE2x5Right;: %f\n",RelE2x5Right);
  printf("RelE5x5: %f\n",RelE5x5);
  
  printf("EtaWidth: %f\n",EtaWidth);
  printf("PhiWidth: %f\n",PhiWidth);
  printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
  printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
  
  if (!isbarrel) {
    printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
    }*/
  
  return bdt;
}

Float_t MVATools::GetMVAbdtValue_2011(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els,Bool_t applyElectronVeto) {
  
  //get the variables used to compute MVA variables
  ecalIso3 = p->EcalRecHitIsoDr03();
  ecalIso4 = p->EcalRecHitIsoDr04();
  hcalIso4 = p->HcalTowerSumEtDr04();
  
  wVtxInd = 0;
  
  trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
    
  // track iso worst vtx
  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
  
  combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
  combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
  
  RawEnergy = p->SCluster()->RawEnergy();
  
  ScEta = p->SCluster()->Eta();
  
  //mva varialbes v1 and v2
  tIso1 = (combIso1) *50./p->Et();
  tIso3 = (trackIso1)*50./p->Et();
  tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
  RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
  RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();

  //compute mva variables for v3
  HoE = p->HadOverEm();
  covIEtaIEta = p->CoviEtaiEta();
  tIso1abs = combIso1;
  tIso3abs = trackIso1;
  tIso2abs = combIso2;
  R9 = p->R9();

  absIsoEcal=(ecalIso3-0.17*_tRho);
  absIsoHcal=(hcalIso4-0.17*_tRho);
  RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
  RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
  RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
  RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
  RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
  RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
  RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
  RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
  RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
  RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
  RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
  
  EtaWidth=p->EtaWidth();
  PhiWidth=p->PhiWidth();
  CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi(); 
  CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();

  RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
  NVertexes=vtxCol->GetEntries();
  
  //spectator variables
  Pt_MVA=p->Pt();
  ScEta_MVA=p->SCluster()->Eta();

  //

  isbarrel = (fabs(ScEta_MVA)<1.4442);
  
  if (isbarrel) {
    reader = fReaderBarrel;
  }
  else {
    reader = fReaderEndcap;
  }
  
  assert(reader); 

  bdt = reader->EvaluateMVA("BDT method");

  /* printf("HoE: %f\n",HoE);
  printf("covIEtaIEta: %f\n",covIEtaIEta);
  printf("tIso1abs: %f\n",tIso1abs);
  printf("tIso3abs: %f\n",tIso3abs);
  printf("tIso2abs: %f\n",tIso2abs);
  
  printf("absIsoEcal: %f\n",absIsoEcal);
  printf("absIsoHcal: %f\n",absIsoHcal);
  printf("RelEMax: %f\n",RelEMax);
  printf("RelETop: %f\n",RelETop);
  printf("RelEBottom: %f\n",RelEBottom);
  printf("RelELeft: %f\n",RelELeft);
  printf("RelERight: %f\n",RelERight);
  printf("RelE2x5Max: %f\n",RelE2x5Max);
  printf("RelE2x5Top: %f\n",RelE2x5Top);
  printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
  printf("RelE2x5Left: %f\n",RelE2x5Left);
  printf("RelE2x5Right;: %f\n",RelE2x5Right);
  printf("RelE5x5: %f\n",RelE5x5);
  
  printf("EtaWidth: %f\n",EtaWidth);
  printf("PhiWidth: %f\n",PhiWidth);
  printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
  printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
  
  if (!isbarrel) {
    printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
    }*/
  
  return bdt;
}
Revision:	1.14
Committed:	Sun Jun 10 16:28:40 2012 UTC (12 years, 10 months ago) by mingyang
Content type:	text/plain
Branch:	MAIN
Changes since 1.13:	+2 -12 lines
Log Message:	* empty log message *
#	Content
1	// $Id: MVATools.cc,v 1.13 2012/05/29 21:20:12 mingyang Exp $
2
3	#include "MitPhysics/Utils/interface/PhotonTools.h"
4	#include "MitPhysics/Utils/interface/MVATools.h"
5	#include "MitPhysics/Utils/interface/ElectronTools.h"
6	#include "MitPhysics/Utils/interface/IsolationTools.h"
7	#include "MitAna/DataTree/interface/PFCandidateCol.h"
8	#include "MitAna/DataTree/interface/StableData.h"
9	#include <TMath.h>
10	#include <TFile.h>
11	#include <TRandom3.h>
12	#include "TMVA/Tools.h"//MVA
13	#include "TMVA/Reader.h"//MVA
14
15
16	ClassImp(mithep::MVATools)
17
18	using namespace mithep;
19
20
21	//--------------------------------------------------------------------------------------------------
22	MVATools::MVATools():
23	fReaderEndcap(0),
24	fReaderBarrel(0),
25
26	//MVA Variables v4
27	HoE(0),
28	covIEtaIEta(0),
29	tIso1abs(0),
30	tIso3abs(0),
31	tIso2abs(0),
32	R9(0),
33
34	absIsoEcal(0),
35	absIsoHcal(0),
36	RelEMax(0),
37	RelETop(0),
38	RelEBottom(0),
39	RelELeft(0),
40	RelERight(0),
41	RelE2x5Max(0),
42	RelE2x5Top(0),
43	RelE2x5Bottom(0),
44	RelE2x5Left(0),
45	RelE2x5Right(0),
46	RelE5x5(0),
47
48	EtaWidth(0),
49	PhiWidth(0),
50	CoviEtaiPhi(0),
51	CoviPhiiPhi(0),
52
53	NVertexes(0),
54	RelPreshowerEnergy(0),
55
56	//MVA v2 and v1
57	RelIsoEcal(0),
58	RelIsoHcal(0),
59	tIso1(0),
60	tIso3(0),
61	tIso2(0),
62	ScEta(0.)
63	{
64	// Constructor.
65	}
66
67	//--------------------------------------------------------------------------------------------------
68	void MVATools::InitializeMVA(int VariableType, TString EndcapWeights,TString BarrelWeights) {
69
70	if (fReaderEndcap) delete fReaderEndcap;
71	if (fReaderBarrel) delete fReaderBarrel;
72
73	fReaderEndcap = new TMVA::Reader( "!Color:!Silent:Error" );
74	fReaderBarrel = new TMVA::Reader( "!Color:!Silent:Error" );
75
76	TMVA::Reader *readers[2];
77	readers[0] = fReaderEndcap;
78	readers[1] = fReaderBarrel;
79
80	for (UInt_t i=0; i<2; ++i) {
81
82	if(VariableType==0\|\|VariableType==1\|\|VariableType==2){
83	readers[i]->AddVariable( "HoE", &HoE );
84	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
85	readers[i]->AddVariable( "tIso1", &tIso1 );
86	readers[i]->AddVariable( "tIso3", &tIso3 );
87	readers[i]->AddVariable( "tIso2", &tIso2 );
88	readers[i]->AddVariable( "R9", &R9 );
89	}
90
91	if(VariableType==3\|\|VariableType==4){
92	readers[i]->AddVariable( "HoE", &HoE );
93	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
94	readers[i]->AddVariable( "tIso1abs", &tIso1abs );
95	readers[i]->AddVariable( "tIso3abs", &tIso3abs );
96	readers[i]->AddVariable( "tIso2abs", &tIso2abs );
97	readers[i]->AddVariable( "R9", &R9 );
98	}
99
100	if(VariableType==1\|\|VariableType==2){
101	readers[i]->AddVariable( "RelIsoEcal", &RelIsoEcal );
102	readers[i]->AddVariable( "RelIsoHcal", &RelIsoHcal );
103	readers[i]->AddVariable( "RelEMax", &RelEMax );
104	readers[i]->AddVariable( "RelETop", &RelETop );
105	readers[i]->AddVariable( "RelEBottom", &RelEBottom );
106	readers[i]->AddVariable( "RelELeft", &RelELeft );
107	readers[i]->AddVariable( "RelERight", &RelERight );
108	readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
109	readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
110	readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
111	readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
112	readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
113	readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
114	}
115
116	if(VariableType==3\|\|VariableType==4){
117	readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
118	readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
119	readers[i]->AddVariable( "RelEMax", &RelEMax );
120	readers[i]->AddVariable( "RelETop", &RelEBottom);
121	readers[i]->AddVariable( "RelEBottom", &RelEBottom );
122	readers[i]->AddVariable( "RelELeft", &RelELeft );
123	readers[i]->AddVariable( "RelERight", &RelERight );
124	readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
125	readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
126	readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
127	readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
128	readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
129	readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
130	}
131
132	if(VariableType==2\|\|VariableType==3\|\|VariableType==4){
133	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
134	readers[i]->AddVariable( "PhiWidth", &PhiWidth );
135	readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
136	readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
137	if(VariableType==4){
138	readers[i]->AddVariable( "NVertexes", &NVertexes );
139	}
140	if(i==0){
141	readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
142	}
143	}
144
145	if(VariableType==6){
146	readers[i]->AddVariable( "HoE", &HoE );
147	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
148	readers[i]->AddVariable( "tIso1abs", &tIso1abs );
149	readers[i]->AddVariable( "tIso3abs", &tIso3abs );
150	readers[i]->AddVariable( "tIso2abs", &tIso2abs );
151	readers[i]->AddVariable( "R9", &R9 );
152	readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
153	readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
154	readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
155	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
156	readers[i]->AddVariable( "PhiWidth", &PhiWidth );
157	readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
158	readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
159	readers[i]->AddVariable( "NVertexes", &NVertexes );
160	if(i==0){
161	readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
162	}
163	}
164
165	if(VariableType==7){
166	readers[i]->AddVariable( "HoE", &HoE );
167	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
168	readers[i]->AddVariable( "tIso1abs", &tIso1abs );
169	readers[i]->AddVariable( "tIso3abs", &tIso3abs );
170	readers[i]->AddVariable( "tIso2abs", &tIso2abs );
171	readers[i]->AddVariable( "R9", &R9 );
172	readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
173	readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
174	readers[i]->AddVariable( "NVertexes", &NVertexes );
175	readers[i]->AddVariable( "ScEta", &ScEta );
176	}
177
178	if(VariableType==10){
179	readers[i]->AddVariable( "HoE", &HoE );
180	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
181	readers[i]->AddVariable( "tIso1abs", &tIso1abs );
182	readers[i]->AddVariable( "tIso3abs", &tIso3abs );
183	readers[i]->AddVariable( "tIso2abs", &tIso2abs );
184	readers[i]->AddVariable( "R9", &R9 );
185	readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
186	readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
187	readers[i]->AddVariable( "NVertexes", &NVertexes );
188	readers[i]->AddVariable( "ScEta", &ScEta );
189	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
190	readers[i]->AddVariable( "PhiWidth", &PhiWidth );
191	}
192
193	if(VariableType==1201){
194	readers[i]->AddVariable( "myphoton_pfchargedisogood03", &myphoton_pfchargedisogood03);
195	readers[i]->AddVariable( "myphoton_pfchargedisobad03", &myphoton_pfchargedisobad03);
196	readers[i]->AddVariable( "myphoton_pfphotoniso03", &myphoton_pfphotoniso03 );
197	readers[i]->AddVariable( "myphoton_sieie", &myphoton_sieie );
198	readers[i]->AddVariable( "myphoton_sieip", &myphoton_sieip );
199	readers[i]->AddVariable( "myphoton_etawidth", &myphoton_etawidth );
200	readers[i]->AddVariable( "myphoton_phiwidth", &myphoton_phiwidth );
201	readers[i]->AddVariable( "myphoton_r9", &myphoton_r9 );
202	readers[i]->AddVariable( "myphoton_s4ratio", &myphoton_s4ratio );
203	readers[i]->AddVariable( "myphoton_SCeta", &myphoton_SCeta );
204	readers[i]->AddVariable( "event_rho", &event_rho );
205	if(i==0){
206	readers[i]->AddVariable( "myphoton_ESEffSigmaRR", &myphoton_ESEffSigmaRR);
207	}
208	}
209
210	}
211
212	fReaderEndcap->BookMVA("BDT method",EndcapWeights);
213	fReaderBarrel->BookMVA("BDT method",BarrelWeights);
214
215	assert(fReaderEndcap);
216	assert(fReaderBarrel);
217
218	}
219
220	Bool_t MVATools::PassMVASelection(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,Float_t bdtCutBarrel, Float_t bdtCutEndcap, const ElectronCol* els, Bool_t applyElectronVeto) {
221
222	//initilize the bool value
223	PassMVA=kFALSE;
224
225	Float_t photon_bdt = MVATools::GetMVAbdtValue_2011(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
226
227	if (isbarrel) {
228	if(bdt>bdtCutBarrel){
229	PassMVA=kTRUE;
230	}
231	}
232	else {
233	if(bdt>bdtCutEndcap){
234	PassMVA=kTRUE;
235	}
236	}
237	return PassMVA;
238	}
239
240	//---------------------------------------------------------------------------------
241	Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
242
243	// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
244	float cic4_allcuts_temp_sublead[] = {
245	3.8, 2.2, 1.77, 1.29,
246	11.7, 3.4, 3.9, 1.84,
247	3.5, 2.2, 2.3, 1.45,
248	0.0106, 0.0097, 0.028, 0.027,
249	0.082, 0.062, 0.065, 0.048,
250	0.94, 0.36, 0.94, 0.32,
251	1., 0.062, 0.97, 0.97,
252	1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
253
254	//initilize the bool value
255	PassElecVetoInt=0;
256
257	dRTrack = PhotonTools::ElectronVetoCiC(p, els);
258
259	ScEta_MVA=p->SCluster()->Eta();
260
261	isbarrel = (fabs(ScEta_MVA)<1.4442);
262
263	R9 = p->R9();
264	//R9 = p->E33()/p->SCluster()->RawEnergy();
265
266	// check which category it is ...
267	_tCat = 1;
268	if ( !isbarrel ) _tCat = 3;
269	if ( R9 < 0.94 ) _tCat++;
270
271	//Electron Veto
272	if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
273	PassElecVetoInt=1;
274	}
275
276	return PassElecVetoInt;
277
278	}
279
280	//--------------------------------------------------------------------------------------------------
281
282	Float_t MVATools::GetMVAbdtValue_2012_globe(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho, const PFCandidateCol fPFCands,const ElectronCol els,Bool_t applyElectronVeto) {
283
284	//get the variables used to compute MVA variables
285	ecalIso3 = p->EcalRecHitIsoDr03();
286	ecalIso4 = p->EcalRecHitIsoDr04();
287	hcalIso4 = p->HcalTowerSumEtDr04();
288
289	wVtxInd = 0;
290
291	trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
292
293	// track iso worst vtx
294	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
295
296	combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
297	combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
298
299	RawEnergy = p->SCluster()->RawEnergy();
300
301	ScEta = p->SCluster()->Eta();
302
303	//mva varialbes v1 and v2
304	tIso1 = (combIso1) *50./p->Et();
305	tIso3 = (trackIso1)*50./p->Et();
306	tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
307	RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
308	RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
309
310	//compute mva variables for v3
311	HoE = p->HadOverEm();
312	covIEtaIEta = p->CoviEtaiEta();
313	tIso1abs = combIso1;
314	tIso3abs = trackIso1;
315	tIso2abs = combIso2;
316	R9 = p->R9();
317
318	absIsoEcal=(ecalIso3-0.17*_tRho);
319	absIsoHcal=(hcalIso4-0.17*_tRho);
320	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
321	RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
322	RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
323	RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
324	RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
325	RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
326	RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
327	RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
328	RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
329	RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
330	RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
331
332	EtaWidth=p->EtaWidth();
333	PhiWidth=p->PhiWidth();
334	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
335	CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
336
337	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
338	NVertexes=vtxCol->GetEntries();
339
340	//spectator variables
341	Pt_MVA=p->Pt();
342	ScEta_MVA=p->SCluster()->Eta();
343
344	//
345
346	isbarrel = (fabs(ScEta_MVA)<1.4442);
347
348	//variable 1201
349	myphoton_pfchargedisogood03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
350	myphoton_pfchargedisobad03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol);
351	myphoton_pfphotoniso03=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
352	myphoton_sieie=covIEtaIEta;
353	myphoton_sieip=CoviEtaiPhi;
354	myphoton_etawidth=EtaWidth;
355	myphoton_phiwidth=PhiWidth;
356	myphoton_r9=R9;
357	myphoton_s4ratio=p->SCluster()->Seed()->E2x2()/p->SCluster()->Seed()->E5x5();
358	myphoton_SCeta=ScEta_MVA;
359	event_rho= _tRho;
360
361	myphoton_ESEffSigmaRR=-99;
362
363	if(!isbarrel){
364	myphoton_ESEffSigmaRR=float(sqrt((p->SCluster()->PsEffWidthSigmaXX())(p->SCluster()->PsEffWidthSigmaXX())+(p->SCluster()->PsEffWidthSigmaYY())(p->SCluster()->PsEffWidthSigmaYY())));
365	}
366
367	if (isbarrel) {
368	reader = fReaderBarrel;
369	}
370	else {
371	reader = fReaderEndcap;
372	}
373
374	assert(reader);
375
376	bdt = reader->EvaluateMVA("BDT method");
377
378	/* printf("HoE: %f\n",HoE);
379	printf("covIEtaIEta: %f\n",covIEtaIEta);
380	printf("tIso1abs: %f\n",tIso1abs);
381	printf("tIso3abs: %f\n",tIso3abs);
382	printf("tIso2abs: %f\n",tIso2abs);
383
384	printf("absIsoEcal: %f\n",absIsoEcal);
385	printf("absIsoHcal: %f\n",absIsoHcal);
386	printf("RelEMax: %f\n",RelEMax);
387	printf("RelETop: %f\n",RelETop);
388	printf("RelEBottom: %f\n",RelEBottom);
389	printf("RelELeft: %f\n",RelELeft);
390	printf("RelERight: %f\n",RelERight);
391	printf("RelE2x5Max: %f\n",RelE2x5Max);
392	printf("RelE2x5Top: %f\n",RelE2x5Top);
393	printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
394	printf("RelE2x5Left: %f\n",RelE2x5Left);
395	printf("RelE2x5Right;: %f\n",RelE2x5Right);
396	printf("RelE5x5: %f\n",RelE5x5);
397
398	printf("EtaWidth: %f\n",EtaWidth);
399	printf("PhiWidth: %f\n",PhiWidth);
400	printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
401	printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
402
403	if (!isbarrel) {
404	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
405	}*/
406
407	return bdt;
408	}
409
410	Float_t MVATools::GetMVAbdtValue_2011(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els,Bool_t applyElectronVeto) {
411
412	//get the variables used to compute MVA variables
413	ecalIso3 = p->EcalRecHitIsoDr03();
414	ecalIso4 = p->EcalRecHitIsoDr04();
415	hcalIso4 = p->HcalTowerSumEtDr04();
416
417	wVtxInd = 0;
418
419	trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
420
421	// track iso worst vtx
422	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
423
424	combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
425	combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
426
427	RawEnergy = p->SCluster()->RawEnergy();
428
429	ScEta = p->SCluster()->Eta();
430
431	//mva varialbes v1 and v2
432	tIso1 = (combIso1) *50./p->Et();
433	tIso3 = (trackIso1)*50./p->Et();
434	tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
435	RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
436	RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
437
438	//compute mva variables for v3
439	HoE = p->HadOverEm();
440	covIEtaIEta = p->CoviEtaiEta();
441	tIso1abs = combIso1;
442	tIso3abs = trackIso1;
443	tIso2abs = combIso2;
444	R9 = p->R9();
445
446	absIsoEcal=(ecalIso3-0.17*_tRho);
447	absIsoHcal=(hcalIso4-0.17*_tRho);
448	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
449	RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
450	RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
451	RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
452	RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
453	RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
454	RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
455	RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
456	RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
457	RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
458	RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
459
460	EtaWidth=p->EtaWidth();
461	PhiWidth=p->PhiWidth();
462	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
463	CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
464
465	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
466	NVertexes=vtxCol->GetEntries();
467
468	//spectator variables
469	Pt_MVA=p->Pt();
470	ScEta_MVA=p->SCluster()->Eta();
471
472	//
473
474	isbarrel = (fabs(ScEta_MVA)<1.4442);
475
476	if (isbarrel) {
477	reader = fReaderBarrel;
478	}
479	else {
480	reader = fReaderEndcap;
481	}
482
483	assert(reader);
484
485	bdt = reader->EvaluateMVA("BDT method");
486
487	/* printf("HoE: %f\n",HoE);
488	printf("covIEtaIEta: %f\n",covIEtaIEta);
489	printf("tIso1abs: %f\n",tIso1abs);
490	printf("tIso3abs: %f\n",tIso3abs);
491	printf("tIso2abs: %f\n",tIso2abs);
492
493	printf("absIsoEcal: %f\n",absIsoEcal);
494	printf("absIsoHcal: %f\n",absIsoHcal);
495	printf("RelEMax: %f\n",RelEMax);
496	printf("RelETop: %f\n",RelETop);
497	printf("RelEBottom: %f\n",RelEBottom);
498	printf("RelELeft: %f\n",RelELeft);
499	printf("RelERight: %f\n",RelERight);
500	printf("RelE2x5Max: %f\n",RelE2x5Max);
501	printf("RelE2x5Top: %f\n",RelE2x5Top);
502	printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
503	printf("RelE2x5Left: %f\n",RelE2x5Left);
504	printf("RelE2x5Right;: %f\n",RelE2x5Right);
505	printf("RelE5x5: %f\n",RelE5x5);
506
507	printf("EtaWidth: %f\n",EtaWidth);
508	printf("PhiWidth: %f\n",PhiWidth);
509	printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
510	printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
511
512	if (!isbarrel) {
513	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
514	}*/
515
516	return bdt;
517	}