ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MuonIDMVA.cc
Revision: 1.9
Committed: Tue Apr 24 12:48:30 2012 UTC (13 years ago) by sixie
Content type: text/plain
Branch: MAIN
Changes since 1.8: +184 -20 lines
Log Message:
update electron and muon mva interfaces. add Andrews's muon ID mva.

File Contents

# User Rev Content
1 sixie 1.1 #include "MitPhysics/Utils/interface/MuonIDMVA.h"
2     #include "MitPhysics/Utils/interface/MuonTools.h"
3     #include "MitPhysics/Utils/interface/IsolationTools.h"
4     #include "MitAna/DataTree/interface/StableData.h"
5     #include <TFile.h>
6     #include <TRandom3.h>
7     #include "TMVA/Tools.h"
8     #include "TMVA/Reader.h"
9    
10    
11     ClassImp(mithep::MuonIDMVA)
12    
13     using namespace mithep;
14    
15     //--------------------------------------------------------------------------------------------------
16     MuonIDMVA::MuonIDMVA() :
17     fMethodname("BDTG method"),
18 sixie 1.8 fIsInitialized(kFALSE),
19     fMVAType(MuonIDMVA::kUninitialized),
20     fUseBinnedVersion(kTRUE),
21     fNMVABins(0)
22 sixie 1.1 {
23     }
24    
25    
26     //--------------------------------------------------------------------------------------------------
27     MuonIDMVA::~MuonIDMVA()
28     {
29 sixie 1.8 for(UInt_t i=0; i<fTMVAReader.size(); ++i) {
30 sixie 1.1 if (fTMVAReader[i]) delete fTMVAReader[i];
31     }
32     }
33    
34     //--------------------------------------------------------------------------------------------------
35 sixie 1.8 void MuonIDMVA::Initialize( std::string methodName,
36     std::string weightsfile,
37     MuonIDMVA::MVAType type)
38     {
39    
40     std::vector<std::string> tempWeightFileVector;
41     tempWeightFileVector.push_back(weightsfile);
42     Initialize(methodName,type,kFALSE,tempWeightFileVector);
43     }
44    
45     //--------------------------------------------------------------------------------------------------
46 sixie 1.1 void MuonIDMVA::Initialize( TString methodName,
47 sixie 1.8 TString Subdet0Pt10To20Weights ,
48     TString Subdet1Pt10To20Weights ,
49     TString Subdet2Pt10To20Weights,
50     TString Subdet0Pt20ToInfWeights,
51     TString Subdet1Pt20ToInfWeights,
52     TString Subdet2Pt20ToInfWeights,
53     MuonIDMVA::MVAType type) {
54    
55     std::vector<std::string> tempWeightFileVector;
56     tempWeightFileVector.push_back(std::string(Subdet0Pt10To20Weights.Data()));
57     tempWeightFileVector.push_back(std::string(Subdet1Pt10To20Weights.Data()));
58     tempWeightFileVector.push_back(std::string(Subdet2Pt10To20Weights.Data()));
59     tempWeightFileVector.push_back(std::string(Subdet0Pt20ToInfWeights.Data()));
60     tempWeightFileVector.push_back(std::string(Subdet1Pt20ToInfWeights.Data()));
61     tempWeightFileVector.push_back(std::string(Subdet2Pt20ToInfWeights.Data()));
62     Initialize(std::string(methodName.Data()),type,kTRUE,tempWeightFileVector);
63    
64     }
65    
66    
67     //--------------------------------------------------------------------------------------------------
68     void MuonIDMVA::Initialize( std::string methodName,
69     MuonIDMVA::MVAType type,
70     Bool_t useBinnedVersion,
71     std::vector<std::string> weightsfiles) {
72 sixie 1.1
73 sixie 1.8 //clean up first
74     for (uint i=0;i<fTMVAReader.size(); ++i) {
75     if (fTMVAReader[i]) delete fTMVAReader[i];
76     }
77     fTMVAReader.clear();
78    
79    
80     //initialize
81 sixie 1.1 fIsInitialized = kTRUE;
82     fMethodname = methodName;
83 sixie 1.8 fMVAType = type;
84     fUseBinnedVersion = useBinnedVersion;
85    
86     //Define expected number of bins
87     UInt_t ExpectedNBins = 0;
88     if (!fUseBinnedVersion) {
89     ExpectedNBins = 1;
90     } else if (type == kV2
91 sixie 1.9 || type == kV3
92     || type == kV8
93     || type == kIDIsoCombinedDetIso
94     || type == kIsoRingsV0
95     || type == kIDV0
96     || type == kIDIsoCombinedIsoRingsV0
97     ) {
98 sixie 1.8 ExpectedNBins = 6;
99 sixie 1.9 }
100 sixie 1.8 fNMVABins = ExpectedNBins;
101    
102     //Check number of weight files given
103     if (fNMVABins != weightsfiles.size() ) {
104     std::cout << "Error: Expected Number of bins = " << fNMVABins << " does not equal to weightsfiles.size() = "
105     << weightsfiles.size() << std::endl;
106     assert(fNMVABins == weightsfiles.size());
107     }
108 sixie 1.1
109 sixie 1.8 for(UInt_t i=0; i<fNMVABins; ++i) {
110     TMVA::Reader *tmpTMVAReader = new TMVA::Reader( "!Color:!Silent:Error" );
111     tmpTMVAReader->SetVerbose(kTRUE);
112 sixie 1.1
113 sixie 1.2 if (type == kV2) {
114 sixie 1.8 tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
115     tmpTMVAReader->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
116     tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
117     tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
118     tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
119     tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
120     tmpTMVAReader->AddVariable( "D0", &fMVAVar_MuD0 );
121     tmpTMVAReader->AddVariable( "IP3d", &fMVAVar_MuIP3d );
122     tmpTMVAReader->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
123     tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
124     tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
125 sixie 1.2 }
126    
127 sixie 1.1 if (type == kV3) {
128 sixie 1.8 tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
129     tmpTMVAReader->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
130     tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
131     tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
132     tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
133     tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
134     tmpTMVAReader->AddVariable( "D0", &fMVAVar_MuD0 );
135     tmpTMVAReader->AddVariable( "IP3d", &fMVAVar_MuIP3d );
136     tmpTMVAReader->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
137     tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
138     tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
139     tmpTMVAReader->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
140     tmpTMVAReader->AddVariable( "HadEnergyOverPt", &fMVAVar_MuHadEnergyOverPt );
141 sixie 1.1 if (i==0 || i==2 || i==4) {
142 sixie 1.8 tmpTMVAReader->AddVariable( "HoEnergyOverPt", &fMVAVar_MuHoEnergyOverPt );
143 sixie 1.1 }
144 sixie 1.8 tmpTMVAReader->AddVariable( "EmEnergyOverPt", &fMVAVar_MuEmEnergyOverPt );
145     tmpTMVAReader->AddVariable( "HadS9EnergyOverPt", &fMVAVar_MuHadS9EnergyOverPt );
146 sixie 1.1 if (i==0 || i==2 || i==4) {
147 sixie 1.8 tmpTMVAReader->AddVariable( "HoS9EnergyOverPt", &fMVAVar_MuHoS9EnergyOverPt );
148 sixie 1.1 }
149 sixie 1.8 tmpTMVAReader->AddVariable( "EmS9EnergyOverPt", &fMVAVar_MuEmS9EnergyOverPt );
150 sixie 1.1 }
151 sixie 1.2
152     if (type == kV8) {
153 sixie 1.8 tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
154     tmpTMVAReader->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
155     tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
156     tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
157     tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
158     tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
159     tmpTMVAReader->AddVariable( "D0", &fMVAVar_MuD0 );
160     tmpTMVAReader->AddVariable( "IP3d", &fMVAVar_MuIP3d );
161     tmpTMVAReader->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
162     tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
163     tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
164     tmpTMVAReader->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
165     tmpTMVAReader->AddVariable( "HadEnergyOverPt", &fMVAVar_MuHadEnergyOverPt );
166     tmpTMVAReader->AddVariable( "EmEnergyOverPt", &fMVAVar_MuEmEnergyOverPt );
167     tmpTMVAReader->AddVariable( "HadS9EnergyOverPt", &fMVAVar_MuHadS9EnergyOverPt );
168     tmpTMVAReader->AddVariable( "EmS9EnergyOverPt", &fMVAVar_MuEmS9EnergyOverPt );
169     tmpTMVAReader->AddVariable( "ChargedIso03OverPt", &fMVAVar_MuChargedIso03OverPt );
170     tmpTMVAReader->AddVariable( "NeutralIso03OverPt", &fMVAVar_MuNeutralIso03OverPt );
171     tmpTMVAReader->AddVariable( "ChargedIso04OverPt", &fMVAVar_MuChargedIso04OverPt );
172     tmpTMVAReader->AddVariable( "NeutralIso04OverPt", &fMVAVar_MuNeutralIso04OverPt );
173 sixie 1.2 }
174 sixie 1.1
175 sixie 1.7 if (type == kIDIsoCombinedDetIso) {
176 sixie 1.8 tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
177     tmpTMVAReader->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
178     tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
179     tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
180     tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
181     tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
182     tmpTMVAReader->AddVariable( "D0", &fMVAVar_MuD0 );
183     tmpTMVAReader->AddVariable( "IP3d", &fMVAVar_MuIP3d );
184     tmpTMVAReader->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
185     tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
186     tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
187     tmpTMVAReader->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
188     tmpTMVAReader->AddVariable( "HadEnergyOverPt", &fMVAVar_MuHadEnergyOverPt );
189     tmpTMVAReader->AddVariable( "EmEnergyOverPt", &fMVAVar_MuEmEnergyOverPt );
190     tmpTMVAReader->AddVariable( "HadS9EnergyOverPt", &fMVAVar_MuHadS9EnergyOverPt );
191     tmpTMVAReader->AddVariable( "EmS9EnergyOverPt", &fMVAVar_MuEmS9EnergyOverPt );
192     tmpTMVAReader->AddVariable( "TrkIso03OverPt", &fMVAVar_MuTrkIso03OverPt );
193     tmpTMVAReader->AddVariable( "EMIso03OverPt", &fMVAVar_MuEMIso03OverPt );
194     tmpTMVAReader->AddVariable( "HadIso03OverPt", &fMVAVar_MuHadIso03OverPt );
195     tmpTMVAReader->AddVariable( "TrkIso05OverPt", &fMVAVar_MuTrkIso05OverPt );
196     tmpTMVAReader->AddVariable( "EMIso05OverPt", &fMVAVar_MuEMIso05OverPt );
197     tmpTMVAReader->AddVariable( "HadIso05OverPt", &fMVAVar_MuHadIso05OverPt );
198     }
199    
200 sixie 1.9 if (type == kIDV0) {
201     tmpTMVAReader->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
202     if (i!=4) tmpTMVAReader->AddVariable( "GlobalNchi2",&fMVAVar_MuGlobalNchi2 );
203     if (i!=4) tmpTMVAReader->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
204     tmpTMVAReader->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
205     tmpTMVAReader->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
206     if (i!=5) tmpTMVAReader->AddVariable( "NMatches", &fMVAVar_MuNMatches );
207     tmpTMVAReader->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
208     tmpTMVAReader->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
209     tmpTMVAReader->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
210     tmpTMVAReader->AddVariable( "HadEnergy", &fMVAVar_MuHadEnergy );
211     tmpTMVAReader->AddVariable( "EmEnergy", &fMVAVar_MuEmEnergy );
212     tmpTMVAReader->AddVariable( "HadS9Energy", &fMVAVar_MuHadS9Energy );
213     tmpTMVAReader->AddVariable( "EmS9Energy", &fMVAVar_MuEmS9Energy );
214     }
215    
216 sixie 1.8 if (type == kIsoRingsV0) {
217     tmpTMVAReader->AddVariable( "ChargedIso_DR0p0To0p1", &fMVAVar_ChargedIso_DR0p0To0p1 );
218     tmpTMVAReader->AddVariable( "ChargedIso_DR0p1To0p2", &fMVAVar_ChargedIso_DR0p1To0p2 );
219     tmpTMVAReader->AddVariable( "ChargedIso_DR0p2To0p3", &fMVAVar_ChargedIso_DR0p2To0p3 );
220     tmpTMVAReader->AddVariable( "ChargedIso_DR0p3To0p4", &fMVAVar_ChargedIso_DR0p3To0p4 );
221     tmpTMVAReader->AddVariable( "ChargedIso_DR0p4To0p5", &fMVAVar_ChargedIso_DR0p4To0p5 );
222     tmpTMVAReader->AddVariable( "GammaIso_DR0p0To0p1", &fMVAVar_GammaIso_DR0p0To0p1 );
223     tmpTMVAReader->AddVariable( "GammaIso_DR0p1To0p2", &fMVAVar_GammaIso_DR0p1To0p2 );
224     tmpTMVAReader->AddVariable( "GammaIso_DR0p2To0p3", &fMVAVar_GammaIso_DR0p2To0p3 );
225     tmpTMVAReader->AddVariable( "GammaIso_DR0p3To0p4", &fMVAVar_GammaIso_DR0p3To0p4 );
226     tmpTMVAReader->AddVariable( "GammaIso_DR0p4To0p5", &fMVAVar_GammaIso_DR0p4To0p5 );
227     tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p0To0p1", &fMVAVar_NeutralHadronIso_DR0p0To0p1 );
228     tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p1To0p2", &fMVAVar_NeutralHadronIso_DR0p1To0p2 );
229     tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p2To0p3", &fMVAVar_NeutralHadronIso_DR0p2To0p3 );
230     tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p3To0p4", &fMVAVar_NeutralHadronIso_DR0p3To0p4 );
231     tmpTMVAReader->AddVariable( "NeutralHadronIso_DR0p4To0p5", &fMVAVar_NeutralHadronIso_DR0p4To0p5 );
232 sixie 1.7 }
233    
234 sixie 1.8 tmpTMVAReader->BookMVA(fMethodname , weightsfiles[i] );
235     std::cout << "MVABin " << i << " : MethodName = " << fMethodname
236     << " , type == " << type << " , "
237     << "Load weights file : " << weightsfiles[i]
238     << std::endl;
239     fTMVAReader.push_back(tmpTMVAReader);
240 sixie 1.1
241     }
242 sixie 1.8 std::cout << "Muon ID MVA Completed\n";
243     }
244    
245     //--------------------------------------------------------------------------------------------------
246     UInt_t MuonIDMVA::GetMVABin( double eta, double pt,
247     Bool_t isGlobal, Bool_t isTrackerMuon) const {
248    
249     //Default is to return the first bin
250     uint bin = 0;
251    
252     //return the first bin if not using binned version
253     if (!fUseBinnedVersion) return 0;
254    
255     if (fMVAType == MuonIDMVA::kV2
256     || fMVAType == MuonIDMVA::kV3
257     || fMVAType == MuonIDMVA::kV8
258     || fMVAType == MuonIDMVA::kIDIsoCombinedDetIso) {
259     if (pt < 14.5 && fabs(eta) < 1.5) bin = 0;
260     if (pt < 14.5 && fabs(eta) >= 1.5) bin = 1;
261     if (pt >= 14.5 && pt < 20 && fabs(eta) < 1.5) bin = 2;
262     if (pt >= 14.5 && pt < 20 && fabs(eta) >= 1.5) bin = 3;
263     if (pt >= 20 && fabs(eta) < 1.5) bin = 4;
264     if (pt >= 20 && fabs(eta) >= 1.5) bin = 5;
265     }
266    
267 sixie 1.9 if (fMVAType == MuonIDMVA::kIsoRingsV0 || fMVAType == MuonIDMVA::kIDV0
268     || fMVAType == MuonIDMVA::kIDIsoCombinedIsoRingsV0) {
269 sixie 1.8 if (isGlobal && isTrackerMuon) {
270     if (pt < 10 && fabs(eta) < 1.479) bin = 0;
271 sixie 1.9 if (pt >= 10 && fabs(eta) < 1.479) bin = 1;
272     if (pt < 10 && fabs(eta) >= 1.479) bin = 2;
273 sixie 1.8 if (pt >= 10 && fabs(eta) >= 1.479) bin = 3;
274     } else if (!isGlobal && isTrackerMuon) {
275     bin = 4;
276 sixie 1.9 } else if (isGlobal && !isTrackerMuon) {
277     bin = 5;
278 sixie 1.8 } else {
279     std::cout << "Warning: Muon is not a tracker muon. Such muons are not supported. \n";
280     bin = 0;
281     }
282     }
283 sixie 1.1
284 sixie 1.8 return bin;
285 sixie 1.1 }
286    
287     //--------------------------------------------------------------------------------------------------
288     Double_t MuonIDMVA::MVAValue(Double_t MuPt , Double_t MuEta,
289 sixie 1.2 Double_t MuTkNchi2,
290     Double_t MuGlobalNchi2,
291     Double_t MuNValidHits,
292     Double_t MuNTrackerHits,
293     Double_t MuNPixelHits,
294     Double_t MuNMatches,
295     Double_t MuD0,
296     Double_t MuIP3d,
297     Double_t MuIP3dSig,
298     Double_t MuTrkKink,
299     Double_t MuSegmentCompatibility,
300     Double_t MuCaloCompatibility,
301     Double_t MuHadEnergyOverPt,
302     Double_t MuHoEnergyOverPt,
303     Double_t MuEmEnergyOverPt,
304     Double_t MuHadS9EnergyOverPt,
305     Double_t MuHoS9EnergyOverPt,
306     Double_t MuEmS9EnergyOverPt,
307     Double_t MuChargedIso03OverPt,
308     Double_t MuNeutralIso03OverPt,
309     Double_t MuChargedIso04OverPt,
310 sixie 1.6 Double_t MuNeutralIso04OverPt,
311     Bool_t printDebug
312 sixie 1.1 ) {
313    
314     if (!fIsInitialized) {
315     std::cout << "Error: MuonIDMVA not properly initialized.\n";
316     return -9999;
317     }
318    
319    
320     //set all input variables
321     fMVAVar_MuTkNchi2 = MuTkNchi2;
322     fMVAVar_MuGlobalNchi2 = MuGlobalNchi2;
323     fMVAVar_MuNValidHits = MuNValidHits;
324     fMVAVar_MuNTrackerHits = MuNTrackerHits;
325     fMVAVar_MuNPixelHits = MuNPixelHits;
326     fMVAVar_MuNMatches = MuNMatches;
327     fMVAVar_MuD0 = MuD0;
328     fMVAVar_MuIP3d = MuIP3d;
329     fMVAVar_MuIP3dSig = MuIP3dSig;
330     fMVAVar_MuTrkKink = MuTrkKink;
331     fMVAVar_MuSegmentCompatibility = MuSegmentCompatibility;
332     fMVAVar_MuCaloCompatibility = MuCaloCompatibility;
333     fMVAVar_MuHadEnergyOverPt = MuHadEnergyOverPt;
334     fMVAVar_MuHoEnergyOverPt = MuHoEnergyOverPt;
335     fMVAVar_MuEmEnergyOverPt = MuEmEnergyOverPt;
336     fMVAVar_MuHadS9EnergyOverPt = MuHadS9EnergyOverPt;
337     fMVAVar_MuHoS9EnergyOverPt = MuHoS9EnergyOverPt;
338     fMVAVar_MuEmS9EnergyOverPt = MuEmS9EnergyOverPt;
339 sixie 1.2 fMVAVar_MuChargedIso03OverPt = MuChargedIso03OverPt;
340     fMVAVar_MuNeutralIso03OverPt = MuNeutralIso03OverPt;
341     fMVAVar_MuChargedIso04OverPt = MuChargedIso04OverPt;
342     fMVAVar_MuNeutralIso04OverPt = MuNeutralIso04OverPt;
343 sixie 1.1
344     Double_t mva = -9999;
345     TMVA::Reader *reader = 0;
346 sixie 1.8 reader = fTMVAReader[GetMVABin(MuEta, MuPt, kTRUE, kTRUE )];
347 sixie 1.1
348     mva = reader->EvaluateMVA( fMethodname );
349    
350 sixie 1.6 if (printDebug) {
351 sixie 1.3 std::cout << "Debug Muon MVA: "
352 sixie 1.8 << MuPt << " " << MuEta << " --> MVABin " << GetMVABin(MuEta, MuPt, kTRUE, kTRUE ) << " : "
353 sixie 1.3 << fMVAVar_MuTkNchi2 << " "
354     << fMVAVar_MuGlobalNchi2 << " "
355     << fMVAVar_MuNValidHits << " "
356     << fMVAVar_MuNTrackerHits << " "
357     << fMVAVar_MuNPixelHits << " "
358     << fMVAVar_MuNMatches << " "
359     << fMVAVar_MuD0 << " "
360     << fMVAVar_MuIP3d << " "
361     << fMVAVar_MuIP3dSig << " "
362     << fMVAVar_MuTrkKink << " "
363     << fMVAVar_MuSegmentCompatibility << " "
364     << fMVAVar_MuCaloCompatibility << " "
365     << fMVAVar_MuHadEnergyOverPt << " "
366     << fMVAVar_MuHoEnergyOverPt << " "
367     << fMVAVar_MuEmEnergyOverPt << " "
368     << fMVAVar_MuHadS9EnergyOverPt << " "
369     << fMVAVar_MuHoS9EnergyOverPt << " "
370     << fMVAVar_MuEmS9EnergyOverPt << " "
371     << fMVAVar_MuChargedIso03OverPt << " "
372     << fMVAVar_MuNeutralIso03OverPt << " "
373     << fMVAVar_MuChargedIso04OverPt << " "
374     << fMVAVar_MuNeutralIso04OverPt << " "
375     << " === : === "
376     << mva
377     << std::endl;
378     }
379    
380 sixie 1.1 return mva;
381     }
382    
383    
384 sixie 1.7 //--------------------------------------------------------------------------------------------------
385     Double_t MuonIDMVA::MVAValue(Double_t MuPt , Double_t MuEta,
386     Double_t MuTkNchi2,
387     Double_t MuGlobalNchi2,
388     Double_t MuNValidHits,
389     Double_t MuNTrackerHits,
390     Double_t MuNPixelHits,
391     Double_t MuNMatches,
392     Double_t MuD0,
393     Double_t MuIP3d,
394     Double_t MuIP3dSig,
395     Double_t MuTrkKink,
396     Double_t MuSegmentCompatibility,
397     Double_t MuCaloCompatibility,
398     Double_t MuHadEnergyOverPt,
399     Double_t MuHoEnergyOverPt,
400     Double_t MuEmEnergyOverPt,
401     Double_t MuHadS9EnergyOverPt,
402     Double_t MuHoS9EnergyOverPt,
403     Double_t MuEmS9EnergyOverPt,
404     Double_t MuTrkIso03OverPt,
405     Double_t MuEMIso03OverPt,
406     Double_t MuHadIso03OverPt,
407     Double_t MuTrkIso05OverPt,
408     Double_t MuEMIso05OverPt,
409     Double_t MuHadIso05OverPt,
410     Bool_t printDebug
411     ) {
412    
413     if (!fIsInitialized) {
414     std::cout << "Error: MuonIDMVA not properly initialized.\n";
415     return -9999;
416     }
417    
418     //set all input variables
419     fMVAVar_MuTkNchi2 = MuTkNchi2;
420     fMVAVar_MuGlobalNchi2 = MuGlobalNchi2;
421     fMVAVar_MuNValidHits = MuNValidHits;
422     fMVAVar_MuNTrackerHits = MuNTrackerHits;
423     fMVAVar_MuNPixelHits = MuNPixelHits;
424     fMVAVar_MuNMatches = MuNMatches;
425     fMVAVar_MuD0 = MuD0;
426     fMVAVar_MuIP3d = MuIP3d;
427     fMVAVar_MuIP3dSig = MuIP3dSig;
428     fMVAVar_MuTrkKink = MuTrkKink;
429     fMVAVar_MuSegmentCompatibility = MuSegmentCompatibility;
430     fMVAVar_MuCaloCompatibility = MuCaloCompatibility;
431     fMVAVar_MuHadEnergyOverPt = MuHadEnergyOverPt;
432     fMVAVar_MuHoEnergyOverPt = MuHoEnergyOverPt;
433     fMVAVar_MuEmEnergyOverPt = MuEmEnergyOverPt;
434     fMVAVar_MuHadS9EnergyOverPt = MuHadS9EnergyOverPt;
435     fMVAVar_MuHoS9EnergyOverPt = MuHoS9EnergyOverPt;
436     fMVAVar_MuEmS9EnergyOverPt = MuEmS9EnergyOverPt;
437     fMVAVar_MuTrkIso03OverPt = MuTrkIso03OverPt;
438     fMVAVar_MuEMIso03OverPt = MuEMIso03OverPt;
439     fMVAVar_MuHadIso03OverPt = MuHadIso03OverPt;
440     fMVAVar_MuTrkIso05OverPt = MuTrkIso05OverPt;
441     fMVAVar_MuEMIso05OverPt = MuEMIso05OverPt;
442     fMVAVar_MuHadIso05OverPt = MuHadIso05OverPt;
443    
444     Double_t mva = -9999;
445     TMVA::Reader *reader = 0;
446 sixie 1.8 reader = fTMVAReader[GetMVABin(MuEta, MuPt, kTRUE, kTRUE )];
447 sixie 1.7
448     mva = reader->EvaluateMVA( fMethodname );
449    
450     if (printDebug) {
451     std::cout << "Debug Muon MVA: "
452 sixie 1.8 << MuPt << " " << MuEta << " --> MVABin " << GetMVABin(MuEta, MuPt, kTRUE, kTRUE ) << " : "
453 sixie 1.7 << fMVAVar_MuTkNchi2 << " "
454     << fMVAVar_MuGlobalNchi2 << " "
455     << fMVAVar_MuNValidHits << " "
456     << fMVAVar_MuNTrackerHits << " "
457     << fMVAVar_MuNPixelHits << " "
458     << fMVAVar_MuNMatches << " "
459     << fMVAVar_MuD0 << " "
460     << fMVAVar_MuIP3d << " "
461     << fMVAVar_MuIP3dSig << " "
462     << fMVAVar_MuTrkKink << " "
463     << fMVAVar_MuSegmentCompatibility << " "
464     << fMVAVar_MuCaloCompatibility << " "
465     << fMVAVar_MuHadEnergyOverPt << " "
466     << fMVAVar_MuHoEnergyOverPt << " "
467     << fMVAVar_MuEmEnergyOverPt << " "
468     << fMVAVar_MuHadS9EnergyOverPt << " "
469     << fMVAVar_MuHoS9EnergyOverPt << " "
470     << fMVAVar_MuEmS9EnergyOverPt << " "
471     << fMVAVar_MuTrkIso03OverPt << " "
472     << fMVAVar_MuEMIso03OverPt << " "
473     << fMVAVar_MuHadIso03OverPt << " "
474     << fMVAVar_MuTrkIso05OverPt << " "
475     << fMVAVar_MuEMIso05OverPt << " "
476     << fMVAVar_MuHadIso05OverPt << " "
477     << " === : === "
478     << mva
479     << std::endl;
480     }
481    
482     return mva;
483     }
484    
485    
486 sixie 1.9 Double_t MuonIDMVA::MVAValue_IsoRings( Double_t MuPt,
487     Double_t MuEta,
488     Double_t ChargedIso_DR0p0To0p1,
489     Double_t ChargedIso_DR0p1To0p2,
490     Double_t ChargedIso_DR0p2To0p3,
491     Double_t ChargedIso_DR0p3To0p4,
492     Double_t ChargedIso_DR0p4To0p5,
493     Double_t GammaIso_DR0p0To0p1,
494     Double_t GammaIso_DR0p1To0p2,
495     Double_t GammaIso_DR0p2To0p3,
496     Double_t GammaIso_DR0p3To0p4,
497     Double_t GammaIso_DR0p4To0p5,
498     Double_t NeutralHadronIso_DR0p0To0p1,
499     Double_t NeutralHadronIso_DR0p1To0p2,
500     Double_t NeutralHadronIso_DR0p2To0p3,
501     Double_t NeutralHadronIso_DR0p3To0p4,
502     Double_t NeutralHadronIso_DR0p4To0p5,
503     Bool_t printDebug) {
504    
505     if (fMVAType != MuonIDMVA::kIsoRingsV0) {
506     std::cout << "Error: This function is only supported for MVAType == kIsoRingsV0.\n" << std::endl;
507     assert(kFALSE);
508     }
509    
510     fMVAVar_MuPt = MuPt;
511     fMVAVar_MuEta = MuEta;
512     fMVAVar_ChargedIso_DR0p0To0p1 = ChargedIso_DR0p0To0p1;
513     fMVAVar_ChargedIso_DR0p1To0p2 = ChargedIso_DR0p1To0p2;
514     fMVAVar_ChargedIso_DR0p2To0p3 = ChargedIso_DR0p2To0p3;
515     fMVAVar_ChargedIso_DR0p3To0p4 = ChargedIso_DR0p3To0p4;
516     fMVAVar_ChargedIso_DR0p4To0p5 = ChargedIso_DR0p4To0p5;
517     fMVAVar_GammaIso_DR0p0To0p1 = GammaIso_DR0p0To0p1;
518     fMVAVar_GammaIso_DR0p1To0p2 = GammaIso_DR0p1To0p2;
519     fMVAVar_GammaIso_DR0p2To0p3 = GammaIso_DR0p2To0p3;
520     fMVAVar_GammaIso_DR0p3To0p4 = GammaIso_DR0p3To0p4;
521     fMVAVar_GammaIso_DR0p4To0p5 = GammaIso_DR0p4To0p5;
522     fMVAVar_NeutralHadronIso_DR0p0To0p1 = NeutralHadronIso_DR0p0To0p1;
523     fMVAVar_NeutralHadronIso_DR0p1To0p2 = NeutralHadronIso_DR0p1To0p2;
524     fMVAVar_NeutralHadronIso_DR0p2To0p3 = NeutralHadronIso_DR0p2To0p3;
525     fMVAVar_NeutralHadronIso_DR0p3To0p4 = NeutralHadronIso_DR0p3To0p4;
526     fMVAVar_NeutralHadronIso_DR0p4To0p5 = NeutralHadronIso_DR0p4To0p5;
527    
528     Double_t mva = -9999;
529     TMVA::Reader *reader = 0;
530    
531     if (printDebug == kTRUE) {
532     std::cout <<" -> BIN: " << fMVAVar_MuEta << " " << fMVAVar_MuPt << " : "
533     << GetMVABin( fMVAVar_MuEta , fMVAVar_MuPt) << std::endl;
534     }
535     reader = fTMVAReader[GetMVABin( fMVAVar_MuEta , fMVAVar_MuPt)];
536     mva = reader->EvaluateMVA( fMethodname );
537    
538     if (printDebug == kTRUE) {
539    
540     std::cout << "Debug Muon MVA: \n";
541     std::cout << fMVAVar_ChargedIso_DR0p0To0p1 << " "
542     << fMVAVar_ChargedIso_DR0p1To0p2 << " "
543     << fMVAVar_ChargedIso_DR0p2To0p3 << " "
544     << fMVAVar_ChargedIso_DR0p3To0p4 << " "
545     << fMVAVar_ChargedIso_DR0p4To0p5 << " "
546     << fMVAVar_GammaIso_DR0p0To0p1 << " "
547     << fMVAVar_GammaIso_DR0p1To0p2 << " "
548     << fMVAVar_GammaIso_DR0p2To0p3 << " "
549     << fMVAVar_GammaIso_DR0p3To0p4 << " "
550     << fMVAVar_GammaIso_DR0p4To0p5 << " "
551     << fMVAVar_NeutralHadronIso_DR0p0To0p1 << " "
552     << fMVAVar_NeutralHadronIso_DR0p1To0p2 << " "
553     << fMVAVar_NeutralHadronIso_DR0p2To0p3 << " "
554     << fMVAVar_NeutralHadronIso_DR0p3To0p4 << " "
555     << fMVAVar_NeutralHadronIso_DR0p4To0p5 << " "
556     << std::endl;
557     std::cout << "MVA: " << mva << " "
558     << std::endl;
559     }
560     return mva;
561     }
562    
563    
564    
565     Double_t MuonIDMVA::MVAValue_ID( Double_t MuPt,
566     Double_t MuEta,
567     Bool_t MuIsGlobal,
568     Bool_t MuIsTracker,
569     Double_t MuTkNchi2,
570     Double_t MuGlobalNchi2,
571     Double_t MuNValidHits,
572     Double_t MuNTrackerHits,
573     Double_t MuNPixelHits,
574     Double_t MuNMatches,
575     Double_t MuTrkKink,
576     Double_t MuSegmentCompatibility,
577     Double_t MuCaloCompatibility,
578     Double_t MuHadEnergy,
579     Double_t MuEmEnergy,
580     Double_t MuHadS9Energy,
581     Double_t MuEmS9Energy,
582     Bool_t printDebug) {
583    
584     if (fMVAType != MuonIDMVA::kIDV0) {
585     std::cout << "Error: This function is only supported for MVAType == kIDV0.\n" << std::endl;
586     assert(kFALSE);
587     }
588    
589     fMVAVar_MuPt = MuPt;
590     fMVAVar_MuEta = MuEta;
591    
592     fMVAVar_MuTkNchi2 = MuTkNchi2;
593     fMVAVar_MuGlobalNchi2 = MuGlobalNchi2;
594     fMVAVar_MuNValidHits = MuNValidHits;
595     fMVAVar_MuNTrackerHits = MuNTrackerHits;
596     fMVAVar_MuNPixelHits = MuNPixelHits;
597     fMVAVar_MuNMatches = MuNMatches;
598     fMVAVar_MuTrkKink = MuTrkKink;
599     fMVAVar_MuSegmentCompatibility = MuSegmentCompatibility;
600     fMVAVar_MuCaloCompatibility = MuCaloCompatibility;
601     fMVAVar_MuHadEnergy = MuHadEnergy;
602     fMVAVar_MuEmEnergy = MuEmEnergy;
603     fMVAVar_MuHadS9Energy = MuHadS9Energy;
604     fMVAVar_MuEmS9Energy = MuEmS9Energy;
605    
606     Double_t mva = -9999;
607     TMVA::Reader *reader = 0;
608    
609     if (printDebug == kTRUE) {
610     std::cout <<" -> BIN: " << fMVAVar_MuEta << " " << fMVAVar_MuPt << " : "
611     << GetMVABin( fMVAVar_MuEta , fMVAVar_MuPt, MuIsGlobal, MuIsTracker) << std::endl;
612     }
613     reader = fTMVAReader[GetMVABin( fMVAVar_MuEta , fMVAVar_MuPt, MuIsGlobal, MuIsTracker)];
614     mva = reader->EvaluateMVA( fMethodname );
615    
616     if (printDebug == kTRUE) {
617    
618     std::cout << "Debug Muon MVA: \n";
619     std::cout << fMVAVar_MuTkNchi2 << " "
620     << fMVAVar_MuGlobalNchi2 << " "
621     << fMVAVar_MuNValidHits << " "
622     << fMVAVar_MuNTrackerHits << " "
623     << fMVAVar_MuNPixelHits << " "
624     << fMVAVar_MuNMatches << " "
625     << fMVAVar_MuTrkKink << " "
626     << fMVAVar_MuSegmentCompatibility << " "
627     << fMVAVar_MuCaloCompatibility << " "
628     << fMVAVar_MuHadEnergy << " "
629     << fMVAVar_MuEmEnergy << " "
630     << fMVAVar_MuHadS9Energy << " "
631     << fMVAVar_MuEmS9Energy << " "
632     << std::endl;
633     std::cout << "MVA: " << mva << " "
634     << std::endl;
635     }
636     return mva;
637     }
638    
639    
640 sixie 1.1
641 sixie 1.4 //--------------------------------------------------------------------------------------------------
642     Double_t MuonIDMVA::MVAValue(const Muon *mu, const Vertex *vertex, MuonTools *fMuonTools,
643     const PFCandidateCol *PFCands,
644 sixie 1.6 const PileupEnergyDensityCol *PileupEnergyDensity,
645     Bool_t printDebug) {
646 sixie 1.1
647 sixie 1.4 if (!fIsInitialized) {
648     std::cout << "Error: MuonIDMVA not properly initialized.\n";
649     return -9999;
650     }
651    
652     const Track *muTrk=0;
653     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
654     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
655 sixie 1.1
656 sixie 1.4 Double_t muNchi2 = 0.0;
657     if(mu->HasGlobalTrk()) { muNchi2 = mu->GlobalTrk()->RChi2(); }
658     else if(mu->HasStandaloneTrk()) { muNchi2 = mu->StandaloneTrk()->RChi2(); }
659     else if(mu->HasTrackerTrk()) { muNchi2 = mu->TrackerTrk()->RChi2(); }
660    
661 sixie 1.7 Double_t ChargedIso03 = 0;
662     Double_t NeutralIso03_05Threshold = 0;
663     Double_t ChargedIso04 = 0;
664     Double_t NeutralIso04_05Threshold = 0;
665     ChargedIso03 = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.1, 99999, 0.3, 0.0, 0.0);
666     NeutralIso03_05Threshold = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.0, 0.5, 0.3, 0.0, 0.0);
667     ChargedIso04 = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.1, 99999, 0.4, 0.0, 0.0);
668     NeutralIso04_05Threshold = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.0, 0.5, 0.4, 0.0, 0.0);
669    
670 sixie 1.4 Double_t Rho = 0;
671     if (!(TMath::IsNaN(PileupEnergyDensity->At(0)->Rho()) || isinf(PileupEnergyDensity->At(0)->Rho()))) Rho = PileupEnergyDensity->At(0)->Rho();
672    
673 sixie 1.8
674 sixie 1.4 //set all input variables
675     fMVAVar_MuTkNchi2 = muTrk->RChi2();
676     fMVAVar_MuGlobalNchi2 = muNchi2;
677     fMVAVar_MuNValidHits = mu->NValidHits();
678     fMVAVar_MuNTrackerHits = muTrk->NHits();
679     fMVAVar_MuNPixelHits = muTrk->NPixelHits();
680     fMVAVar_MuNMatches = mu->NMatches();
681     fMVAVar_MuD0 = muTrk->D0Corrected(*vertex);
682     fMVAVar_MuIP3d = mu->Ip3dPV();
683     fMVAVar_MuIP3dSig = mu->Ip3dPVSignificance();
684     fMVAVar_MuTrkKink = mu->TrkKink();
685     fMVAVar_MuSegmentCompatibility = fMuonTools->GetSegmentCompatability(mu);
686     fMVAVar_MuCaloCompatibility = fMuonTools->GetCaloCompatability(mu, kTRUE, kTRUE);
687     fMVAVar_MuHadEnergyOverPt = (mu->HadEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadEnergy,muTrk->Eta()))/muTrk->Pt();
688     fMVAVar_MuHoEnergyOverPt = (mu->HoEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHoEnergy,muTrk->Eta()))/muTrk->Pt();
689     fMVAVar_MuEmEnergyOverPt = (mu->EmEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEmEnergy,muTrk->Eta()))/muTrk->Pt();
690     fMVAVar_MuHadS9EnergyOverPt = (mu->HadS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadS9Energy,muTrk->Eta()))/muTrk->Pt();
691     fMVAVar_MuHoS9EnergyOverPt = (mu->HoS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHoS9Energy,muTrk->Eta()))/muTrk->Pt();
692     fMVAVar_MuEmS9EnergyOverPt = (mu->EmS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEmS9Energy,muTrk->Eta()))/muTrk->Pt();
693     fMVAVar_MuChargedIso03OverPt = (ChargedIso03 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso03,muTrk->Eta()))/muTrk->Pt();
694     fMVAVar_MuNeutralIso03OverPt = (NeutralIso03_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03,muTrk->Eta()))/muTrk->Pt();
695     fMVAVar_MuChargedIso04OverPt = (ChargedIso04 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso04,muTrk->Eta()))/muTrk->Pt();
696     fMVAVar_MuNeutralIso04OverPt = (NeutralIso04_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso04,muTrk->Eta()))/muTrk->Pt();
697 sixie 1.7 fMVAVar_MuTrkIso03OverPt = (mu->IsoR03SumPt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuTrkIso03,muTrk->Eta()))/muTrk->Pt();
698     fMVAVar_MuEMIso03OverPt = (mu->IsoR03EmEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03,muTrk->Eta()))/muTrk->Pt();
699     fMVAVar_MuHadIso03OverPt = (mu->IsoR03HadEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03,muTrk->Eta()))/muTrk->Pt();
700     fMVAVar_MuTrkIso05OverPt = (mu->IsoR05SumPt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuTrkIso05,muTrk->Eta()))/muTrk->Pt();
701     fMVAVar_MuEMIso05OverPt = (mu->IsoR05EmEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso05,muTrk->Eta()))/muTrk->Pt();
702     fMVAVar_MuHadIso05OverPt = (mu->IsoR05HadEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso05,muTrk->Eta()))/muTrk->Pt();
703 sixie 1.4
704     Double_t mva = -9999;
705     TMVA::Reader *reader = 0;
706 sixie 1.8 reader = fTMVAReader[GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon())];
707 sixie 1.1
708 sixie 1.4 mva = reader->EvaluateMVA( fMethodname );
709    
710 sixie 1.6 if (printDebug) {
711 sixie 1.4 std::cout << "Debug Muon MVA: "
712     << mu->Pt() << " " << mu->Eta() << " " << mu->Phi() << " : "
713 sixie 1.8 << muTrk->Pt() << " " << muTrk->Eta() << " --> MVABin " << GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon()) << " : "
714 sixie 1.4 << fMVAVar_MuTkNchi2 << " "
715     << fMVAVar_MuGlobalNchi2 << " "
716     << fMVAVar_MuNValidHits << " "
717     << fMVAVar_MuNTrackerHits << " "
718     << fMVAVar_MuNPixelHits << " "
719     << fMVAVar_MuNMatches << " "
720     << fMVAVar_MuD0 << " "
721     << fMVAVar_MuIP3d << " "
722     << fMVAVar_MuIP3dSig << " "
723     << fMVAVar_MuTrkKink << " "
724     << fMVAVar_MuSegmentCompatibility << " "
725     << fMVAVar_MuCaloCompatibility << " "
726     << fMVAVar_MuHadEnergyOverPt << " "
727     << fMVAVar_MuHoEnergyOverPt << " "
728     << fMVAVar_MuEmEnergyOverPt << " "
729     << fMVAVar_MuHadS9EnergyOverPt << " "
730     << fMVAVar_MuHoS9EnergyOverPt << " "
731     << fMVAVar_MuEmS9EnergyOverPt << " "
732     << fMVAVar_MuChargedIso03OverPt << " "
733     << fMVAVar_MuNeutralIso03OverPt << " "
734     << fMVAVar_MuChargedIso04OverPt << " "
735     << fMVAVar_MuNeutralIso04OverPt << " "
736 sixie 1.7 << fMVAVar_MuTrkIso03OverPt << " "
737     << fMVAVar_MuEMIso03OverPt << " "
738     << fMVAVar_MuHadIso03OverPt << " "
739     << fMVAVar_MuTrkIso05OverPt << " "
740     << fMVAVar_MuEMIso05OverPt << " "
741     << fMVAVar_MuHadIso05OverPt << " "
742 sixie 1.4 << " === : === "
743     << mva
744     << std::endl;
745     }
746 sixie 1.1
747 sixie 1.4 return mva;
748     }
749 sixie 1.8
750    
751     //--------------------------------------------------------------------------------------------------
752     Double_t MuonIDMVA::MVAValue(const Muon *mu, const Vertex *vertex, MuonTools *fMuonTools,
753     const PFCandidateCol *PFCands,
754     const PileupEnergyDensityCol *PileupEnergyDensity,
755     MuonTools::EMuonEffectiveAreaTarget EffectiveAreaTarget,
756     const ElectronCol *goodElectrons,
757     const MuonCol *goodMuons,
758     Bool_t printDebug) {
759    
760     if (!fIsInitialized) {
761     std::cout << "Error: MuonIDMVA not properly initialized.\n";
762     return -9999;
763     }
764    
765     const Track *muTrk=0;
766     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
767     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
768    
769     Double_t muNchi2 = 0.0;
770     if(mu->HasGlobalTrk()) { muNchi2 = mu->GlobalTrk()->RChi2(); }
771     else if(mu->HasStandaloneTrk()) { muNchi2 = mu->StandaloneTrk()->RChi2(); }
772     else if(mu->HasTrackerTrk()) { muNchi2 = mu->TrackerTrk()->RChi2(); }
773    
774     Double_t ChargedIso03 = 0;
775     Double_t NeutralIso03_05Threshold = 0;
776     Double_t ChargedIso04 = 0;
777     Double_t NeutralIso04_05Threshold = 0;
778     ChargedIso03 = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.1, 99999, 0.3, 0.0, 0.0);
779     NeutralIso03_05Threshold = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.0, 0.5, 0.3, 0.0, 0.0);
780     ChargedIso04 = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.1, 99999, 0.4, 0.0, 0.0);
781     NeutralIso04_05Threshold = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.0, 0.5, 0.4, 0.0, 0.0);
782    
783     Double_t Rho = 0;
784     if (!(TMath::IsNaN(PileupEnergyDensity->At(0)->Rho()) || isinf(PileupEnergyDensity->At(0)->Rho()))) Rho = PileupEnergyDensity->At(0)->Rho();
785    
786    
787     //set all input variables
788     fMVAVar_MuPt = muTrk->Pt();
789     fMVAVar_MuEta = muTrk->Eta();
790     fMVAVar_MuTkNchi2 = muTrk->RChi2();
791     fMVAVar_MuGlobalNchi2 = muNchi2;
792     fMVAVar_MuNValidHits = mu->NValidHits();
793     fMVAVar_MuNTrackerHits = muTrk->NHits();
794     fMVAVar_MuNPixelHits = muTrk->NPixelHits();
795     fMVAVar_MuNMatches = mu->NMatches();
796     fMVAVar_MuD0 = muTrk->D0Corrected(*vertex);
797     fMVAVar_MuIP3d = mu->Ip3dPV();
798     fMVAVar_MuIP3dSig = mu->Ip3dPVSignificance();
799     fMVAVar_MuTrkKink = mu->TrkKink();
800     fMVAVar_MuSegmentCompatibility = fMuonTools->GetSegmentCompatability(mu);
801     fMVAVar_MuCaloCompatibility = fMuonTools->GetCaloCompatability(mu, kTRUE, kTRUE);
802     fMVAVar_MuHadEnergy = mu->HadEnergy() ;
803     fMVAVar_MuEmEnergy = mu->EmEnergy();
804     fMVAVar_MuHadS9Energy = mu->HadS9Energy();
805     fMVAVar_MuEmS9Energy = mu->EmS9Energy();
806     fMVAVar_MuChargedIso03OverPt = (ChargedIso03 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso03,muTrk->Eta()))/muTrk->Pt();
807     fMVAVar_MuNeutralIso03OverPt = (NeutralIso03_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03,muTrk->Eta()))/muTrk->Pt();
808     fMVAVar_MuChargedIso04OverPt = (ChargedIso04 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso04,muTrk->Eta()))/muTrk->Pt();
809     fMVAVar_MuNeutralIso04OverPt = (NeutralIso04_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso04,muTrk->Eta()))/muTrk->Pt();
810     fMVAVar_MuTrkIso03OverPt = (mu->IsoR03SumPt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuTrkIso03,muTrk->Eta()))/muTrk->Pt();
811     fMVAVar_MuEMIso03OverPt = (mu->IsoR03EmEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03,muTrk->Eta()))/muTrk->Pt();
812     fMVAVar_MuHadIso03OverPt = (mu->IsoR03HadEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03,muTrk->Eta()))/muTrk->Pt();
813     fMVAVar_MuTrkIso05OverPt = (mu->IsoR05SumPt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuTrkIso05,muTrk->Eta()))/muTrk->Pt();
814     fMVAVar_MuEMIso05OverPt = (mu->IsoR05EmEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso05,muTrk->Eta()))/muTrk->Pt();
815     fMVAVar_MuHadIso05OverPt = (mu->IsoR05HadEt() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso05,muTrk->Eta()))/muTrk->Pt();
816    
817    
818     Double_t tmpChargedIso_DR0p0To0p1 = 0;
819     Double_t tmpChargedIso_DR0p1To0p2 = 0;
820     Double_t tmpChargedIso_DR0p2To0p3 = 0;
821     Double_t tmpChargedIso_DR0p3To0p4 = 0;
822     Double_t tmpChargedIso_DR0p4To0p5 = 0;
823     Double_t tmpGammaIso_DR0p0To0p1 = 0;
824     Double_t tmpGammaIso_DR0p1To0p2 = 0;
825     Double_t tmpGammaIso_DR0p2To0p3 = 0;
826     Double_t tmpGammaIso_DR0p3To0p4 = 0;
827     Double_t tmpGammaIso_DR0p4To0p5 = 0;
828     Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
829     Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
830     Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
831     Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
832     Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
833    
834     for (UInt_t p=0; p<PFCands->GetEntries();p++) {
835     const PFCandidate *pf = PFCands->At(p);
836    
837     //exclude the muon itself
838     if(pf->TrackerTrk() && mu->TrackerTrk() &&
839     pf->TrackerTrk() == mu->TrackerTrk()) continue;
840    
841     //************************************************************
842     // New Isolation Calculations
843     //************************************************************
844     Double_t dr = MathUtils::DeltaR(mu->Mom(), pf->Mom());
845    
846     if (dr < 0.5) {
847     Bool_t IsLeptonFootprint = kFALSE;
848     //************************************************************
849     // Lepton Footprint Removal
850     //************************************************************
851     for (UInt_t q=0; q < goodElectrons->GetEntries() ; ++q) {
852     //if pf candidate matches an electron passing ID cuts, then veto it
853     if(pf->GsfTrk() && goodElectrons->At(q)->GsfTrk() &&
854     pf->GsfTrk() == goodElectrons->At(q)->GsfTrk()) IsLeptonFootprint = kTRUE;
855     if(pf->TrackerTrk() && goodElectrons->At(q)->TrackerTrk() &&
856     pf->TrackerTrk() == goodElectrons->At(q)->TrackerTrk()) IsLeptonFootprint = kTRUE;
857     //if pf candidate lies in veto regions of electron passing ID cuts, then veto it
858     if(pf->BestTrk() && fabs(goodElectrons->At(q)->SCluster()->Eta()) >= 1.479
859     && MathUtils::DeltaR(goodElectrons->At(q)->Mom(), pf->Mom()) < 0.015) IsLeptonFootprint = kTRUE;
860     if(pf->PFType() == PFCandidate::eGamma && fabs(goodElectrons->At(q)->SCluster()->Eta()) >= 1.479 &&
861     MathUtils::DeltaR(goodElectrons->At(q)->Mom(), pf->Mom()) < 0.08) IsLeptonFootprint = kTRUE;
862     }
863     for (UInt_t q=0; q < goodMuons->GetEntries() ; ++q) {
864     //if pf candidate matches an muon passing ID cuts, then veto it
865     if(pf->TrackerTrk() && goodMuons->At(q)->TrackerTrk() &&
866     pf->TrackerTrk() == goodMuons->At(q)->TrackerTrk()) IsLeptonFootprint = kTRUE;
867     //if pf candidate lies in veto regions of muon passing ID cuts, then veto it
868     if(pf->BestTrk() && MathUtils::DeltaR(goodMuons->At(q)->Mom(), pf->Mom()) < 0.01) IsLeptonFootprint = kTRUE;
869     }
870    
871     if (!IsLeptonFootprint) {
872     Bool_t passVeto = kTRUE;
873     //Charged
874     if(pf->BestTrk()) {
875     if (!(fabs(pf->BestTrk()->DzCorrected(*vertex) - mu->BestTrk()->DzCorrected(*vertex)) < 0.2)) passVeto = kFALSE;
876     //************************************************************
877     // Veto any PFmuon, or PFEle
878     if (pf->PFType() == PFCandidate::eElectron || pf->PFType() == PFCandidate::eMuon) passVeto = kFALSE;
879     //************************************************************
880     //************************************************************
881     // Footprint Veto
882     if (dr < 0.01) passVeto = kFALSE;
883     //************************************************************
884     if (passVeto) {
885     if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
886     if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
887     if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
888     if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
889     if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
890     } //pass veto
891     }
892     //Gamma
893     else if (pf->PFType() == PFCandidate::eGamma) {
894     if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
895     if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
896     if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
897     if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
898     if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
899     }
900     //NeutralHadron
901     else {
902     if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
903     if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
904     if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
905     if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
906     if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
907     }
908     } //not lepton footprint
909     } //in 1.0 dr cone
910     } //loop over PF candidates
911    
912     Double_t fMVAVar_ChargedIso_DR0p0To0p1 = 0;
913     Double_t fMVAVar_ChargedIso_DR0p1To0p2 = 0;
914     Double_t fMVAVar_ChargedIso_DR0p2To0p3 = 0;
915     Double_t fMVAVar_ChargedIso_DR0p3To0p4 = 0;
916     Double_t fMVAVar_ChargedIso_DR0p4To0p5 = 0;
917     Double_t fMVAVar_GammaIso_DR0p0To0p1 = 0;
918     Double_t fMVAVar_GammaIso_DR0p1To0p2 = 0;
919     Double_t fMVAVar_GammaIso_DR0p2To0p3 = 0;
920     Double_t fMVAVar_GammaIso_DR0p3To0p4 = 0;
921     Double_t fMVAVar_GammaIso_DR0p4To0p5 = 0;
922     Double_t fMVAVar_NeutralHadronIso_DR0p0To0p1 = 0;
923     Double_t fMVAVar_NeutralHadronIso_DR0p1To0p2 = 0;
924     Double_t fMVAVar_NeutralHadronIso_DR0p2To0p3 = 0;
925     Double_t fMVAVar_NeutralHadronIso_DR0p3To0p4 = 0;
926     Double_t fMVAVar_NeutralHadronIso_DR0p4To0p5 = 0;
927    
928     fMVAVar_ChargedIso_DR0p0To0p1 = TMath::Min((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
929     fMVAVar_ChargedIso_DR0p1To0p2 = TMath::Min((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
930     fMVAVar_ChargedIso_DR0p2To0p3 = TMath::Min((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
931     fMVAVar_ChargedIso_DR0p3To0p4 = TMath::Min((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
932     fMVAVar_ChargedIso_DR0p4To0p5 = TMath::Min((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);
933     fMVAVar_GammaIso_DR0p0To0p1 = TMath::Max(TMath::Min((tmpGammaIso_DR0p0To0p1 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuGammaIsoDR0p0To0p1, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
934     fMVAVar_GammaIso_DR0p1To0p2 = TMath::Max(TMath::Min((tmpGammaIso_DR0p1To0p2 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuGammaIsoDR0p1To0p2, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
935     fMVAVar_GammaIso_DR0p2To0p3 = TMath::Max(TMath::Min((tmpGammaIso_DR0p2To0p3 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuGammaIsoDR0p2To0p3, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
936     fMVAVar_GammaIso_DR0p3To0p4 = TMath::Max(TMath::Min((tmpGammaIso_DR0p3To0p4 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuGammaIsoDR0p3To0p4, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
937     fMVAVar_GammaIso_DR0p4To0p5 = TMath::Max(TMath::Min((tmpGammaIso_DR0p4To0p5 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuGammaIsoDR0p4To0p5, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
938     fMVAVar_NeutralHadronIso_DR0p0To0p1 = TMath::Max(TMath::Min((tmpNeutralHadronIso_DR0p0To0p1 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralHadronIsoDR0p0To0p1, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
939     fMVAVar_NeutralHadronIso_DR0p1To0p2 = TMath::Max(TMath::Min((tmpNeutralHadronIso_DR0p1To0p2 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralHadronIsoDR0p1To0p2, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
940     fMVAVar_NeutralHadronIso_DR0p2To0p3 = TMath::Max(TMath::Min((tmpNeutralHadronIso_DR0p2To0p3 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralHadronIsoDR0p2To0p3, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
941     fMVAVar_NeutralHadronIso_DR0p3To0p4 = TMath::Max(TMath::Min((tmpNeutralHadronIso_DR0p3To0p4 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralHadronIsoDR0p3To0p4, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
942     fMVAVar_NeutralHadronIso_DR0p4To0p5 = TMath::Max(TMath::Min((tmpNeutralHadronIso_DR0p4To0p5 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralHadronIsoDR0p4To0p5, mu->Eta(), EffectiveAreaTarget))/mu->Pt(), 2.5), 0.0);
943    
944    
945    
946     Double_t mva = -9999;
947     TMVA::Reader *reader = 0;
948    
949     if (printDebug) {
950     std::cout <<" -> BIN: " << fMVAVar_MuEta << " " << fMVAVar_MuPt << " : "
951     << GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon() )
952     << std::endl;
953     }
954    
955     reader = fTMVAReader[GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon() )];
956    
957     mva = reader->EvaluateMVA( fMethodname );
958    
959     if (printDebug) {
960     std::cout << "Debug Muon MVA: \n";
961     std::cout << " MuTkNchi2 " << fMVAVar_MuTkNchi2
962     << " MuGlobalNchi2 " << fMVAVar_MuGlobalNchi2
963     << " MuNValidHits " << fMVAVar_MuNValidHits
964     << " MuNTrackerHits " << fMVAVar_MuNTrackerHits
965     << " MuNPixelHits " << fMVAVar_MuNPixelHits
966     << " MuNMatches " << fMVAVar_MuNMatches
967     << " MuD0 " << fMVAVar_MuD0
968     << " MuIP3d " << fMVAVar_MuIP3d
969     << " MuIP3dSig " << fMVAVar_MuIP3dSig
970     << " MuTrkKink " << fMVAVar_MuTrkKink
971     << " MuSegmentCompatibility " << fMVAVar_MuSegmentCompatibility
972     << " MuCaloCompatibility " << fMVAVar_MuCaloCompatibility
973     << " MuHadEnergy " << fMVAVar_MuHadEnergy
974     << " MuEmEnergy " << fMVAVar_MuEmEnergy
975     << " MuHadS9Energy " << fMVAVar_MuHadS9Energy
976     << " MuEmS9Energy " << fMVAVar_MuEmS9Energy
977     << " eta " << fMVAVar_MuEta
978     << " pt " << fMVAVar_MuPt << std::endl;
979    
980     std::cout << fMVAVar_ChargedIso_DR0p0To0p1 << " "
981     << fMVAVar_ChargedIso_DR0p1To0p2 << " "
982     << fMVAVar_ChargedIso_DR0p2To0p3 << " "
983     << fMVAVar_ChargedIso_DR0p3To0p4 << " "
984     << fMVAVar_ChargedIso_DR0p4To0p5 << " "
985     << fMVAVar_GammaIso_DR0p0To0p1 << " "
986     << fMVAVar_GammaIso_DR0p1To0p2 << " "
987     << fMVAVar_GammaIso_DR0p2To0p3 << " "
988     << fMVAVar_GammaIso_DR0p3To0p4 << " "
989     << fMVAVar_GammaIso_DR0p4To0p5 << " "
990     << fMVAVar_NeutralHadronIso_DR0p0To0p1 << " "
991     << fMVAVar_NeutralHadronIso_DR0p1To0p2 << " "
992     << fMVAVar_NeutralHadronIso_DR0p2To0p3 << " "
993     << fMVAVar_NeutralHadronIso_DR0p3To0p4 << " "
994     << fMVAVar_NeutralHadronIso_DR0p4To0p5 << " "
995     << std::endl;
996     std::cout << "MVA: " << mva
997     << std::endl;
998     }
999    
1000     return mva;
1001     }