Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.8 2011/12/13 21:13:23 bendavid Exp $

#include "MitPhysics/Utils/interface/PhotonTools.h"
#include "MitPhysics/Utils/interface/MVATools.h"
#include "MitPhysics/Utils/interface/ElectronTools.h"
#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitAna/DataTree/interface/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),
  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 );
    }    

  }
  
  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(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
  
  if (isbarrel) {
    if(bdt>bdtCutBarrel){
      PassMVA=kTRUE;    
    }
  }
  else {
    if(bdt>bdtCutEndcap){
      PassMVA=kTRUE;    
    }
  } 
  return PassMVA;
}

//---------------------------------------------------------------------------------
Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
  
  // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
  float cic4_allcuts_temp_sublead[] = { 
    3.8,         2.2,         1.77,        1.29,
    11.7,        3.4,         3.9,         1.84,
    3.5,         2.2,         2.3,         1.45,
    0.0106,      0.0097,      0.028,       0.027,
    0.082,       0.062,       0.065,       0.048,
    0.94,        0.36,        0.94,        0.32,
    1.,          0.062,       0.97,        0.97,
    1.5,         1.5,         1.5,         1.5 };  // the last line is PixelmatchVeto and un-used
  
  //initilize the bool value
  PassElecVetoInt=0;
  
  dRTrack = PhotonTools::ElectronVetoCiC(p, els);
  
  ScEta_MVA=p->SCluster()->Eta();
  
  isbarrel = (fabs(ScEta_MVA)<1.4442);

  R9 = p->R9();
  //R9 = p->E33()/p->SCluster()->RawEnergy();
  
  // check which category it is ...
  _tCat = 1;
  if ( !isbarrel ) _tCat = 3;
  if ( R9 < 0.94 ) _tCat++;
  
  //Electron Veto
  if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
    PassElecVetoInt=1;
  }
  
  return  PassElecVetoInt;
  
}

//--------------------------------------------------------------------------------------------------

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