| 1 |
// $Id: MVATools.cc,v 1.4 2011/12/05 00:47:46 mingyang Exp $
|
| 2 |
|
| 3 |
#include "MitPhysics/Utils/interface/PhotonTools.h"
|
| 4 |
#include "MitPhysics/Utils/interface/MVATools.h"
|
| 5 |
#include "MitPhysics/Utils/interface/ElectronTools.h"
|
| 6 |
#include "MitPhysics/Utils/interface/IsolationTools.h"
|
| 7 |
#include "MitAna/DataTree/interface/StableData.h"
|
| 8 |
#include <TFile.h>
|
| 9 |
#include <TRandom3.h>
|
| 10 |
#include "TMVA/Tools.h"//MVA
|
| 11 |
#include "TMVA/Reader.h"//MVA
|
| 12 |
|
| 13 |
ClassImp(mithep::MVATools)
|
| 14 |
|
| 15 |
using namespace mithep;
|
| 16 |
|
| 17 |
|
| 18 |
//--------------------------------------------------------------------------------------------------
|
| 19 |
MVATools::MVATools():
|
| 20 |
fReaderEndcap(0),
|
| 21 |
fReaderBarrel(0),
|
| 22 |
|
| 23 |
//MVA Variables 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 |
{
|
| 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 |
|
| 78 |
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 |
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 |
if(VariableType==1||VariableType==2){
|
| 97 |
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 |
readers[i]->AddVariable( "RelELeft", &RelELeft );
|
| 119 |
readers[i]->AddVariable( "RelERight", &RelERight );
|
| 120 |
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 |
if(VariableType==2||VariableType==3||VariableType==4){
|
| 129 |
readers[i]->AddVariable( "EtaWidth", &EtaWidth );
|
| 130 |
readers[i]->AddVariable( "PhiWidth", &PhiWidth );
|
| 131 |
readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
|
| 132 |
readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
|
| 133 |
if(VariableType==4){
|
| 134 |
readers[i]->AddVariable( "NVertexes", &NVertexes );
|
| 135 |
}
|
| 136 |
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 |
|
| 145 |
assert(fReaderEndcap);
|
| 146 |
assert(fReaderBarrel);
|
| 147 |
|
| 148 |
}
|
| 149 |
|
| 150 |
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 |
|
| 152 |
//initilize the bool value
|
| 153 |
PassMVA=kFALSE;
|
| 154 |
|
| 155 |
Float_t photon_bdt = MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho);
|
| 156 |
|
| 157 |
if (isbarrel) {
|
| 158 |
if(bdt>bdtCutBarrel){
|
| 159 |
PassMVA=kTRUE;
|
| 160 |
}
|
| 161 |
}
|
| 162 |
else {
|
| 163 |
if(bdt>bdtCutEndcap){
|
| 164 |
PassMVA=kTRUE;
|
| 165 |
}
|
| 166 |
}
|
| 167 |
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 |
Float_t MVATools::GetMVAbdtValue(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho) {
|
| 212 |
|
| 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 |
//mva varialbes v1 and v2
|
| 234 |
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 |
|
| 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 |
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 |
|
| 267 |
RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
|
| 268 |
NVertexes=vtxCol->GetEntries();
|
| 269 |
|
| 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 |
/* printf("HoE: %f\n",HoE);
|
| 288 |
printf("covIEtaIEta: %f\n",covIEtaIEta);
|
| 289 |
printf("tIso1abs: %f\n",tIso1abs);
|
| 290 |
printf("tIso3abs: %f\n",tIso3abs);
|
| 291 |
printf("tIso2abs: %f\n",tIso2abs);
|
| 292 |
|
| 293 |
printf("absIsoEcal: %f\n",absIsoEcal);
|
| 294 |
printf("absIsoHcal: %f\n",absIsoHcal);
|
| 295 |
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 |
if (!isbarrel) {
|
| 313 |
printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
|
| 314 |
}*/
|
| 315 |
|
| 316 |
return bdt;
|
| 317 |
}
|
