Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.10 2011/12/19 23:45:00 bendavid 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::GetMVAbdtValueOld(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(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::GetMVAbdtValueOld(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.11
Committed:	Sun May 27 16:39:19 2012 UTC (12 years, 11 months ago) by mingyang
Content type:	text/plain
Branch:	MAIN
Changes since 1.10:	+153 -3 lines
Log Message:	2012 id added
#	User	Rev	Content
1	mingyang	1.11	// $Id: MVATools.cc,v 1.10 2011/12/19 23:45:00 bendavid Exp $
2	mingyang	1.1
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	mingyang	1.11	#include "MitAna/DataTree/interface/PFCandidateCol.h"
8	mingyang	1.1	#include "MitAna/DataTree/interface/StableData.h"
9	mingyang	1.11	#include <TMath.h>
10	mingyang	1.1	#include <TFile.h>
11			#include <TRandom3.h>
12			#include "TMVA/Tools.h"//MVA
13			#include "TMVA/Reader.h"//MVA
14
15	mingyang	1.11
16	mingyang	1.1	ClassImp(mithep::MVATools)
17
18			using namespace mithep;
19
20
21			//--------------------------------------------------------------------------------------------------
22			MVATools::MVATools():
23			fReaderEndcap(0),
24			fReaderBarrel(0),
25
26	mingyang	1.4	//MVA Variables v4
27	mingyang	1.1	HoE(0),
28			covIEtaIEta(0),
29	mingyang	1.4	tIso1abs(0),
30			tIso3abs(0),
31			tIso2abs(0),
32	mingyang	1.1	R9(0),
33
34	mingyang	1.4	absIsoEcal(0),
35			absIsoHcal(0),
36	mingyang	1.1	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	mingyang	1.4	CoviEtaiPhi(0),
51	mingyang	1.1	CoviPhiiPhi(0),
52
53	mingyang	1.4	NVertexes(0),
54			RelPreshowerEnergy(0),
55
56			//MVA v2 and v1
57			RelIsoEcal(0),
58			RelIsoHcal(0),
59			tIso1(0),
60			tIso3(0),
61	bendavid	1.9	tIso2(0),
62			ScEta(0.)
63	mingyang	1.1	{
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	mingyang	1.4
82	mingyang	1.1	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	mingyang	1.4	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	mingyang	1.1	if(VariableType==1\|\|VariableType==2){
101	mingyang	1.4	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	mingyang	1.1	readers[i]->AddVariable( "RelELeft", &RelELeft );
123	mingyang	1.4	readers[i]->AddVariable( "RelERight", &RelERight );
124	mingyang	1.1	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	mingyang	1.4	if(VariableType==2\|\|VariableType==3\|\|VariableType==4){
133	mingyang	1.1	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
134			readers[i]->AddVariable( "PhiWidth", &PhiWidth );
135			readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
136			readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
137	mingyang	1.4	if(VariableType==4){
138			readers[i]->AddVariable( "NVertexes", &NVertexes );
139			}
140	mingyang	1.1	if(i==0){
141			readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
142			}
143			}
144	mingyang	1.6
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	bendavid	1.9
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	mingyang	1.6
178	bendavid	1.10	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	mingyang	1.11	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	mingyang	1.1	}
211
212			fReaderEndcap->BookMVA("BDT method",EndcapWeights);
213			fReaderBarrel->BookMVA("BDT method",BarrelWeights);
214	mingyang	1.4
215	mingyang	1.1	assert(fReaderEndcap);
216			assert(fReaderBarrel);
217
218			}
219
220	bendavid	1.7	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	mingyang	1.1
222			//initilize the bool value
223			PassMVA=kFALSE;
224
225	mingyang	1.11	Float_t photon_bdt = MVATools::GetMVAbdtValueOld(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
226	mingyang	1.1
227	mingyang	1.4	if (isbarrel) {
228			if(bdt>bdtCutBarrel){
229			PassMVA=kTRUE;
230			}
231	mingyang	1.1	}
232	mingyang	1.4	else {
233			if(bdt>bdtCutEndcap){
234			PassMVA=kTRUE;
235	mingyang	1.1	}
236	mingyang	1.4	}
237	mingyang	1.1	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	bendavid	1.9	R9 = p->R9();
264			//R9 = p->E33()/p->SCluster()->RawEnergy();
265	mingyang	1.1
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	mingyang	1.11	Float_t MVATools::GetMVAbdtValue(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	mingyang	1.1
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	bendavid	1.7	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	bendavid	1.8
293	mingyang	1.1	// track iso worst vtx
294	bendavid	1.7	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
295	mingyang	1.1
296			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
297			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
298
299			RawEnergy = p->SCluster()->RawEnergy();
300
301	bendavid	1.9	ScEta = p->SCluster()->Eta();
302
303	mingyang	1.4	//mva varialbes v1 and v2
304	mingyang	1.1	tIso1 = (combIso1) *50./p->Et();
305	bendavid	1.8	tIso3 = (trackIso1)*50./p->Et();
306	mingyang	1.1	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	mingyang	1.4
310			//compute mva variables for v3
311			HoE = p->HadOverEm();
312			covIEtaIEta = p->CoviEtaiEta();
313			tIso1abs = combIso1;
314	bendavid	1.8	tIso3abs = trackIso1;
315	mingyang	1.4	tIso2abs = combIso2;
316			R9 = p->R9();
317
318			absIsoEcal=(ecalIso3-0.17*_tRho);
319			absIsoHcal=(hcalIso4-0.17*_tRho);
320	mingyang	1.1	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	bendavid	1.10	EtaWidth=p->EtaWidth();
333			PhiWidth=p->PhiWidth();
334	mingyang	1.1	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
335			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
336	mingyang	1.4
337	mingyang	1.1	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
338	mingyang	1.4	NVertexes=vtxCol->GetEntries();
339	mingyang	1.1
340			//spectator variables
341			Pt_MVA=p->Pt();
342			ScEta_MVA=p->SCluster()->Eta();
343	mingyang	1.11
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	mingyang	1.1
398	mingyang	1.11	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::GetMVAbdtValueOld(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	mingyang	1.1	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	mingyang	1.5	/* printf("HoE: %f\n",HoE);
488	mingyang	1.1	printf("covIEtaIEta: %f\n",covIEtaIEta);
489	mingyang	1.4	printf("tIso1abs: %f\n",tIso1abs);
490			printf("tIso3abs: %f\n",tIso3abs);
491			printf("tIso2abs: %f\n",tIso2abs);
492	mingyang	1.1
493	mingyang	1.4	printf("absIsoEcal: %f\n",absIsoEcal);
494			printf("absIsoHcal: %f\n",absIsoHcal);
495	mingyang	1.1	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	mingyang	1.4	if (!isbarrel) {
513	mingyang	1.1	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
514	mingyang	1.5	}*/
515	mingyang	1.1
516			return bdt;
517			}