ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVATools.cc
Revision: 1.17
Committed: Thu Aug 2 12:30:55 2012 UTC (12 years, 9 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.16: +722 -451 lines
Log Message: 
updated MVATools

File Contents

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