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
#	Content
1	// $Id: MVATools.cc,v 1.9 2011/12/17 20:00:40 bendavid 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 v4
24	HoE(0),
25	covIEtaIEta(0),
26	tIso1abs(0),
27	tIso3abs(0),
28	tIso2abs(0),
29	R9(0),
30
31	absIsoEcal(0),
32	absIsoHcal(0),
33	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	CoviEtaiPhi(0),
48	CoviPhiiPhi(0),
49
50	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	ScEta(0.)
60	{
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
79	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	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	if(VariableType==1\|\|VariableType==2){
98	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	readers[i]->AddVariable( "RelELeft", &RelELeft );
120	readers[i]->AddVariable( "RelERight", &RelERight );
121	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	if(VariableType==2\|\|VariableType==3\|\|VariableType==4){
130	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
131	readers[i]->AddVariable( "PhiWidth", &PhiWidth );
132	readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
133	readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
134	if(VariableType==4){
135	readers[i]->AddVariable( "NVertexes", &NVertexes );
136	}
137	if(i==0){
138	readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
139	}
140	}
141
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
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
175	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	}
191
192	fReaderEndcap->BookMVA("BDT method",EndcapWeights);
193	fReaderBarrel->BookMVA("BDT method",BarrelWeights);
194
195	assert(fReaderEndcap);
196	assert(fReaderBarrel);
197
198	}
199
200	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
202	//initilize the bool value
203	PassMVA=kFALSE;
204
205	Float_t photon_bdt = MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
206
207	if (isbarrel) {
208	if(bdt>bdtCutBarrel){
209	PassMVA=kTRUE;
210	}
211	}
212	else {
213	if(bdt>bdtCutEndcap){
214	PassMVA=kTRUE;
215	}
216	}
217	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	R9 = p->R9();
244	//R9 = p->E33()/p->SCluster()->RawEnergy();
245
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	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
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	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
273	// track iso worst vtx
274	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
275
276	combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
277	combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
278
279	RawEnergy = p->SCluster()->RawEnergy();
280
281	ScEta = p->SCluster()->Eta();
282
283	//mva varialbes v1 and v2
284	tIso1 = (combIso1) *50./p->Et();
285	tIso3 = (trackIso1)*50./p->Et();
286	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
290	//compute mva variables for v3
291	HoE = p->HadOverEm();
292	covIEtaIEta = p->CoviEtaiEta();
293	tIso1abs = combIso1;
294	tIso3abs = trackIso1;
295	tIso2abs = combIso2;
296	R9 = p->R9();
297
298	absIsoEcal=(ecalIso3-0.17*_tRho);
299	absIsoHcal=(hcalIso4-0.17*_tRho);
300	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	EtaWidth=p->EtaWidth();
313	PhiWidth=p->PhiWidth();
314	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
315	CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
316
317	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
318	NVertexes=vtxCol->GetEntries();
319
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	/* printf("HoE: %f\n",HoE);
338	printf("covIEtaIEta: %f\n",covIEtaIEta);
339	printf("tIso1abs: %f\n",tIso1abs);
340	printf("tIso3abs: %f\n",tIso3abs);
341	printf("tIso2abs: %f\n",tIso2abs);
342
343	printf("absIsoEcal: %f\n",absIsoEcal);
344	printf("absIsoHcal: %f\n",absIsoHcal);
345	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	if (!isbarrel) {
363	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
364	}*/
365
366	return bdt;
367	}