Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.12 2012/05/27 19:41:00 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,Bool_t applyNewScale,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(applyNewScale){
    if(isbarrel){
      myphoton_s4ratio=1.0055*myphoton_s4ratio;
    }
    if(!isbarrel){
      myphoton_s4ratio=1.0085*myphoton_s4ratio;
      myphoton_ESEffSigmaRR=1.04*myphoton_ESEffSigmaRR;
    }
  }
  
  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.13
Committed:	Tue May 29 21:20:12 2012 UTC (12 years, 11 months ago) by mingyang
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_028a
Changes since 1.12:	+13 -3 lines
Log Message:	add new scale
#	Content
1	// $Id: MVATools.cc,v 1.12 2012/05/27 19:41:00 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,Bool_t applyNewScale,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(applyNewScale){
368	if(isbarrel){
369	myphoton_s4ratio=1.0055*myphoton_s4ratio;
370	}
371	if(!isbarrel){
372	myphoton_s4ratio=1.0085*myphoton_s4ratio;
373	myphoton_ESEffSigmaRR=1.04*myphoton_ESEffSigmaRR;
374	}
375	}
376
377	if (isbarrel) {
378	reader = fReaderBarrel;
379	}
380	else {
381	reader = fReaderEndcap;
382	}
383
384	assert(reader);
385
386	bdt = reader->EvaluateMVA("BDT method");
387
388	/* printf("HoE: %f\n",HoE);
389	printf("covIEtaIEta: %f\n",covIEtaIEta);
390	printf("tIso1abs: %f\n",tIso1abs);
391	printf("tIso3abs: %f\n",tIso3abs);
392	printf("tIso2abs: %f\n",tIso2abs);
393
394	printf("absIsoEcal: %f\n",absIsoEcal);
395	printf("absIsoHcal: %f\n",absIsoHcal);
396	printf("RelEMax: %f\n",RelEMax);
397	printf("RelETop: %f\n",RelETop);
398	printf("RelEBottom: %f\n",RelEBottom);
399	printf("RelELeft: %f\n",RelELeft);
400	printf("RelERight: %f\n",RelERight);
401	printf("RelE2x5Max: %f\n",RelE2x5Max);
402	printf("RelE2x5Top: %f\n",RelE2x5Top);
403	printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
404	printf("RelE2x5Left: %f\n",RelE2x5Left);
405	printf("RelE2x5Right;: %f\n",RelE2x5Right);
406	printf("RelE5x5: %f\n",RelE5x5);
407
408	printf("EtaWidth: %f\n",EtaWidth);
409	printf("PhiWidth: %f\n",PhiWidth);
410	printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
411	printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
412
413	if (!isbarrel) {
414	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
415	}*/
416
417	return bdt;
418	}
419
420	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) {
421
422	//get the variables used to compute MVA variables
423	ecalIso3 = p->EcalRecHitIsoDr03();
424	ecalIso4 = p->EcalRecHitIsoDr04();
425	hcalIso4 = p->HcalTowerSumEtDr04();
426
427	wVtxInd = 0;
428
429	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?
430
431	// track iso worst vtx
432	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
433
434	combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
435	combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
436
437	RawEnergy = p->SCluster()->RawEnergy();
438
439	ScEta = p->SCluster()->Eta();
440
441	//mva varialbes v1 and v2
442	tIso1 = (combIso1) *50./p->Et();
443	tIso3 = (trackIso1)*50./p->Et();
444	tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
445	RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
446	RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
447
448	//compute mva variables for v3
449	HoE = p->HadOverEm();
450	covIEtaIEta = p->CoviEtaiEta();
451	tIso1abs = combIso1;
452	tIso3abs = trackIso1;
453	tIso2abs = combIso2;
454	R9 = p->R9();
455
456	absIsoEcal=(ecalIso3-0.17*_tRho);
457	absIsoHcal=(hcalIso4-0.17*_tRho);
458	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
459	RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
460	RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
461	RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
462	RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
463	RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
464	RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
465	RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
466	RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
467	RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
468	RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
469
470	EtaWidth=p->EtaWidth();
471	PhiWidth=p->PhiWidth();
472	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
473	CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
474
475	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
476	NVertexes=vtxCol->GetEntries();
477
478	//spectator variables
479	Pt_MVA=p->Pt();
480	ScEta_MVA=p->SCluster()->Eta();
481
482	//
483
484	isbarrel = (fabs(ScEta_MVA)<1.4442);
485
486	if (isbarrel) {
487	reader = fReaderBarrel;
488	}
489	else {
490	reader = fReaderEndcap;
491	}
492
493	assert(reader);
494
495	bdt = reader->EvaluateMVA("BDT method");
496
497	/* printf("HoE: %f\n",HoE);
498	printf("covIEtaIEta: %f\n",covIEtaIEta);
499	printf("tIso1abs: %f\n",tIso1abs);
500	printf("tIso3abs: %f\n",tIso3abs);
501	printf("tIso2abs: %f\n",tIso2abs);
502
503	printf("absIsoEcal: %f\n",absIsoEcal);
504	printf("absIsoHcal: %f\n",absIsoHcal);
505	printf("RelEMax: %f\n",RelEMax);
506	printf("RelETop: %f\n",RelETop);
507	printf("RelEBottom: %f\n",RelEBottom);
508	printf("RelELeft: %f\n",RelELeft);
509	printf("RelERight: %f\n",RelERight);
510	printf("RelE2x5Max: %f\n",RelE2x5Max);
511	printf("RelE2x5Top: %f\n",RelE2x5Top);
512	printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
513	printf("RelE2x5Left: %f\n",RelE2x5Left);
514	printf("RelE2x5Right;: %f\n",RelE2x5Right);
515	printf("RelE5x5: %f\n",RelE5x5);
516
517	printf("EtaWidth: %f\n",EtaWidth);
518	printf("PhiWidth: %f\n",PhiWidth);
519	printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
520	printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
521
522	if (!isbarrel) {
523	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
524	}*/
525
526	return bdt;
527	}