Utils/src/MVATools.cc

// $Id: MVATools.cc,v 1.9 2011/12/17 20:00:40 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 );
    }    

    if(VariableType==10){
      readers[i]->AddVariable( "HoE", &HoE );
      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
      readers[i]->AddVariable( "tIso1abs", &tIso1abs );
      readers[i]->AddVariable( "tIso3abs", &tIso3abs );
      readers[i]->AddVariable( "tIso2abs", &tIso2abs );
      readers[i]->AddVariable( "R9", &R9 );
      readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
      readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
      readers[i]->AddVariable( "NVertexes", &NVertexes );
      readers[i]->AddVariable( "ScEta", &ScEta );
      readers[i]->AddVariable( "EtaWidth", &EtaWidth );
      readers[i]->AddVariable( "PhiWidth", &PhiWidth );      
    }    

  }
  
  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->EtaWidth();
  PhiWidth=p->PhiWidth();
  CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi(); 
  CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();

  RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
  NVertexes=vtxCol->GetEntries();
  
  //spectator variables
  Pt_MVA=p->Pt();
  ScEta_MVA=p->SCluster()->Eta();
  
  isbarrel = (fabs(ScEta_MVA)<1.4442);
  
  if (isbarrel) {
    reader = fReaderBarrel;
  }
  else {
    reader = fReaderEndcap;
  }
  
  assert(reader); 

  bdt = reader->EvaluateMVA("BDT method");

  /* printf("HoE: %f\n",HoE);
  printf("covIEtaIEta: %f\n",covIEtaIEta);
  printf("tIso1abs: %f\n",tIso1abs);
  printf("tIso3abs: %f\n",tIso3abs);
  printf("tIso2abs: %f\n",tIso2abs);
  
  printf("absIsoEcal: %f\n",absIsoEcal);
  printf("absIsoHcal: %f\n",absIsoHcal);
  printf("RelEMax: %f\n",RelEMax);
  printf("RelETop: %f\n",RelETop);
  printf("RelEBottom: %f\n",RelEBottom);
  printf("RelELeft: %f\n",RelELeft);
  printf("RelERight: %f\n",RelERight);
  printf("RelE2x5Max: %f\n",RelE2x5Max);
  printf("RelE2x5Top: %f\n",RelE2x5Top);
  printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
  printf("RelE2x5Left: %f\n",RelE2x5Left);
  printf("RelE2x5Right;: %f\n",RelE2x5Right);
  printf("RelE5x5: %f\n",RelE5x5);
  
  printf("EtaWidth: %f\n",EtaWidth);
  printf("PhiWidth: %f\n",PhiWidth);
  printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
  printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
  
  if (!isbarrel) {
    printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
    }*/
  
  return bdt;
}
Revision:	1.10
Committed:	Mon Dec 19 23:45:00 2011 UTC (13 years, 4 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_028, Mit_027, Mit_027a, Mit_025e, Mit_025d
Changes since 1.9:	+18 -3 lines
Log Message:	updated id mva
#	User	Rev	Content
1	bendavid	1.10	// $Id: MVATools.cc,v 1.9 2011/12/17 20:00:40 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	bendavid	1.10	if(VariableType==10){
176			readers[i]->AddVariable( "HoE", &HoE );
177			readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
178			readers[i]->AddVariable( "tIso1abs", &tIso1abs );
179			readers[i]->AddVariable( "tIso3abs", &tIso3abs );
180			readers[i]->AddVariable( "tIso2abs", &tIso2abs );
181			readers[i]->AddVariable( "R9", &R9 );
182			readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
183			readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
184			readers[i]->AddVariable( "NVertexes", &NVertexes );
185			readers[i]->AddVariable( "ScEta", &ScEta );
186			readers[i]->AddVariable( "EtaWidth", &EtaWidth );
187			readers[i]->AddVariable( "PhiWidth", &PhiWidth );
188			}
189
190	mingyang	1.1	}
191
192			fReaderEndcap->BookMVA("BDT method",EndcapWeights);
193			fReaderBarrel->BookMVA("BDT method",BarrelWeights);
194	mingyang	1.4
195	mingyang	1.1	assert(fReaderEndcap);
196			assert(fReaderBarrel);
197
198			}
199
200	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) {
201	mingyang	1.1
202			//initilize the bool value
203			PassMVA=kFALSE;
204
205	bendavid	1.7	Float_t photon_bdt = MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
206	mingyang	1.1
207	mingyang	1.4	if (isbarrel) {
208			if(bdt>bdtCutBarrel){
209			PassMVA=kTRUE;
210			}
211	mingyang	1.1	}
212	mingyang	1.4	else {
213			if(bdt>bdtCutEndcap){
214			PassMVA=kTRUE;
215	mingyang	1.1	}
216	mingyang	1.4	}
217	mingyang	1.1	return PassMVA;
218			}
219
220			//---------------------------------------------------------------------------------
221			Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
222
223			// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
224			float cic4_allcuts_temp_sublead[] = {
225			3.8, 2.2, 1.77, 1.29,
226			11.7, 3.4, 3.9, 1.84,
227			3.5, 2.2, 2.3, 1.45,
228			0.0106, 0.0097, 0.028, 0.027,
229			0.082, 0.062, 0.065, 0.048,
230			0.94, 0.36, 0.94, 0.32,
231			1., 0.062, 0.97, 0.97,
232			1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
233
234			//initilize the bool value
235			PassElecVetoInt=0;
236
237			dRTrack = PhotonTools::ElectronVetoCiC(p, els);
238
239			ScEta_MVA=p->SCluster()->Eta();
240
241			isbarrel = (fabs(ScEta_MVA)<1.4442);
242
243	bendavid	1.9	R9 = p->R9();
244			//R9 = p->E33()/p->SCluster()->RawEnergy();
245	mingyang	1.1
246			// check which category it is ...
247			_tCat = 1;
248			if ( !isbarrel ) _tCat = 3;
249			if ( R9 < 0.94 ) _tCat++;
250
251			//Electron Veto
252			if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
253			PassElecVetoInt=1;
254			}
255
256			return PassElecVetoInt;
257
258			}
259
260			//--------------------------------------------------------------------------------------------------
261
262	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) {
263	mingyang	1.1
264			//get the variables used to compute MVA variables
265			ecalIso3 = p->EcalRecHitIsoDr03();
266			ecalIso4 = p->EcalRecHitIsoDr04();
267			hcalIso4 = p->HcalTowerSumEtDr04();
268
269			wVtxInd = 0;
270
271	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?
272	bendavid	1.8
273	mingyang	1.1	// track iso worst vtx
274	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) );
275	mingyang	1.1
276			combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
277			combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
278
279			RawEnergy = p->SCluster()->RawEnergy();
280
281	bendavid	1.9	ScEta = p->SCluster()->Eta();
282
283	mingyang	1.4	//mva varialbes v1 and v2
284	mingyang	1.1	tIso1 = (combIso1) *50./p->Et();
285	bendavid	1.8	tIso3 = (trackIso1)*50./p->Et();
286	mingyang	1.1	tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
287			RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
288			RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
289	mingyang	1.4
290			//compute mva variables for v3
291			HoE = p->HadOverEm();
292			covIEtaIEta = p->CoviEtaiEta();
293			tIso1abs = combIso1;
294	bendavid	1.8	tIso3abs = trackIso1;
295	mingyang	1.4	tIso2abs = combIso2;
296			R9 = p->R9();
297
298			absIsoEcal=(ecalIso3-0.17*_tRho);
299			absIsoHcal=(hcalIso4-0.17*_tRho);
300	mingyang	1.1	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
301			RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
302			RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
303			RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
304			RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
305			RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
306			RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
307			RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
308			RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
309			RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
310			RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
311
312	bendavid	1.10	EtaWidth=p->EtaWidth();
313			PhiWidth=p->PhiWidth();
314	mingyang	1.1	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
315			CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
316	mingyang	1.4
317	mingyang	1.1	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
318	mingyang	1.4	NVertexes=vtxCol->GetEntries();
319	mingyang	1.1
320			//spectator variables
321			Pt_MVA=p->Pt();
322			ScEta_MVA=p->SCluster()->Eta();
323
324			isbarrel = (fabs(ScEta_MVA)<1.4442);
325
326			if (isbarrel) {
327			reader = fReaderBarrel;
328			}
329			else {
330			reader = fReaderEndcap;
331			}
332
333			assert(reader);
334
335			bdt = reader->EvaluateMVA("BDT method");
336
337	mingyang	1.5	/* printf("HoE: %f\n",HoE);
338	mingyang	1.1	printf("covIEtaIEta: %f\n",covIEtaIEta);
339	mingyang	1.4	printf("tIso1abs: %f\n",tIso1abs);
340			printf("tIso3abs: %f\n",tIso3abs);
341			printf("tIso2abs: %f\n",tIso2abs);
342	mingyang	1.1
343	mingyang	1.4	printf("absIsoEcal: %f\n",absIsoEcal);
344			printf("absIsoHcal: %f\n",absIsoHcal);
345	mingyang	1.1	printf("RelEMax: %f\n",RelEMax);
346			printf("RelETop: %f\n",RelETop);
347			printf("RelEBottom: %f\n",RelEBottom);
348			printf("RelELeft: %f\n",RelELeft);
349			printf("RelERight: %f\n",RelERight);
350			printf("RelE2x5Max: %f\n",RelE2x5Max);
351			printf("RelE2x5Top: %f\n",RelE2x5Top);
352			printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
353			printf("RelE2x5Left: %f\n",RelE2x5Left);
354			printf("RelE2x5Right;: %f\n",RelE2x5Right);
355			printf("RelE5x5: %f\n",RelE5x5);
356
357			printf("EtaWidth: %f\n",EtaWidth);
358			printf("PhiWidth: %f\n",PhiWidth);
359			printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
360			printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
361
362	mingyang	1.4	if (!isbarrel) {
363	mingyang	1.1	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
364	mingyang	1.5	}*/
365	mingyang	1.1
366			return bdt;
367			}