Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.4 2011/12/05 00:47:46 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 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)
{
  // 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 );
      }
    }
  }
  
  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) {
  
  //initilize the bool value
  PassMVA=kFALSE;
  
  Float_t photon_bdt =  MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho);
  
  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) {
  
  //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();
  
  //mva varialbes v1 and v2
  tIso1 = (combIso1) *50./p->Et();
  tIso3 = (trackIso3)*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 = trackIso3;
  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->SCluster()->EtaWidth();
  PhiWidth=p->SCluster()->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.5
Committed:	Mon Dec 5 01:39:48 2011 UTC (13 years, 5 months ago) by mingyang
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+6 -8 lines
Log Message:	change the getbdt and passmva function
#	User	Rev	Content
1	mingyang	1.5	// $Id: MVATools.cc,v 1.4 2011/12/05 00:47:46 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	mingyang	1.4	//MVA Variables v4
24	mingyang	1.1	HoE(0),
25			covIEtaIEta(0),
26	mingyang	1.4	tIso1abs(0),
27			tIso3abs(0),
28			tIso2abs(0),
29	mingyang	1.1	R9(0),
30
31	mingyang	1.4	absIsoEcal(0),
32			absIsoHcal(0),
33	mingyang	1.1	RelEMax(0),
34			RelETop(0),
35			RelEBottom(0),
36			RelELeft(0),
37			RelERight(0),
38			RelE2x5Max(0),
39			RelE2x5Top(0),
40			RelE2x5Bottom(0),
41			RelE2x5Left(0),
42			RelE2x5Right(0),
43			RelE5x5(0),
44
45			EtaWidth(0),
46			PhiWidth(0),
47	mingyang	1.4	CoviEtaiPhi(0),
48	mingyang	1.1	CoviPhiiPhi(0),
49
50	mingyang	1.4	NVertexes(0),
51			RelPreshowerEnergy(0),
52
53			//MVA v2 and v1
54			RelIsoEcal(0),
55			RelIsoHcal(0),
56			tIso1(0),
57			tIso3(0),
58			tIso2(0)
59	mingyang	1.1	{
60			// Constructor.
61			}
62
63			//--------------------------------------------------------------------------------------------------
64			void MVATools::InitializeMVA(int VariableType, TString EndcapWeights,TString BarrelWeights) {
65
66			if (fReaderEndcap) delete fReaderEndcap;
67			if (fReaderBarrel) delete fReaderBarrel;
68
69			fReaderEndcap = new TMVA::Reader( "!Color:!Silent:Error" );
70			fReaderBarrel = new TMVA::Reader( "!Color:!Silent:Error" );
71
72			TMVA::Reader *readers[2];
73			readers[0] = fReaderEndcap;
74			readers[1] = fReaderBarrel;
75
76			for (UInt_t i=0; i<2; ++i) {
77	mingyang	1.4
78	mingyang	1.1	if(VariableType==0\|\|VariableType==1\|\|VariableType==2){
79			readers[i]->AddVariable( "HoE", &HoE );
80			readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
81			readers[i]->AddVariable( "tIso1", &tIso1 );
82			readers[i]->AddVariable( "tIso3", &tIso3 );
83			readers[i]->AddVariable( "tIso2", &tIso2 );
84			readers[i]->AddVariable( "R9", &R9 );
85			}
86
87	mingyang	1.4	if(VariableType==3\|\|VariableType==4){
88			readers[i]->AddVariable( "HoE", &HoE );
89			readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
90			readers[i]->AddVariable( "tIso1abs", &tIso1abs );
91			readers[i]->AddVariable( "tIso3abs", &tIso3abs );
92			readers[i]->AddVariable( "tIso2abs", &tIso2abs );
93			readers[i]->AddVariable( "R9", &R9 );
94			}
95
96	mingyang	1.1	if(VariableType==1\|\|VariableType==2){
97	mingyang	1.4	readers[i]->AddVariable( "RelIsoEcal", &RelIsoEcal );
98			readers[i]->AddVariable( "RelIsoHcal", &RelIsoHcal );
99			readers[i]->AddVariable( "RelEMax", &RelEMax );
100			readers[i]->AddVariable( "RelETop", &RelETop );
101			readers[i]->AddVariable( "RelEBottom", &RelEBottom );
102			readers[i]->AddVariable( "RelELeft", &RelELeft );
103			readers[i]->AddVariable( "RelERight", &RelERight );
104			readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
105			readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
106			readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
107			readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
108			readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
109			readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
110			}
111
112			if(VariableType==3\|\|VariableType==4){
113			readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
114			readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
115			readers[i]->AddVariable( "RelEMax", &RelEMax );
116			readers[i]->AddVariable( "RelETop", &RelEBottom);
117			readers[i]->AddVariable( "RelEBottom", &RelEBottom );
118	mingyang	1.1	readers[i]->AddVariable( "RelELeft", &RelELeft );
119	mingyang	1.4	readers[i]->AddVariable( "RelERight", &RelERight );
120	mingyang	1.1	readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
121			readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
122			readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
123			readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
124			readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
125			readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
126			}
127
128	mingyang	1.4	if(VariableType==2\|\|VariableType==3\|\|VariableType==4){
129	mingyang	1.1	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
130			readers[i]->AddVariable( "PhiWidth", &PhiWidth );
131			readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
132			readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
133	mingyang	1.4	if(VariableType==4){
134			readers[i]->AddVariable( "NVertexes", &NVertexes );
135			}
136	mingyang	1.1	if(i==0){
137			readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
138			}
139			}
140			}
141
142			fReaderEndcap->BookMVA("BDT method",EndcapWeights);
143			fReaderBarrel->BookMVA("BDT method",BarrelWeights);
144	mingyang	1.4
145	mingyang	1.1	assert(fReaderEndcap);
146			assert(fReaderBarrel);
147
148			}
149
150	mingyang	1.5	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) {
151	mingyang	1.1
152			//initilize the bool value
153			PassMVA=kFALSE;
154
155	mingyang	1.5	Float_t photon_bdt = MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho);
156	mingyang	1.1
157	mingyang	1.4	if (isbarrel) {
158			if(bdt>bdtCutBarrel){
159			PassMVA=kTRUE;
160			}
161	mingyang	1.1	}
162	mingyang	1.4	else {
163			if(bdt>bdtCutEndcap){
164			PassMVA=kTRUE;
165	mingyang	1.1	}
166	mingyang	1.4	}
167	mingyang	1.1	return PassMVA;
168			}
169
170			//---------------------------------------------------------------------------------
171			Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
172
173			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
174			float cic4_allcuts_temp_sublead[] = {
175			3.8, 2.2, 1.77, 1.29,
176			11.7, 3.4, 3.9, 1.84,
177			3.5, 2.2, 2.3, 1.45,
178			0.0106, 0.0097, 0.028, 0.027,
179			0.082, 0.062, 0.065, 0.048,
180			0.94, 0.36, 0.94, 0.32,
181			1., 0.062, 0.97, 0.97,
182			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
183
184			//initilize the bool value
185			PassElecVetoInt=0;
186
187			dRTrack = PhotonTools::ElectronVetoCiC(p, els);
188
189			ScEta_MVA=p->SCluster()->Eta();
190
191			isbarrel = (fabs(ScEta_MVA)<1.4442);
192
193			R9 = p->R9();
194
195			// check which category it is ...
196			_tCat = 1;
197			if ( !isbarrel ) _tCat = 3;
198			if ( R9 < 0.94 ) _tCat++;
199
200			//Electron Veto
201			if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
202			PassElecVetoInt=1;
203			}
204
205			return PassElecVetoInt;
206
207			}
208
209			//--------------------------------------------------------------------------------------------------
210
211	mingyang	1.5	Float_t MVATools::GetMVAbdtValue(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho) {
212	mingyang	1.1
213			//get the variables used to compute MVA variables
214			ecalIso3 = p->EcalRecHitIsoDr03();
215			ecalIso4 = p->EcalRecHitIsoDr04();
216			hcalIso4 = p->HcalTowerSumEtDr04();
217
218			wVtxInd = 0;
219
220			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?
221
222			// track iso only
223			trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
224
225			// track iso worst vtx
226			trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol);
227
228			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
229			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
230
231			RawEnergy = p->SCluster()->RawEnergy();
232
233	mingyang	1.4	//mva varialbes v1 and v2
234	mingyang	1.1	tIso1 = (combIso1) *50./p->Et();
235			tIso3 = (trackIso3)*50./p->Et();
236			tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
237			RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
238			RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
239	mingyang	1.4
240			//compute mva variables for v3
241			HoE = p->HadOverEm();
242			covIEtaIEta = p->CoviEtaiEta();
243			tIso1abs = combIso1;
244			tIso3abs = trackIso3;
245			tIso2abs = combIso2;
246			R9 = p->R9();
247
248			absIsoEcal=(ecalIso3-0.17*_tRho);
249			absIsoHcal=(hcalIso4-0.17*_tRho);
250	mingyang	1.1	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
251			RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
252			RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
253			RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
254			RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
255			RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
256			RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
257			RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
258			RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
259			RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
260			RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
261
262			EtaWidth=p->SCluster()->EtaWidth();
263			PhiWidth=p->SCluster()->PhiWidth();
264			CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
265			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
266	mingyang	1.4
267	mingyang	1.1	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
268	mingyang	1.4	NVertexes=vtxCol->GetEntries();
269	mingyang	1.1
270			//spectator variables
271			Pt_MVA=p->Pt();
272			ScEta_MVA=p->SCluster()->Eta();
273
274			isbarrel = (fabs(ScEta_MVA)<1.4442);
275
276			if (isbarrel) {
277			reader = fReaderBarrel;
278			}
279			else {
280			reader = fReaderEndcap;
281			}
282
283			assert(reader);
284
285			bdt = reader->EvaluateMVA("BDT method");
286
287	mingyang	1.5	/* printf("HoE: %f\n",HoE);
288	mingyang	1.1	printf("covIEtaIEta: %f\n",covIEtaIEta);
289	mingyang	1.4	printf("tIso1abs: %f\n",tIso1abs);
290			printf("tIso3abs: %f\n",tIso3abs);
291			printf("tIso2abs: %f\n",tIso2abs);
292	mingyang	1.1
293	mingyang	1.4	printf("absIsoEcal: %f\n",absIsoEcal);
294			printf("absIsoHcal: %f\n",absIsoHcal);
295	mingyang	1.1	printf("RelEMax: %f\n",RelEMax);
296			printf("RelETop: %f\n",RelETop);
297			printf("RelEBottom: %f\n",RelEBottom);
298			printf("RelELeft: %f\n",RelELeft);
299			printf("RelERight: %f\n",RelERight);
300			printf("RelE2x5Max: %f\n",RelE2x5Max);
301			printf("RelE2x5Top: %f\n",RelE2x5Top);
302			printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
303			printf("RelE2x5Left: %f\n",RelE2x5Left);
304			printf("RelE2x5Right;: %f\n",RelE2x5Right);
305			printf("RelE5x5: %f\n",RelE5x5);
306
307			printf("EtaWidth: %f\n",EtaWidth);
308			printf("PhiWidth: %f\n",PhiWidth);
309			printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
310			printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
311
312	mingyang	1.4	if (!isbarrel) {
313	mingyang	1.1	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
314	mingyang	1.5	}*/
315	mingyang	1.1
316			return bdt;
317			}