ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVATools.cc
Revision: 1.22
Committed: Sat Nov 9 15:06:27 2013 UTC (11 years, 5 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.21: +43 -6 lines
Error occurred while calculating annotation data.
Log Message:
*** empty log message ***

File Contents

# Content
1 // $Id: MVATools.cc,v 1.21 2013/07/30 21:08:50 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/PFCandidateCol.h"
8 #include "MitAna/DataTree/interface/StableData.h"
9 #include <TMath.h>
10 #include <TFile.h>
11 #include <TRandom3.h>
12 #include <TSystem.h>
13 #include "TMVA/Tools.h"//MVA
14 #include "TMVA/Reader.h"//MVA
15
16
17 ClassImp(mithep::MVATools)
18
19 using namespace mithep;
20
21
22 //--------------------------------------------------------------------------------------------------
23 MVATools::MVATools():
24 fReaderEndcap(0),
25 fReaderBarrel(0),
26
27 fMVAType (MVATools::kNone)
28
29 // ------------------------------------------------------------------------------
30 // fab: These guys should all go away....
31 //MVA Variables v4
32 // HoE(0),
33 // covIEtaIEta(0),
34 // tIso1abs(0),
35 // tIso3abs(0),
36 // tIso2abs(0),
37 // R9(0),
38
39 // absIsoEcal(0),
40 // absIsoHcal(0),
41 // RelEMax(0),
42 // RelETop(0),
43 // RelEBottom(0),
44 // RelELeft(0),
45 // RelERight(0),
46 // RelE2x5Max(0),
47 // RelE2x5Top(0),
48 // RelE2x5Bottom(0),
49 // RelE2x5Left(0),
50 // RelE2x5Right(0),
51 // RelE5x5(0),
52
53 // EtaWidth(0),
54 // PhiWidth(0),
55 // CoviEtaiPhi(0),
56 // CoviPhiiPhi(0),
57
58 // NVertexes(0),
59 // RelPreshowerEnergy(0),
60
61 // //MVA v2 and v1
62 // RelIsoEcal(0),
63 // RelIsoHcal(0),
64 // tIso1(0),
65 // tIso3(0),
66 // tIso2(0),
67 // ScEta(0.)
68
69 {
70 // Constructor.
71 }
72
73 //--------------------------------------------------------------------------------------------------
74 void MVATools::InitializeMVA(MVATools::IdMVAType type) {
75
76 fMVAType = type;
77
78 if (fReaderEndcap) delete fReaderEndcap;
79 if (fReaderBarrel) delete fReaderBarrel;
80
81 // no MVA needed if none requested
82 if( type == kNone ) {
83 return;
84 }
85
86 fReaderEndcap = new TMVA::Reader( "!Color:!Silent:Error" );
87 fReaderBarrel = new TMVA::Reader( "!Color:!Silent:Error" );
88
89 TString BarrelWeights;
90 TString EndcapWeights;
91
92 varNames.resize(0);
93
94 switch (type) {
95
96 case k2011IdMVA_HZg:
97
98 EndcapWeights = (gSystem->Getenv("CMSSW_BASE")+
99 TString("/src/MitPhysics/data/")+
100 TString("PhotonId_lowPt_EE_BDTG.")+
101 TString("weights.xml"));
102 BarrelWeights = (gSystem->Getenv("CMSSW_BASE")+
103 TString("/src/MitPhysics/data/")+
104 TString("PhotonId_lowPt_EB_BDTG.")+
105 TString("weights.xml"));
106
107
108 varNames.push_back("sigieie");
109 varNames.push_back("covieip");
110 varNames.push_back("s4r" );
111 varNames.push_back("r9" );
112 varNames.push_back("sigeta" );
113 varNames.push_back("sigphi" );
114 varNames.push_back("pfgios" );
115 varNames.push_back("pfciso" );
116 varNames.push_back("rho" );
117 varNames.push_back("sceta" );
118 varNames.push_back("rawe" );
119
120 mvaVars.resize(varNames.size());
121
122 std::cout<<" Adding stuff here.... "<<std::endl;
123
124 for( unsigned int iV = 0; iV < mvaVars.size(); ++iV) {
125 mvaVarMapEB.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
126 mvaVarMapEE.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
127 }
128
129 std::cout<<" ... done "<<std::endl;
130
131 break;
132
133 case k2011IdMVA:
134
135 EndcapWeights = (gSystem->Getenv("CMSSW_BASE")+
136 TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
137 TString("Endcap_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
138 TString("weights.xml"));
139 BarrelWeights = (gSystem->Getenv("CMSSW_BASE")+
140 TString("/src/MitPhysics/data/TMVAClassificationPhotonID_")+
141 TString("Barrel_PassPreSel_Variable_10_BDTnCuts2000_BDT.")+
142 TString("weights.xml"));
143
144 // set up the variable names
145 mvaVars.resize(12);
146 varNames.push_back("HoE" );
147 varNames.push_back("covIEtaIEta");
148 varNames.push_back("tIso1abs" );
149 varNames.push_back("tIso3abs" );
150 varNames.push_back("tIso2abs" );
151 varNames.push_back("R9" );
152 varNames.push_back("absIsoEcal" );
153 varNames.push_back("absIsoHcal" );
154 varNames.push_back("NVertexes" );
155 varNames.push_back("ScEta" );
156 varNames.push_back("EtaWidth" );
157 varNames.push_back("PhiWidth" );
158
159 for( unsigned int iV = 0; iV < mvaVars.size(); ++iV) {
160 mvaVarMapEB.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
161 mvaVarMapEE.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
162 }
163
164 break;
165
166 case k2012IdMVA:
167 case k2012IdMVA_globe:
168
169 EndcapWeights = (gSystem->Getenv("CMSSW_BASE")+
170 TString("/src/MitPhysics/data/")+
171 TString("2012ICHEP_PhotonID_Endcap_BDT.")+
172 TString("weights_PSCorr.xml"));
173 BarrelWeights = (gSystem->Getenv("CMSSW_BASE")+
174 TString("/src/MitPhysics/data/")+
175 TString("2012ICHEP_PhotonID_Barrel_BDT.")+
176 TString("weights.xml"));
177
178 mvaVars.resize(12);
179 varNames.push_back("ph.r9" );
180 varNames.push_back("ph.sigietaieta" );
181 varNames.push_back("ph.scetawidth" );
182 varNames.push_back("ph.scphiwidth" );
183 varNames.push_back("ph.idmva_CoviEtaiPhi" );
184 varNames.push_back("ph.idmva_s4ratio" );
185 varNames.push_back("ph.idmva_GammaIso" );
186 varNames.push_back("ph.idmva_ChargedIso_selvtx" );
187 varNames.push_back("ph.idmva_ChargedIso_worstvtx" );
188 varNames.push_back("ph.sceta" );
189 varNames.push_back("rho" );
190 varNames.push_back("ph.idmva_PsEffWidthSigmaRR" );
191
192 for( unsigned int iV = 0; iV < mvaVars.size() - 1; ++iV) {
193 mvaVarMapEB.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
194 mvaVarMapEE.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
195 }
196
197 // pre-shower only used for Endcaps
198 mvaVarMapEE.insert( std::pair<std::string,unsigned int> ( varNames[mvaVars.size() - 1], mvaVars.size() - 1) );
199
200 break;
201
202 case k2013FinalIdMVA_8TeV:
203
204 EndcapWeights = (gSystem->Getenv("CMSSW_BASE")+
205 TString("/src/MitPhysics/data/2013FinalPaper_PhotonID_Weight/8TeV/")+
206 TString("2013FinalPaper_PhotonID_Endcap_BDT_TrainRangePT15_8TeV.")+
207 TString("weights.xml"));
208 BarrelWeights = (gSystem->Getenv("CMSSW_BASE")+
209 TString("/src/MitPhysics/data/2013FinalPaper_PhotonID_Weight/8TeV/")+
210 TString("2013FinalPaper_PhotonID_Barrel_BDT_TrainRangePT15_8TeV.")+
211 TString("weights.xml"));
212
213 mvaVars.resize(13);
214 varNames.push_back("ph.scrawe" );
215 varNames.push_back("ph.r9" );
216 varNames.push_back("ph.sigietaieta" );
217 varNames.push_back("ph.scetawidth" );
218 varNames.push_back("ph.scphiwidth" );
219 varNames.push_back("ph.idmva_CoviEtaiPhi" );
220 varNames.push_back("ph.idmva_s4ratio" );
221 varNames.push_back("ph.idmva_GammaIso" );
222 varNames.push_back("ph.idmva_ChargedIso_selvtx" );
223 varNames.push_back("ph.idmva_ChargedIso_worstvtx" );
224 varNames.push_back("ph.sceta" );
225 varNames.push_back("rho" );
226 varNames.push_back("ph.idmva_PsEffWidthSigmaRR" );
227
228 for( unsigned int iV = 0; iV < mvaVars.size() - 1; ++iV) {
229 mvaVarMapEB.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
230 mvaVarMapEE.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
231 }
232
233 // pre-shower only used for Endcaps
234 mvaVarMapEE.insert( std::pair<std::string,unsigned int> ( varNames[mvaVars.size() - 1], mvaVars.size() - 1) );
235
236 break;
237
238 case k2013FinalIdMVA_7TeV:
239
240 EndcapWeights = (gSystem->Getenv("CMSSW_BASE")+
241 TString("/src/MitPhysics/data/2013FinalPaper_PhotonID_Weight/7TeV/")+
242 TString("2013FinalPaper_PhotonID_Endcap_BDT_TrainRangePT15_7TeV.")+
243 TString("weights.xml"));
244 BarrelWeights = (gSystem->Getenv("CMSSW_BASE")+
245 TString("/src/MitPhysics/data/2013FinalPaper_PhotonID_Weight/7TeV/")+
246 TString("2013FinalPaper_PhotonID_Barrel_BDT_TrainRangePT15_7TeV.")+
247 TString("weights.xml"));
248
249 mvaVars.resize(13);
250 varNames.push_back("ph.scrawe" );
251 varNames.push_back("ph.r9" );
252 varNames.push_back("ph.sigietaieta" );
253 varNames.push_back("ph.scetawidth" );
254 varNames.push_back("ph.scphiwidth" );
255 varNames.push_back("ph.idmva_CoviEtaiPhi" );
256 varNames.push_back("ph.idmva_s4ratio" );
257 varNames.push_back("ph.idmva_GammaIso" );
258 varNames.push_back("ph.idmva_ChargedIso_selvtx" );
259 varNames.push_back("ph.idmva_ChargedIso_worstvtx" );
260 varNames.push_back("ph.sceta" );
261 varNames.push_back("rho" );
262 varNames.push_back("ph.idmva_PsEffWidthSigmaRR" );
263
264 for( unsigned int iV = 0; iV < mvaVars.size() - 1; ++iV) {
265 mvaVarMapEB.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
266 mvaVarMapEE.insert( std::pair<std::string,unsigned int>(varNames[iV], iV) );
267 }
268
269 // pre-shower only used for Endcaps
270 mvaVarMapEE.insert( std::pair<std::string,unsigned int> ( varNames[mvaVars.size() - 1], mvaVars.size() - 1) );
271
272 break;
273
274 default:
275 // no variables... better never called..
276 std::cerr<<" MVATools: Trying to initialize with unknown type."<<std::endl;
277 break;
278 }
279
280
281 // looping over both maps and adding Vars to BDT readers
282 for( unsigned int iV = 0; iV < varNames.size(); ++iV ){
283 std::map<std::string,unsigned int>::const_iterator it = mvaVarMapEB.find(varNames[iV]);
284 if ( it != mvaVarMapEB.end() )
285 fReaderBarrel -> AddVariable( (it->first).c_str(), &(mvaVars[it->second]));
286 it = mvaVarMapEE.find(varNames[iV]);
287 if ( it != mvaVarMapEE.end() )
288 fReaderEndcap -> AddVariable( (it->first).c_str(), &(mvaVars[it->second]));
289 }
290
291 fReaderEndcap->BookMVA("BDT method",EndcapWeights);
292 fReaderBarrel->BookMVA("BDT method",BarrelWeights);
293
294 assert(fReaderEndcap);
295 assert(fReaderBarrel);
296
297 }
298
299 // //--------------------------------------------------------------------------------------------------
300 // void MVATools::InitializeMVA(int VariableType, TString EndcapWeights,TString BarrelWeights) {
301
302 // if (fReaderEndcap) delete fReaderEndcap;
303 // if (fReaderBarrel) delete fReaderBarrel;
304
305 // fReaderEndcap = new TMVA::Reader( "!Color:!Silent:Error" );
306 // fReaderBarrel = new TMVA::Reader( "!Color:!Silent:Error" );
307
308 // TMVA::Reader *readers[2];
309 // readers[0] = fReaderEndcap;
310 // readers[1] = fReaderBarrel;
311
312 // for (UInt_t i=0; i<2; ++i) {
313
314 // if(VariableType==0||VariableType==1||VariableType==2){
315 // readers[i]->AddVariable( "HoE", &HoE );
316 // readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
317 // readers[i]->AddVariable( "tIso1", &tIso1 );
318 // readers[i]->AddVariable( "tIso3", &tIso3 );
319 // readers[i]->AddVariable( "tIso2", &tIso2 );
320 // readers[i]->AddVariable( "R9", &R9 );
321 // }
322
323 // if(VariableType==3||VariableType==4){
324 // readers[i]->AddVariable( "HoE", &HoE );
325 // readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
326 // readers[i]->AddVariable( "tIso1abs", &tIso1abs );
327 // readers[i]->AddVariable( "tIso3abs", &tIso3abs );
328 // readers[i]->AddVariable( "tIso2abs", &tIso2abs );
329 // readers[i]->AddVariable( "R9", &R9 );
330 // }
331
332 // if(VariableType==1||VariableType==2){
333 // readers[i]->AddVariable( "RelIsoEcal", &RelIsoEcal );
334 // readers[i]->AddVariable( "RelIsoHcal", &RelIsoHcal );
335 // readers[i]->AddVariable( "RelEMax", &RelEMax );
336 // readers[i]->AddVariable( "RelETop", &RelETop );
337 // readers[i]->AddVariable( "RelEBottom", &RelEBottom );
338 // readers[i]->AddVariable( "RelELeft", &RelELeft );
339 // readers[i]->AddVariable( "RelERight", &RelERight );
340 // readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
341 // readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
342 // readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
343 // readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
344 // readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
345 // readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
346 // }
347
348 // if(VariableType==3||VariableType==4){
349 // readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
350 // readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
351 // readers[i]->AddVariable( "RelEMax", &RelEMax );
352 // readers[i]->AddVariable( "RelETop", &RelEBottom);
353 // readers[i]->AddVariable( "RelEBottom", &RelEBottom );
354 // readers[i]->AddVariable( "RelELeft", &RelELeft );
355 // readers[i]->AddVariable( "RelERight", &RelERight );
356 // readers[i]->AddVariable( "RelE2x5Max", &RelE2x5Max );
357 // readers[i]->AddVariable( "RelE2x5Top", &RelE2x5Top );
358 // readers[i]->AddVariable( "RelE2x5Bottom", &RelE2x5Bottom );
359 // readers[i]->AddVariable( "RelE2x5Left", &RelE2x5Left );
360 // readers[i]->AddVariable( "RelE2x5Right", &RelE2x5Right );
361 // readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
362 // }
363
364 // if(VariableType==2||VariableType==3||VariableType==4){
365 // readers[i]->AddVariable( "EtaWidth", &EtaWidth );
366 // readers[i]->AddVariable( "PhiWidth", &PhiWidth );
367 // readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
368 // readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
369 // if(VariableType==4){
370 // readers[i]->AddVariable( "NVertexes", &NVertexes );
371 // }
372 // if(i==0){
373 // readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
374 // }
375 // }
376
377 // if(VariableType==6){
378 // readers[i]->AddVariable( "HoE", &HoE );
379 // readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
380 // readers[i]->AddVariable( "tIso1abs", &tIso1abs );
381 // readers[i]->AddVariable( "tIso3abs", &tIso3abs );
382 // readers[i]->AddVariable( "tIso2abs", &tIso2abs );
383 // readers[i]->AddVariable( "R9", &R9 );
384 // readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
385 // readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
386 // readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
387 // readers[i]->AddVariable( "EtaWidth", &EtaWidth );
388 // readers[i]->AddVariable( "PhiWidth", &PhiWidth );
389 // readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
390 // readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
391 // readers[i]->AddVariable( "NVertexes", &NVertexes );
392 // if(i==0){
393 // readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
394 // }
395 // }
396
397 // if(VariableType==7){
398 // readers[i]->AddVariable( "HoE", &HoE );
399 // readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
400 // readers[i]->AddVariable( "tIso1abs", &tIso1abs );
401 // readers[i]->AddVariable( "tIso3abs", &tIso3abs );
402 // readers[i]->AddVariable( "tIso2abs", &tIso2abs );
403 // readers[i]->AddVariable( "R9", &R9 );
404 // readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
405 // readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
406 // readers[i]->AddVariable( "NVertexes", &NVertexes );
407 // readers[i]->AddVariable( "ScEta", &ScEta );
408 // }
409
410 // if(VariableType==10){
411 // readers[i]->AddVariable( "HoE", &HoE );
412 // readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
413 // readers[i]->AddVariable( "tIso1abs", &tIso1abs );
414 // readers[i]->AddVariable( "tIso3abs", &tIso3abs );
415 // readers[i]->AddVariable( "tIso2abs", &tIso2abs );
416 // readers[i]->AddVariable( "R9", &R9 );
417 // readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
418 // readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
419 // readers[i]->AddVariable( "NVertexes", &NVertexes );
420 // readers[i]->AddVariable( "ScEta", &ScEta );
421 // readers[i]->AddVariable( "EtaWidth", &EtaWidth );
422 // readers[i]->AddVariable( "PhiWidth", &PhiWidth );
423 // }
424
425 // if(VariableType==1201){
426 // /*readers[i]->AddVariable( "myphoton_pfchargedisogood03", &myphoton_pfchargedisogood03);
427 // readers[i]->AddVariable( "myphoton_pfchargedisobad03", &myphoton_pfchargedisobad03);
428 // readers[i]->AddVariable( "myphoton_pfphotoniso03", &myphoton_pfphotoniso03 );
429 // readers[i]->AddVariable( "myphoton_sieie", &myphoton_sieie );
430 // readers[i]->AddVariable( "myphoton_sieip", &myphoton_sieip );
431 // readers[i]->AddVariable( "myphoton_etawidth", &myphoton_etawidth );
432 // readers[i]->AddVariable( "myphoton_phiwidth", &myphoton_phiwidth );
433 // readers[i]->AddVariable( "myphoton_r9", &myphoton_r9 );
434 // readers[i]->AddVariable( "myphoton_s4ratio", &myphoton_s4ratio );
435 // readers[i]->AddVariable( "myphoton_SCeta", &myphoton_SCeta );
436 // readers[i]->AddVariable( "event_rho", &event_rho );
437 // if(i==0){
438 // readers[i]->AddVariable( "myphoton_ESEffSigmaRR", &myphoton_ESEffSigmaRR);
439 // }*/
440 // readers[i]->AddVariable( "ph.r9", &myphoton_r9 );
441 // readers[i]->AddVariable( "ph.sigietaieta", &myphoton_sieie );
442 // readers[i]->AddVariable( "ph.scetawidth", &myphoton_etawidth );
443 // readers[i]->AddVariable( "ph.scphiwidth", &myphoton_phiwidth );
444 // readers[i]->AddVariable( "ph.idmva_CoviEtaiPhi", &myphoton_sieip );
445 // readers[i]->AddVariable( "ph.idmva_s4ratio", &myphoton_s4ratio );
446 // readers[i]->AddVariable( "ph.idmva_GammaIso", &myphoton_pfphotoniso03 );
447 // readers[i]->AddVariable( "ph.idmva_ChargedIso_selvtx", &myphoton_pfchargedisogood03);
448 // readers[i]->AddVariable( "ph.idmva_ChargedIso_worstvtx", &myphoton_pfchargedisobad03);
449 // readers[i]->AddVariable( "ph.sceta", &myphoton_SCeta );
450 // readers[i]->AddVariable( "rho", &event_rho );
451 // if(i==0){
452 // //readers[i]->AddVariable( "1.00023*ph.idmva_PsEffWidthSigmaRR + 0.0913", &myphoton_ESEffSigmaRR);
453 // readers[i]->AddVariable( "ph.idmva_PsEffWidthSigmaRR", &myphoton_ESEffSigmaRR);
454 // }
455 // }
456
457 // }
458
459 // fReaderEndcap->BookMVA("BDT method",EndcapWeights);
460 // fReaderBarrel->BookMVA("BDT method",BarrelWeights);
461
462 // assert(fReaderEndcap);
463 // assert(fReaderBarrel);
464
465 // }
466
467 // *** REMOVED THIS COMPLETELY. If a module wants to cut on BDT, it should 1,) compute the BDT value (using GetMVAbdtValue(...) ) and then make the cut itself...
468
469 // 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) {
470
471 // //initilize the bool value
472 // PassMVA=kFALSE;
473
474 // Float_t photon_bdt = MVATools::GetMVAbdtValue_2011(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
475
476 // if (isbarrel) {
477 // if(photon_bdt>bdtCutBarrel){
478 // PassMVA=kTRUE;
479 // }
480 // }
481 // else {
482 // if(photon_bdt>bdtCutEndcap){
483 // PassMVA=kTRUE;
484 // }
485 // }
486 // return PassMVA;
487 // }
488
489 // //---------------------------------------------------------------------------------
490 // Int_t MVATools::PassElectronVetoInt(const Photon* p, const ElectronCol* els) {
491
492 // // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
493 // float cic4_allcuts_temp_sublead[] = {
494 // 3.8, 2.2, 1.77, 1.29,
495 // 11.7, 3.4, 3.9, 1.84,
496 // 3.5, 2.2, 2.3, 1.45,
497 // 0.0106, 0.0097, 0.028, 0.027,
498 // 0.082, 0.062, 0.065, 0.048,
499 // 0.94, 0.36, 0.94, 0.32,
500 // 1., 0.062, 0.97, 0.97,
501 // 1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
502
503 // //initilize the bool value
504 // PassElecVetoInt=0;
505
506 // dRTrack = PhotonTools::ElectronVetoCiC(p, els);
507
508 // ScEta_MVA=p->SCluster()->Eta();
509
510 // isbarrel = (fabs(ScEta_MVA)<1.4442);
511
512 // R9 = p->R9();
513 // //R9 = p->E33()/p->SCluster()->RawEnergy();
514
515 // // check which category it is ...
516 // _tCat = 1;
517 // if ( !isbarrel ) _tCat = 3;
518 // if ( R9 < 0.94 ) _tCat++;
519
520 // //Electron Veto
521 // if(dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4]){
522 // PassElecVetoInt=1;
523 // }
524
525 // return PassElecVetoInt;
526
527 // }
528
529 //--------------------------------------------------------------------------------------------------
530
531
532 Double_t MVATools::GetMVAbdtValue(const Photon* p, const Vertex* vtx, const TrackCol* trackCol, const VertexCol* vtxCol, Double_t _tRho, const PFCandidateCol *fPFCands, const ElectronCol* els, Bool_t applyElectronVeto) {
533
534 // if there's no reader, or the type is kNone, return the default values of -99.
535 //if( ( !fReaderBarrel || !fReaderEndcap ) ) return -199.;
536 if( ( !fReaderBarrel || !fReaderEndcap ) || fMVAType == kNone ) return -99.;
537
538 // we compute the variable names... make sure no confilcts when adding new variables...
539
540 // check if it's a Barrel or EE photon
541 bool isBarrel = ( p->SCluster()->AbsEta() < 1.5 );
542
543 std::map<std::string,unsigned int>* theVars = ( isBarrel ? &mvaVarMapEB : &mvaVarMapEE );
544
545 // loop over all the variables in the map... and keep count (to make sure we have filled all variables)
546 unsigned int varCounter = 0;
547 for( std::map<std::string,unsigned int>::const_iterator iV = theVars->begin(); iV != theVars->end(); ++iV ) {
548
549 TString theVarName = TString(iV->first);
550 float* theVarValue = &(mvaVars[iV->second]); // pointer to the variable...
551
552 if(
553 !theVarName.CompareTo("HoE")
554 ) {
555 (*theVarValue) = p->HadOverEm();
556 varCounter++;
557 } else if (
558 !theVarName.CompareTo("covIEtaIEta") || !theVarName.CompareTo("ph.sigietaieta") || !theVarName.CompareTo("sigieie")
559 ) {
560 (*theVarValue) = p->CoviEtaiEta();
561 varCounter++;
562 } else if (
563 !theVarName.CompareTo("R9") || !theVarName.CompareTo("ph.r9") || !theVarName.CompareTo("r9")
564 ) {
565 (*theVarValue) = p->R9();
566 varCounter++;
567 } else if (
568 !theVarName.CompareTo("ScEta") || !theVarName.CompareTo("ph.sceta") || !theVarName.CompareTo("sceta")
569 ) {
570 (*theVarValue) = p->SCluster()->Eta();
571 varCounter++;
572 } else if (
573 !theVarName.CompareTo("rho")
574 ) {
575 (*theVarValue) = _tRho;
576 varCounter++;
577 } else if (
578 !theVarName.CompareTo("tIso1abs")
579 ) {
580 double _ecalIso3 = p->EcalRecHitIsoDr03();
581 double _hcalIso4 = p->HcalTowerSumEtDr04();
582 double _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?
583 (*theVarValue) = _ecalIso3+_hcalIso4+_trackIso1 - 0.17*_tRho;
584 varCounter++;
585 } else if (
586 !theVarName.CompareTo("tIso2abs")
587 ) {
588 double _ecalIso4 = p->EcalRecHitIsoDr04();
589 double _hcalIso4 = p->HcalTowerSumEtDr04();
590 unsigned int wVtxInd = 0;
591 double _trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
592 (*theVarValue) = _ecalIso4+_hcalIso4+_trackIso2 - 0.52*_tRho;
593 varCounter++;
594 } else if (
595 !theVarName.CompareTo("tIso3abs")
596 ) {
597 (*theVarValue) = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );
598 varCounter++;
599 } else if (
600 !theVarName.CompareTo("absIsoEcal")
601 ) {
602 double _ecalIso3 = p->EcalRecHitIsoDr03();
603 (*theVarValue) = (_ecalIso3-0.17*_tRho);
604 varCounter++;
605 } else if (
606 !theVarName.CompareTo("absIsoHcal")
607 ) {
608 double _hcalIso4 = p->HcalTowerSumEtDr04();
609 (*theVarValue) = (_hcalIso4-0.17*_tRho);
610 varCounter++;
611 } else if (
612 !theVarName.CompareTo("NVertexes")
613 ) {
614
615 (*theVarValue) = vtxCol->GetEntries();
616 varCounter++;
617 } else if (
618 !theVarName.CompareTo("EtaWidth") || !theVarName.CompareTo("ph.scetawidth") || !theVarName.CompareTo("sigeta")
619 ) {
620 (*theVarValue) = p->EtaWidth();
621 varCounter++;
622 } else if (
623 !theVarName.CompareTo("PhiWidth") || !theVarName.CompareTo("ph.scphiwidth") || !theVarName.CompareTo("sigphi")
624 ) {
625 (*theVarValue) = p->PhiWidth();
626 varCounter++;
627
628 } else if (
629 !theVarName.CompareTo("ph.idmva_CoviEtaiPhi") || !theVarName.CompareTo("covieip")
630 ) {
631 (*theVarValue) = p->SCluster()->Seed()->CoviEtaiPhi();
632 varCounter++;
633 } else if (
634 !theVarName.CompareTo("ph.idmva_s4ratio") || !theVarName.CompareTo("s4r")
635 ) {
636 (*theVarValue) = p->S4Ratio();
637 //(*theVarValue) = p->SCluster()->Seed()->E2x2()/p->E55();
638 varCounter++;
639 } else if (
640 !theVarName.CompareTo("ph.idmva_GammaIso") || !theVarName.CompareTo("pfgiso")
641 ) {
642 (*theVarValue) = IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
643 varCounter++;
644 } else if (
645 !theVarName.CompareTo("ph.idmva_ChargedIso_selvtx") || !theVarName.CompareTo("pfciso")
646 ) {
647 (*theVarValue) = IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
648 varCounter++;
649 } else if (
650 !theVarName.CompareTo("ph.idmva_ChargedIso_worstvtx")
651 ) {
652 unsigned int wVtxInd = 0;
653 (*theVarValue) = IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol);
654 varCounter++;
655 } else if (
656 !theVarName.CompareTo("ph.idmva_PsEffWidthSigmaRR")
657 ) {
658 (*theVarValue) = p->EffSigmaRR();
659 varCounter++;
660 } else if (
661 !theVarName.CompareTo("rawe")|| !theVarName.CompareTo("ph.scrawe")
662 ) {
663 (*theVarValue) = p->SCluster()->RawEnergy();
664 varCounter++;
665 } else {
666 // a variable is not know... copmplain!
667 std::cerr<<" ERROR: MVA Evaluation called with unknown variable name >"<<theVarName<<">."<<std::endl;
668 }
669 }
670
671 // now all the variables should be filled... check!
672 if( varCounter != theVars->size() )
673 std::cerr<<" ERROR: MVA Evaludation called and not all variables are filled."<<std::endl;
674
675 // we're ready to compute the MVA value
676 TMVA::Reader* reader = NULL;
677 if (isBarrel)
678 reader = fReaderBarrel;
679 else
680 reader = fReaderEndcap;
681
682 assert(reader);
683
684 return (reader->EvaluateMVA("BDT method"));
685 }
686
687 // Float_t MVATools::GetMVAbdtValue_2012_globe(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho, const PFCandidateCol *fPFCands,const ElectronCol* els,Bool_t applyElectronVeto) {
688
689 // //get the variables used to compute MVA variables
690 // ecalIso3 = p->EcalRecHitIsoDr03();
691 // ecalIso4 = p->EcalRecHitIsoDr04();
692 // hcalIso4 = p->HcalTowerSumEtDr04();
693
694 // wVtxInd = 0;
695
696 // 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?
697
698 // // track iso worst vtx
699 // trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
700
701 // combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
702 // combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
703
704 // RawEnergy = p->SCluster()->RawEnergy();
705
706 // ScEta = p->SCluster()->Eta();
707
708 // //mva varialbes v1 and v2
709 // tIso1 = (combIso1) *50./p->Et();
710 // tIso3 = (trackIso1)*50./p->Et();
711 // tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
712 // RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
713 // RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
714
715 // //compute mva variables for v3
716 // HoE = p->HadOverEm();
717 // covIEtaIEta = p->CoviEtaiEta();
718 // tIso1abs = combIso1;
719 // tIso3abs = trackIso1;
720 // tIso2abs = combIso2;
721 // R9 = p->R9();
722
723 // absIsoEcal=(ecalIso3-0.17*_tRho);
724 // absIsoHcal=(hcalIso4-0.17*_tRho);
725 // RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
726 // RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
727 // RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
728 // RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
729 // RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
730 // RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
731 // RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
732 // RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
733 // RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
734 // RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
735 // RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
736
737 // EtaWidth=p->EtaWidth();
738 // PhiWidth=p->PhiWidth();
739 // CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
740 // CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
741
742 // RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
743 // NVertexes=vtxCol->GetEntries();
744
745 // //spectator variables
746 // Pt_MVA=p->Pt();
747 // ScEta_MVA=p->SCluster()->Eta();
748
749 // //
750
751 // isbarrel = (fabs(ScEta_MVA)<1.4442);
752
753 // //variable 1201
754 // myphoton_pfchargedisogood03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
755 // myphoton_pfchargedisobad03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol);
756 // myphoton_pfphotoniso03=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
757 // myphoton_sieie=covIEtaIEta;
758 // myphoton_sieip=CoviEtaiPhi;
759 // myphoton_etawidth=EtaWidth;
760 // myphoton_phiwidth=PhiWidth;
761 // myphoton_r9=R9;
762 // myphoton_s4ratio=p->S4Ratio();
763 // myphoton_SCeta=ScEta_MVA;
764 // event_rho= _tRho;
765
766 // myphoton_ESEffSigmaRR=-99;
767
768 // if(!isbarrel){
769 // myphoton_ESEffSigmaRR=p->EffSigmaRR();
770 // }
771
772 // if (isbarrel) {
773 // reader = fReaderBarrel;
774 // }
775 // else {
776 // reader = fReaderEndcap;
777 // }
778
779 // assert(reader);
780
781 // double bdt = reader->EvaluateMVA("BDT method");
782
783 // /* printf("HoE: %f\n",HoE);
784 // printf("covIEtaIEta: %f\n",covIEtaIEta);
785 // printf("tIso1abs: %f\n",tIso1abs);
786 // printf("tIso3abs: %f\n",tIso3abs);
787 // printf("tIso2abs: %f\n",tIso2abs);
788
789 // printf("absIsoEcal: %f\n",absIsoEcal);
790 // printf("absIsoHcal: %f\n",absIsoHcal);
791 // printf("RelEMax: %f\n",RelEMax);
792 // printf("RelETop: %f\n",RelETop);
793 // printf("RelEBottom: %f\n",RelEBottom);
794 // printf("RelELeft: %f\n",RelELeft);
795 // printf("RelERight: %f\n",RelERight);
796 // printf("RelE2x5Max: %f\n",RelE2x5Max);
797 // printf("RelE2x5Top: %f\n",RelE2x5Top);
798 // printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
799 // printf("RelE2x5Left: %f\n",RelE2x5Left);
800 // printf("RelE2x5Right;: %f\n",RelE2x5Right);
801 // printf("RelE5x5: %f\n",RelE5x5);
802
803 // printf("EtaWidth: %f\n",EtaWidth);
804 // printf("PhiWidth: %f\n",PhiWidth);
805 // printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
806 // printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
807
808 // if (!isbarrel) {
809 // printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
810 // }*/
811
812 // return bdt;
813 // }
814
815 // Float_t MVATools::GetMVAbdtValue_2011(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els,Bool_t applyElectronVeto) {
816
817 // //get the variables used to compute MVA variables
818 // ecalIso3 = p->EcalRecHitIsoDr03();
819 // ecalIso4 = p->EcalRecHitIsoDr04();
820 // hcalIso4 = p->HcalTowerSumEtDr04();
821
822 // wVtxInd = 0;
823
824 // 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?
825
826 // // track iso worst vtx
827 // trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
828
829 // combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
830 // combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
831
832 // RawEnergy = p->SCluster()->RawEnergy();
833
834 // ScEta = p->SCluster()->Eta();
835
836 // //mva varialbes v1 and v2
837 // tIso1 = (combIso1) *50./p->Et();
838 // tIso3 = (trackIso1)*50./p->Et();
839 // tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
840 // RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
841 // RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
842
843 // //compute mva variables for v3
844 // HoE = p->HadOverEm();
845 // covIEtaIEta = p->CoviEtaiEta();
846 // tIso1abs = combIso1;
847 // tIso3abs = trackIso1;
848 // tIso2abs = combIso2;
849 // R9 = p->R9();
850
851 // absIsoEcal=(ecalIso3-0.17*_tRho);
852 // absIsoHcal=(hcalIso4-0.17*_tRho);
853 // RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
854 // RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
855 // RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
856 // RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
857 // RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
858 // RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
859 // RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
860 // RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
861 // RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
862 // RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
863 // RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
864
865 // EtaWidth=p->EtaWidth();
866 // PhiWidth=p->PhiWidth();
867 // CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
868 // CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
869
870 // RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
871 // NVertexes=vtxCol->GetEntries();
872
873 // //spectator variables
874 // Pt_MVA=p->Pt();
875 // ScEta_MVA=p->SCluster()->Eta();
876
877 // //
878
879 // isbarrel = (fabs(ScEta_MVA)<1.4442);
880
881 // if (isbarrel) {
882 // reader = fReaderBarrel;
883 // }
884 // else {
885 // reader = fReaderEndcap;
886 // }
887
888 // assert(reader);
889
890 // double bdt = reader->EvaluateMVA("BDT method");
891
892 // /* printf("HoE: %f\n",HoE);
893 // printf("covIEtaIEta: %f\n",covIEtaIEta);
894 // printf("tIso1abs: %f\n",tIso1abs);
895 // printf("tIso3abs: %f\n",tIso3abs);
896 // printf("tIso2abs: %f\n",tIso2abs);
897
898 // printf("absIsoEcal: %f\n",absIsoEcal);
899 // printf("absIsoHcal: %f\n",absIsoHcal);
900 // printf("RelEMax: %f\n",RelEMax);
901 // printf("RelETop: %f\n",RelETop);
902 // printf("RelEBottom: %f\n",RelEBottom);
903 // printf("RelELeft: %f\n",RelELeft);
904 // printf("RelERight: %f\n",RelERight);
905 // printf("RelE2x5Max: %f\n",RelE2x5Max);
906 // printf("RelE2x5Top: %f\n",RelE2x5Top);
907 // printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
908 // printf("RelE2x5Left: %f\n",RelE2x5Left);
909 // printf("RelE2x5Right;: %f\n",RelE2x5Right);
910 // printf("RelE5x5: %f\n",RelE5x5);
911
912 // printf("EtaWidth: %f\n",EtaWidth);
913 // printf("PhiWidth: %f\n",PhiWidth);
914 // printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
915 // printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
916
917 // if (!isbarrel) {
918 // printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
919 // }*/
920
921 // return bdt;
922 // }