Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.1 2011/10/13 17:28:29 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) {
  
  // 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);//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);
  
  // 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]){
    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.2
Committed:	Thu Oct 13 18:34:50 2011 UTC (13 years, 6 months ago) by mingyang
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+4 -4 lines
Log Message:	bdtcut as a parameter
#	Content
1	// $Id: MVATools.cc,v 1.1 2011/10/13 17:28:29 mingyang Exp $
2
3	#include "MitPhysics/Utils/interface/PhotonTools.h"
4	#include "MitPhysics/Utils/interface/MVATools.h"
5	#include "MitPhysics/Utils/interface/ElectronTools.h"
6	#include "MitPhysics/Utils/interface/IsolationTools.h"
7	#include "MitAna/DataTree/interface/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	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) {
120
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	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?
147
148	// track iso only
149	trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
150
151	// track iso worst vtx
152	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol);
153
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	//Electron Veto
204	if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
205	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	if(bdt>bdtCutBarrel){
221	PassMVA=kTRUE;
222	}
223	}
224	else {
225	if(bdt>bdtCutEndcap){
226	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	}