Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.2 2011/10/13 18:34:50 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/StableData.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 0
  HoE(0),
  covIEtaIEta(0),
  tIso1(0),
  tIso3(0),
  tIso2(0),
  R9(0),
  
  //MVA variables 1
  RelIsoEcal(0),
  RelIsoHcal(0),
  
  RelEMax(0),
  RelETop(0),
  RelEBottom(0),
  RelELeft(0),
  RelERight(0),
  RelE2x5Max(0),
  RelE2x5Top(0),
  RelE2x5Bottom(0),
  RelE2x5Left(0),
  RelE2x5Right(0),
  RelE5x5(0),
  
  //MVA variables 2
  EtaWidth(0),
  PhiWidth(0),
  CoviEtaiPhi(0), 
  CoviPhiiPhi(0),
  RelPreshowerEnergy(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==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==2){
      readers[i]->AddVariable( "EtaWidth", &EtaWidth );
      readers[i]->AddVariable( "PhiWidth", &PhiWidth );
      readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
      readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
      if(i==0){
        readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
      }
    }
  }
  
  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,const ElectronCol* els,double MVAPtMin, Float_t bdtCutBarrel, Float_t bdtCutEndcap, Bool_t applyEleVeto) {
  
  // 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
  PassMVA=kFALSE;
  PassElecVeto=kFALSE;
  //Bool_t PassElecVeto=PhotonTools::PassElectronVeto(p,els);
  
  
  //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, (!applyEleVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
  
  // track iso only
  trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyEleVeto ? els : NULL));
  
  // track iso worst vtx
  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyEleVeto ? els : NULL));
  
  combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
  combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
  
  RawEnergy = p->SCluster()->RawEnergy();
  
  dRTrack = PhotonTools::ElectronVetoCiC(p, els);
  
  //compute MVA variables 0
  HoE = p->HadOverEm();
  covIEtaIEta = p->CoviEtaiEta();
  tIso1 = (combIso1) *50./p->Et();
  tIso3 = (trackIso3)*50./p->Et();
  tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
  R9 = p->R9();
  
  //newly added MVA variables 1
  RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
  RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
  
  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;
  
  //newly added MVA variables 2
  EtaWidth=p->SCluster()->EtaWidth();
  PhiWidth=p->SCluster()->PhiWidth();
  CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi(); 
  CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
  RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
  
  //spectator variables
  Pt_MVA=p->Pt();
  ScEta_MVA=p->SCluster()->Eta();
  
  isbarrel = (fabs(ScEta_MVA)<1.4442);
  
  // check which category it is ...
  _tCat = 1;
  if ( !isbarrel ) _tCat = 3;
  if ( R9 < 0.94 ) _tCat++;
  
  //Electron Veto (made electron Veto optinal (Fabian) )
  if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] || !applyEleVeto ){
    PassElecVeto=kTRUE;
  }
  
  if(PassElecVeto && Pt_MVA>MVAPtMin && ((fabs(ScEta_MVA)<1.4442)||(fabs(ScEta_MVA)>1.566 && fabs(ScEta_MVA)<2.5)) && tIso1<250 && tIso2<250 && tIso3<250){
    
    if (isbarrel) {
      reader = fReaderBarrel;
    }
    else {
      reader = fReaderEndcap;
    }
    
    bdt = reader->EvaluateMVA( "BDT method" );
    
    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();
  
  // 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 ElectronCol* els) {
  
  //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);//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
  
  // track iso only
  trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
  
  // track iso worst vtx
  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol);
  
  combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
  combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
  
  RawEnergy = p->SCluster()->RawEnergy();
  
  dRTrack = PhotonTools::ElectronVetoCiC(p, els);
  
  //compute MVA variables 0
  HoE = p->HadOverEm();
  covIEtaIEta = p->CoviEtaiEta();
  tIso1 = (combIso1) *50./p->Et();
  tIso3 = (trackIso3)*50./p->Et();
  tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
  R9 = p->R9();
  
  //newly added MVA variables 1
  RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
  RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
  
  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;
  
  //newly added MVA variables 2
  EtaWidth=p->SCluster()->EtaWidth();
  PhiWidth=p->SCluster()->PhiWidth();
  CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi(); 
  CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
  RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
  
  //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("tIso1: %f\n",tIso1);
  printf("tIso3: %f\n",tIso3);
  printf("tIso2: %f\n",tIso2);
  printf("tIso1: %f\n",tIso1);
  printf("R9: %f\n",R9);
  printf("RelIsoEcal: %f\n",RelIsoEcal);
  printf("RelIsoHcal: %f\n",RelIsoHcal);
  
  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.3
Committed:	Tue Oct 18 11:27:19 2011 UTC (13 years, 6 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c, Mit_025b, Mit_025a
Changes since 1.2:	+7 -7 lines
Log Message:	Added electron Veto inversion to CiCTrackIso & MVA Tools
#	User	Rev	Content
1	fabstoec	1.3	// $Id: MVATools.cc,v 1.2 2011/10/13 18:34:50 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			#include "MitAna/DataTree/interface/StableData.h"
8			#include <TFile.h>
9			#include <TRandom3.h>
10			#include "TMVA/Tools.h"//MVA
11			#include "TMVA/Reader.h"//MVA
12
13			ClassImp(mithep::MVATools)
14
15			using namespace mithep;
16
17
18			//--------------------------------------------------------------------------------------------------
19			MVATools::MVATools():
20			fReaderEndcap(0),
21			fReaderBarrel(0),
22
23			//MVA variables 0
24			HoE(0),
25			covIEtaIEta(0),
26			tIso1(0),
27			tIso3(0),
28			tIso2(0),
29			R9(0),
30
31			//MVA variables 1
32			RelIsoEcal(0),
33			RelIsoHcal(0),
34
35			RelEMax(0),
36			RelETop(0),
37			RelEBottom(0),
38			RelELeft(0),
39			RelERight(0),
40			RelE2x5Max(0),
41			RelE2x5Top(0),
42			RelE2x5Bottom(0),
43			RelE2x5Left(0),
44			RelE2x5Right(0),
45			RelE5x5(0),
46
47			//MVA variables 2
48			EtaWidth(0),
49			PhiWidth(0),
50			CoviEtaiPhi(0),
51			CoviPhiiPhi(0),
52			RelPreshowerEnergy(0)
53
54			{
55			// Constructor.
56			}
57
58			//--------------------------------------------------------------------------------------------------
59			void MVATools::InitializeMVA(int VariableType, TString EndcapWeights,TString BarrelWeights) {
60
61			if (fReaderEndcap) delete fReaderEndcap;
62			if (fReaderBarrel) delete fReaderBarrel;
63
64			fReaderEndcap = new TMVA::Reader( "!Color:!Silent:Error" );
65			fReaderBarrel = new TMVA::Reader( "!Color:!Silent:Error" );
66
67			TMVA::Reader *readers[2];
68			readers[0] = fReaderEndcap;
69			readers[1] = fReaderBarrel;
70
71			for (UInt_t i=0; i<2; ++i) {
72
73			if(VariableType==0\|\|VariableType==1\|\|VariableType==2){
74			readers[i]->AddVariable( "HoE", &HoE );
75			readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
76			readers[i]->AddVariable( "tIso1", &tIso1 );
77			readers[i]->AddVariable( "tIso3", &tIso3 );
78			readers[i]->AddVariable( "tIso2", &tIso2 );
79			readers[i]->AddVariable( "R9", &R9 );
80			}
81
82			if(VariableType==1\|\|VariableType==2){
83			readers[i]->AddVariable( "RelIsoEcal", &RelIsoEcal);
84			readers[i]->AddVariable( "RelIsoHcal", &RelIsoHcal);
85			readers[i]->AddVariable( "RelEMax", &RelEMax);
86			readers[i]->AddVariable( "RelETop",&RelETop);
87			readers[i]->AddVariable( "RelEBottom", &RelEBottom);
88			readers[i]->AddVariable( "RelELeft", &RelELeft );
89			readers[i]->AddVariable( "RelERight", &RelERight);
90			readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
91			readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
92			readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
93			readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
94			readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
95			readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
96			}
97
98			if(VariableType==2){
99			readers[i]->AddVariable( "EtaWidth", &EtaWidth );
100			readers[i]->AddVariable( "PhiWidth", &PhiWidth );
101			readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
102			readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
103			if(i==0){
104			readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
105			}
106			}
107			}
108
109			fReaderEndcap->BookMVA("BDT method",EndcapWeights);
110			fReaderBarrel->BookMVA("BDT method",BarrelWeights);
111
112			assert(fReaderEndcap);
113			assert(fReaderBarrel);
114
115			}
116
117			//--------------------------------------------------------------------------------------------------
118
119	fabstoec	1.3	Bool_t MVATools::PassMVASelection(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els,double MVAPtMin, Float_t bdtCutBarrel, Float_t bdtCutEndcap, Bool_t applyEleVeto) {
120	mingyang	1.1
121			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
122			float cic4_allcuts_temp_sublead[] = {
123			3.8, 2.2, 1.77, 1.29,
124			11.7, 3.4, 3.9, 1.84,
125			3.5, 2.2, 2.3, 1.45,
126			0.0106, 0.0097, 0.028, 0.027,
127			0.082, 0.062, 0.065, 0.048,
128			0.94, 0.36, 0.94, 0.32,
129			1., 0.062, 0.97, 0.97,
130			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
131
132
133			//initilize the bool value
134			PassMVA=kFALSE;
135			PassElecVeto=kFALSE;
136			//Bool_t PassElecVeto=PhotonTools::PassElectronVeto(p,els);
137
138
139			//get the variables used to compute MVA variables
140			ecalIso3 = p->EcalRecHitIsoDr03();
141			ecalIso4 = p->EcalRecHitIsoDr04();
142			hcalIso4 = p->HcalTowerSumEtDr04();
143
144			wVtxInd = 0;
145
146	fabstoec	1.3	trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyEleVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
147	mingyang	1.1
148			// track iso only
149	fabstoec	1.3	trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyEleVeto ? els : NULL));
150	mingyang	1.1
151			// track iso worst vtx
152	fabstoec	1.3	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyEleVeto ? els : NULL));
153	mingyang	1.1
154			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
155			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
156
157			RawEnergy = p->SCluster()->RawEnergy();
158
159			dRTrack = PhotonTools::ElectronVetoCiC(p, els);
160
161			//compute MVA variables 0
162			HoE = p->HadOverEm();
163			covIEtaIEta = p->CoviEtaiEta();
164			tIso1 = (combIso1) *50./p->Et();
165			tIso3 = (trackIso3)*50./p->Et();
166			tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
167			R9 = p->R9();
168
169			//newly added MVA variables 1
170			RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
171			RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
172
173			RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
174			RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
175			RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
176			RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
177			RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
178			RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
179			RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
180			RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
181			RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
182			RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
183			RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
184
185			//newly added MVA variables 2
186			EtaWidth=p->SCluster()->EtaWidth();
187			PhiWidth=p->SCluster()->PhiWidth();
188			CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
189			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
190			RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
191
192			//spectator variables
193			Pt_MVA=p->Pt();
194			ScEta_MVA=p->SCluster()->Eta();
195
196			isbarrel = (fabs(ScEta_MVA)<1.4442);
197
198			// check which category it is ...
199			_tCat = 1;
200			if ( !isbarrel ) _tCat = 3;
201			if ( R9 < 0.94 ) _tCat++;
202
203	fabstoec	1.3	//Electron Veto (made electron Veto optinal (Fabian) )
204			if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] \|\| !applyEleVeto ){
205	mingyang	1.1	PassElecVeto=kTRUE;
206			}
207
208			if(PassElecVeto && Pt_MVA>MVAPtMin && ((fabs(ScEta_MVA)<1.4442)\|\|(fabs(ScEta_MVA)>1.566 && fabs(ScEta_MVA)<2.5)) && tIso1<250 && tIso2<250 && tIso3<250){
209
210			if (isbarrel) {
211			reader = fReaderBarrel;
212			}
213			else {
214			reader = fReaderEndcap;
215			}
216
217			bdt = reader->EvaluateMVA( "BDT method" );
218
219			if (isbarrel) {
220	mingyang	1.2	if(bdt>bdtCutBarrel){
221	mingyang	1.1	PassMVA=kTRUE;
222			}
223			}
224			else {
225	mingyang	1.2	if(bdt>bdtCutEndcap){
226	mingyang	1.1	PassMVA=kTRUE;
227			}
228			}
229
230			}
231
232			return PassMVA;
233			}
234
235			//---------------------------------------------------------------------------------
236			Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
237
238			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
239			float cic4_allcuts_temp_sublead[] = {
240			3.8, 2.2, 1.77, 1.29,
241			11.7, 3.4, 3.9, 1.84,
242			3.5, 2.2, 2.3, 1.45,
243			0.0106, 0.0097, 0.028, 0.027,
244			0.082, 0.062, 0.065, 0.048,
245			0.94, 0.36, 0.94, 0.32,
246			1., 0.062, 0.97, 0.97,
247			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
248
249			//initilize the bool value
250			PassElecVetoInt=0;
251
252			dRTrack = PhotonTools::ElectronVetoCiC(p, els);
253
254			ScEta_MVA=p->SCluster()->Eta();
255
256			isbarrel = (fabs(ScEta_MVA)<1.4442);
257
258			R9 = p->R9();
259
260			// check which category it is ...
261			_tCat = 1;
262			if ( !isbarrel ) _tCat = 3;
263			if ( R9 < 0.94 ) _tCat++;
264
265			//Electron Veto
266			if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
267			PassElecVetoInt=1;
268			}
269
270			return PassElecVetoInt;
271
272			}
273
274			//--------------------------------------------------------------------------------------------------
275
276			Float_t MVATools::GetMVAbdtValue(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els) {
277
278			//get the variables used to compute MVA variables
279			ecalIso3 = p->EcalRecHitIsoDr03();
280			ecalIso4 = p->EcalRecHitIsoDr04();
281			hcalIso4 = p->HcalTowerSumEtDr04();
282
283			wVtxInd = 0;
284
285			trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
286
287			// track iso only
288			trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
289
290			// track iso worst vtx
291			trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol);
292
293			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
294			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
295
296			RawEnergy = p->SCluster()->RawEnergy();
297
298			dRTrack = PhotonTools::ElectronVetoCiC(p, els);
299
300			//compute MVA variables 0
301			HoE = p->HadOverEm();
302			covIEtaIEta = p->CoviEtaiEta();
303			tIso1 = (combIso1) *50./p->Et();
304			tIso3 = (trackIso3)*50./p->Et();
305			tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
306			R9 = p->R9();
307
308			//newly added MVA variables 1
309			RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
310			RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
311
312			RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
313			RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
314			RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
315			RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
316			RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
317			RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
318			RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
319			RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
320			RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
321			RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
322			RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
323
324			//newly added MVA variables 2
325			EtaWidth=p->SCluster()->EtaWidth();
326			PhiWidth=p->SCluster()->PhiWidth();
327			CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
328			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
329			RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
330
331			//spectator variables
332			Pt_MVA=p->Pt();
333			ScEta_MVA=p->SCluster()->Eta();
334
335			isbarrel = (fabs(ScEta_MVA)<1.4442);
336
337			if (isbarrel) {
338			reader = fReaderBarrel;
339			}
340			else {
341			reader = fReaderEndcap;
342			}
343
344			assert(reader);
345
346			bdt = reader->EvaluateMVA("BDT method");
347
348			/* printf("HoE: %f\n",HoE);
349			printf("covIEtaIEta: %f\n",covIEtaIEta);
350			printf("tIso1: %f\n",tIso1);
351			printf("tIso3: %f\n",tIso3);
352			printf("tIso2: %f\n",tIso2);
353			printf("tIso1: %f\n",tIso1);
354			printf("R9: %f\n",R9);
355			printf("RelIsoEcal: %f\n",RelIsoEcal);
356			printf("RelIsoHcal: %f\n",RelIsoHcal);
357
358			printf("RelEMax: %f\n",RelEMax);
359			printf("RelETop: %f\n",RelETop);
360			printf("RelEBottom: %f\n",RelEBottom);
361			printf("RelELeft: %f\n",RelELeft);
362			printf("RelERight: %f\n",RelERight);
363			printf("RelE2x5Max: %f\n",RelE2x5Max);
364			printf("RelE2x5Top: %f\n",RelE2x5Top);
365			printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
366			printf("RelE2x5Left: %f\n",RelE2x5Left);
367			printf("RelE2x5Right;: %f\n",RelE2x5Right);
368			printf("RelE5x5: %f\n",RelE5x5);
369
370			printf("EtaWidth: %f\n",EtaWidth);
371			printf("PhiWidth: %f\n",PhiWidth);
372			printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
373			printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
374
375			if (isbarrel) {
376			printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
377			}*/
378
379			return bdt;
380			}