Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.14 2012/06/10 16:28:40 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);
        }*/
      readers[i]->AddVariable( "ph.r9", &myphoton_r9 ); 
      readers[i]->AddVariable( "ph.sigietaieta", &myphoton_sieie ); 
      readers[i]->AddVariable( "ph.scetawidth", &myphoton_etawidth ); 
      readers[i]->AddVariable( "ph.scphiwidth", &myphoton_phiwidth );
      readers[i]->AddVariable( "ph.idmva_CoviEtaiPhi", &myphoton_sieip );
      readers[i]->AddVariable( "ph.idmva_s4ratio", &myphoton_s4ratio );
      readers[i]->AddVariable( "ph.idmva_GammaIso", &myphoton_pfphotoniso03 );
      readers[i]->AddVariable( "ph.idmva_ChargedIso_selvtx", &myphoton_pfchargedisogood03);
      readers[i]->AddVariable( "ph.idmva_ChargedIso_worstvtx", &myphoton_pfchargedisobad03); 
      readers[i]->AddVariable( "ph.sceta", &myphoton_SCeta ); 
      readers[i]->AddVariable( "rho", &event_rho );
      if(i==0){
        readers[i]->AddVariable( "1.00023*ph.idmva_PsEffWidthSigmaRR + 0.0913", &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->S4Ratio();
  myphoton_SCeta=ScEta_MVA; 
  event_rho= _tRho;

  myphoton_ESEffSigmaRR=-99;

  if(!isbarrel){
    myphoton_ESEffSigmaRR=p->EffSigmaRR();
  }
  
  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.15
Committed:	Mon Jun 11 09:04:33 2012 UTC (12 years, 10 months ago) by mingyang
Content type:	text/plain
Branch:	MAIN
Changes since 1.14:	+19 -5 lines
Log Message:	newly trained mit mva for ichep added
#	User	Rev	Content
1	mingyang	1.15	// $Id: MVATools.cc,v 1.14 2012/06/10 16:28:40 mingyang 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	mingyang	1.15	/*readers[i]->AddVariable( "myphoton_pfchargedisogood03", &myphoton_pfchargedisogood03);
195	mingyang	1.11	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	mingyang	1.15	}*/
208			readers[i]->AddVariable( "ph.r9", &myphoton_r9 );
209			readers[i]->AddVariable( "ph.sigietaieta", &myphoton_sieie );
210			readers[i]->AddVariable( "ph.scetawidth", &myphoton_etawidth );
211			readers[i]->AddVariable( "ph.scphiwidth", &myphoton_phiwidth );
212			readers[i]->AddVariable( "ph.idmva_CoviEtaiPhi", &myphoton_sieip );
213			readers[i]->AddVariable( "ph.idmva_s4ratio", &myphoton_s4ratio );
214			readers[i]->AddVariable( "ph.idmva_GammaIso", &myphoton_pfphotoniso03 );
215			readers[i]->AddVariable( "ph.idmva_ChargedIso_selvtx", &myphoton_pfchargedisogood03);
216			readers[i]->AddVariable( "ph.idmva_ChargedIso_worstvtx", &myphoton_pfchargedisobad03);
217			readers[i]->AddVariable( "ph.sceta", &myphoton_SCeta );
218			readers[i]->AddVariable( "rho", &event_rho );
219			if(i==0){
220			readers[i]->AddVariable( "1.00023*ph.idmva_PsEffWidthSigmaRR + 0.0913", &myphoton_ESEffSigmaRR);
221			}
222	mingyang	1.11	}
223
224	mingyang	1.1	}
225
226			fReaderEndcap->BookMVA("BDT method",EndcapWeights);
227			fReaderBarrel->BookMVA("BDT method",BarrelWeights);
228	mingyang	1.4
229	mingyang	1.1	assert(fReaderEndcap);
230			assert(fReaderBarrel);
231
232			}
233
234	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) {
235	mingyang	1.1
236			//initilize the bool value
237			PassMVA=kFALSE;
238
239	mingyang	1.12	Float_t photon_bdt = MVATools::GetMVAbdtValue_2011(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
240	mingyang	1.1
241	mingyang	1.4	if (isbarrel) {
242			if(bdt>bdtCutBarrel){
243			PassMVA=kTRUE;
244			}
245	mingyang	1.1	}
246	mingyang	1.4	else {
247			if(bdt>bdtCutEndcap){
248			PassMVA=kTRUE;
249	mingyang	1.1	}
250	mingyang	1.4	}
251	mingyang	1.1	return PassMVA;
252			}
253
254			//---------------------------------------------------------------------------------
255			Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
256
257			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
258			float cic4_allcuts_temp_sublead[] = {
259			3.8, 2.2, 1.77, 1.29,
260			11.7, 3.4, 3.9, 1.84,
261			3.5, 2.2, 2.3, 1.45,
262			0.0106, 0.0097, 0.028, 0.027,
263			0.082, 0.062, 0.065, 0.048,
264			0.94, 0.36, 0.94, 0.32,
265			1., 0.062, 0.97, 0.97,
266			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
267
268			//initilize the bool value
269			PassElecVetoInt=0;
270
271			dRTrack = PhotonTools::ElectronVetoCiC(p, els);
272
273			ScEta_MVA=p->SCluster()->Eta();
274
275			isbarrel = (fabs(ScEta_MVA)<1.4442);
276
277	bendavid	1.9	R9 = p->R9();
278			//R9 = p->E33()/p->SCluster()->RawEnergy();
279	mingyang	1.1
280			// check which category it is ...
281			_tCat = 1;
282			if ( !isbarrel ) _tCat = 3;
283			if ( R9 < 0.94 ) _tCat++;
284
285			//Electron Veto
286			if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
287			PassElecVetoInt=1;
288			}
289
290			return PassElecVetoInt;
291
292			}
293
294			//--------------------------------------------------------------------------------------------------
295
296	mingyang	1.14	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) {
297	mingyang	1.1
298			//get the variables used to compute MVA variables
299			ecalIso3 = p->EcalRecHitIsoDr03();
300			ecalIso4 = p->EcalRecHitIsoDr04();
301			hcalIso4 = p->HcalTowerSumEtDr04();
302
303			wVtxInd = 0;
304
305	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?
306	bendavid	1.8
307	mingyang	1.1	// track iso worst vtx
308	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) );
309	mingyang	1.1
310			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
311			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
312
313			RawEnergy = p->SCluster()->RawEnergy();
314
315	bendavid	1.9	ScEta = p->SCluster()->Eta();
316
317	mingyang	1.4	//mva varialbes v1 and v2
318	mingyang	1.1	tIso1 = (combIso1) *50./p->Et();
319	bendavid	1.8	tIso3 = (trackIso1)*50./p->Et();
320	mingyang	1.1	tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
321			RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
322			RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
323	mingyang	1.4
324			//compute mva variables for v3
325			HoE = p->HadOverEm();
326			covIEtaIEta = p->CoviEtaiEta();
327			tIso1abs = combIso1;
328	bendavid	1.8	tIso3abs = trackIso1;
329	mingyang	1.4	tIso2abs = combIso2;
330			R9 = p->R9();
331
332			absIsoEcal=(ecalIso3-0.17*_tRho);
333			absIsoHcal=(hcalIso4-0.17*_tRho);
334	mingyang	1.1	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
335			RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
336			RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
337			RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
338			RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
339			RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
340			RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
341			RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
342			RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
343			RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
344			RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
345
346	bendavid	1.10	EtaWidth=p->EtaWidth();
347			PhiWidth=p->PhiWidth();
348	mingyang	1.1	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
349			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
350	mingyang	1.4
351	mingyang	1.1	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
352	mingyang	1.4	NVertexes=vtxCol->GetEntries();
353	mingyang	1.1
354			//spectator variables
355			Pt_MVA=p->Pt();
356			ScEta_MVA=p->SCluster()->Eta();
357	mingyang	1.11
358			//
359
360			isbarrel = (fabs(ScEta_MVA)<1.4442);
361
362			//variable 1201
363			myphoton_pfchargedisogood03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
364			myphoton_pfchargedisobad03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol);
365			myphoton_pfphotoniso03=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
366			myphoton_sieie=covIEtaIEta;
367			myphoton_sieip=CoviEtaiPhi;
368			myphoton_etawidth=EtaWidth;
369			myphoton_phiwidth=PhiWidth;
370			myphoton_r9=R9;
371	mingyang	1.15	myphoton_s4ratio=p->S4Ratio();
372	mingyang	1.11	myphoton_SCeta=ScEta_MVA;
373			event_rho= _tRho;
374
375			myphoton_ESEffSigmaRR=-99;
376
377			if(!isbarrel){
378	mingyang	1.15	myphoton_ESEffSigmaRR=p->EffSigmaRR();
379	mingyang	1.11	}
380	mingyang	1.13
381	mingyang	1.11	if (isbarrel) {
382			reader = fReaderBarrel;
383			}
384			else {
385			reader = fReaderEndcap;
386			}
387
388			assert(reader);
389
390			bdt = reader->EvaluateMVA("BDT method");
391
392			/* printf("HoE: %f\n",HoE);
393			printf("covIEtaIEta: %f\n",covIEtaIEta);
394			printf("tIso1abs: %f\n",tIso1abs);
395			printf("tIso3abs: %f\n",tIso3abs);
396			printf("tIso2abs: %f\n",tIso2abs);
397
398			printf("absIsoEcal: %f\n",absIsoEcal);
399			printf("absIsoHcal: %f\n",absIsoHcal);
400			printf("RelEMax: %f\n",RelEMax);
401			printf("RelETop: %f\n",RelETop);
402			printf("RelEBottom: %f\n",RelEBottom);
403			printf("RelELeft: %f\n",RelELeft);
404			printf("RelERight: %f\n",RelERight);
405			printf("RelE2x5Max: %f\n",RelE2x5Max);
406			printf("RelE2x5Top: %f\n",RelE2x5Top);
407			printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
408			printf("RelE2x5Left: %f\n",RelE2x5Left);
409			printf("RelE2x5Right;: %f\n",RelE2x5Right);
410			printf("RelE5x5: %f\n",RelE5x5);
411	mingyang	1.1
412	mingyang	1.11	printf("EtaWidth: %f\n",EtaWidth);
413			printf("PhiWidth: %f\n",PhiWidth);
414			printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
415			printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
416
417			if (!isbarrel) {
418			printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
419			}*/
420
421			return bdt;
422			}
423
424	mingyang	1.12	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) {
425	mingyang	1.11
426			//get the variables used to compute MVA variables
427			ecalIso3 = p->EcalRecHitIsoDr03();
428			ecalIso4 = p->EcalRecHitIsoDr04();
429			hcalIso4 = p->HcalTowerSumEtDr04();
430
431			wVtxInd = 0;
432
433			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?
434
435			// track iso worst vtx
436			trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
437
438			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
439			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
440
441			RawEnergy = p->SCluster()->RawEnergy();
442
443			ScEta = p->SCluster()->Eta();
444
445			//mva varialbes v1 and v2
446			tIso1 = (combIso1) *50./p->Et();
447			tIso3 = (trackIso1)*50./p->Et();
448			tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
449			RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
450			RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
451
452			//compute mva variables for v3
453			HoE = p->HadOverEm();
454			covIEtaIEta = p->CoviEtaiEta();
455			tIso1abs = combIso1;
456			tIso3abs = trackIso1;
457			tIso2abs = combIso2;
458			R9 = p->R9();
459
460			absIsoEcal=(ecalIso3-0.17*_tRho);
461			absIsoHcal=(hcalIso4-0.17*_tRho);
462			RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
463			RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
464			RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
465			RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
466			RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
467			RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
468			RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
469			RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
470			RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
471			RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
472			RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
473
474			EtaWidth=p->EtaWidth();
475			PhiWidth=p->PhiWidth();
476			CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
477			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
478
479			RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
480			NVertexes=vtxCol->GetEntries();
481
482			//spectator variables
483			Pt_MVA=p->Pt();
484			ScEta_MVA=p->SCluster()->Eta();
485
486			//
487
488	mingyang	1.1	isbarrel = (fabs(ScEta_MVA)<1.4442);
489
490			if (isbarrel) {
491			reader = fReaderBarrel;
492			}
493			else {
494			reader = fReaderEndcap;
495			}
496
497			assert(reader);
498
499			bdt = reader->EvaluateMVA("BDT method");
500
501	mingyang	1.5	/* printf("HoE: %f\n",HoE);
502	mingyang	1.1	printf("covIEtaIEta: %f\n",covIEtaIEta);
503	mingyang	1.4	printf("tIso1abs: %f\n",tIso1abs);
504			printf("tIso3abs: %f\n",tIso3abs);
505			printf("tIso2abs: %f\n",tIso2abs);
506	mingyang	1.1
507	mingyang	1.4	printf("absIsoEcal: %f\n",absIsoEcal);
508			printf("absIsoHcal: %f\n",absIsoHcal);
509	mingyang	1.1	printf("RelEMax: %f\n",RelEMax);
510			printf("RelETop: %f\n",RelETop);
511			printf("RelEBottom: %f\n",RelEBottom);
512			printf("RelELeft: %f\n",RelELeft);
513			printf("RelERight: %f\n",RelERight);
514			printf("RelE2x5Max: %f\n",RelE2x5Max);
515			printf("RelE2x5Top: %f\n",RelE2x5Top);
516			printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
517			printf("RelE2x5Left: %f\n",RelE2x5Left);
518			printf("RelE2x5Right;: %f\n",RelE2x5Right);
519			printf("RelE5x5: %f\n",RelE5x5);
520
521			printf("EtaWidth: %f\n",EtaWidth);
522			printf("PhiWidth: %f\n",PhiWidth);
523			printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
524			printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
525
526	mingyang	1.4	if (!isbarrel) {
527	mingyang	1.1	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
528	mingyang	1.5	}*/
529	mingyang	1.1
530			return bdt;
531			}