ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MuonIDMVA.cc
Revision: 1.5
Committed: Wed Jan 4 14:39:16 2012 UTC (13 years, 4 months ago) by sixie
Content type: text/plain
Branch: MAIN
Changes since 1.4: +2 -2 lines
Log Message:
turn off debug printout

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     fIsInitialized(kFALSE)
19     {
20     // Constructor.
21     for(UInt_t i=0; i<6; ++i) {
22     fTMVAReader[i] = 0;
23     }
24     }
25    
26    
27     //--------------------------------------------------------------------------------------------------
28     MuonIDMVA::~MuonIDMVA()
29     {
30     for(UInt_t i=0; i<6; ++i) {
31     if (fTMVAReader[i]) delete fTMVAReader[i];
32     }
33     }
34    
35     //--------------------------------------------------------------------------------------------------
36     void MuonIDMVA::Initialize( TString methodName,
37     TString Subdet0Pt10To14p5Weights ,
38     TString Subdet1Pt10To14p5Weights ,
39     TString Subdet0Pt14p5To20Weights,
40     TString Subdet1Pt14p5To20Weights,
41     TString Subdet0Pt20ToInfWeights,
42     TString Subdet1Pt20ToInfWeights,
43     MuonIDMVA::MVAType type) {
44    
45     fIsInitialized = kTRUE;
46    
47     fMethodname = methodName;
48    
49     for(UInt_t i=0; i<6; ++i) {
50     if (fTMVAReader[i]) delete fTMVAReader[i];
51    
52     fTMVAReader[i] = new TMVA::Reader( "!Color:!Silent:Error" );
53     fTMVAReader[i]->SetVerbose(kTRUE);
54    
55 sixie 1.2 if (type == kV2) {
56     fTMVAReader[i]->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
57     fTMVAReader[i]->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
58     fTMVAReader[i]->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
59     fTMVAReader[i]->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
60     fTMVAReader[i]->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
61     fTMVAReader[i]->AddVariable( "NMatches", &fMVAVar_MuNMatches );
62     fTMVAReader[i]->AddVariable( "D0", &fMVAVar_MuD0 );
63     fTMVAReader[i]->AddVariable( "IP3d", &fMVAVar_MuIP3d );
64     fTMVAReader[i]->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
65     fTMVAReader[i]->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
66     fTMVAReader[i]->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
67     }
68    
69 sixie 1.1 if (type == kV3) {
70     fTMVAReader[i]->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
71     fTMVAReader[i]->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
72     fTMVAReader[i]->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
73     fTMVAReader[i]->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
74     fTMVAReader[i]->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
75     fTMVAReader[i]->AddVariable( "NMatches", &fMVAVar_MuNMatches );
76     fTMVAReader[i]->AddVariable( "D0", &fMVAVar_MuD0 );
77     fTMVAReader[i]->AddVariable( "IP3d", &fMVAVar_MuIP3d );
78     fTMVAReader[i]->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
79     fTMVAReader[i]->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
80     fTMVAReader[i]->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
81     fTMVAReader[i]->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
82     fTMVAReader[i]->AddVariable( "HadEnergyOverPt", &fMVAVar_MuHadEnergyOverPt );
83     if (i==0 || i==2 || i==4) {
84     fTMVAReader[i]->AddVariable( "HoEnergyOverPt", &fMVAVar_MuHoEnergyOverPt );
85     }
86     fTMVAReader[i]->AddVariable( "EmEnergyOverPt", &fMVAVar_MuEmEnergyOverPt );
87     fTMVAReader[i]->AddVariable( "HadS9EnergyOverPt", &fMVAVar_MuHadS9EnergyOverPt );
88     if (i==0 || i==2 || i==4) {
89     fTMVAReader[i]->AddVariable( "HoS9EnergyOverPt", &fMVAVar_MuHoS9EnergyOverPt );
90     }
91     fTMVAReader[i]->AddVariable( "EmS9EnergyOverPt", &fMVAVar_MuEmS9EnergyOverPt );
92     }
93 sixie 1.2
94     if (type == kV8) {
95     fTMVAReader[i]->AddVariable( "TkNchi2", &fMVAVar_MuTkNchi2 );
96     fTMVAReader[i]->AddVariable( "GlobalNchi2", &fMVAVar_MuGlobalNchi2 );
97     fTMVAReader[i]->AddVariable( "NValidHits", &fMVAVar_MuNValidHits );
98     fTMVAReader[i]->AddVariable( "NTrackerHits", &fMVAVar_MuNTrackerHits );
99     fTMVAReader[i]->AddVariable( "NPixelHits", &fMVAVar_MuNPixelHits );
100     fTMVAReader[i]->AddVariable( "NMatches", &fMVAVar_MuNMatches );
101     fTMVAReader[i]->AddVariable( "D0", &fMVAVar_MuD0 );
102     fTMVAReader[i]->AddVariable( "IP3d", &fMVAVar_MuIP3d );
103     fTMVAReader[i]->AddVariable( "IP3dSig", &fMVAVar_MuIP3dSig );
104     fTMVAReader[i]->AddVariable( "TrkKink", &fMVAVar_MuTrkKink );
105     fTMVAReader[i]->AddVariable( "SegmentCompatibility", &fMVAVar_MuSegmentCompatibility );
106     fTMVAReader[i]->AddVariable( "CaloCompatibility", &fMVAVar_MuCaloCompatibility );
107     fTMVAReader[i]->AddVariable( "HadEnergyOverPt", &fMVAVar_MuHadEnergyOverPt );
108     fTMVAReader[i]->AddVariable( "EmEnergyOverPt", &fMVAVar_MuEmEnergyOverPt );
109     fTMVAReader[i]->AddVariable( "HadS9EnergyOverPt", &fMVAVar_MuHadS9EnergyOverPt );
110     fTMVAReader[i]->AddVariable( "EmS9EnergyOverPt", &fMVAVar_MuEmS9EnergyOverPt );
111     fTMVAReader[i]->AddVariable( "ChargedIso03OverPt", &fMVAVar_MuChargedIso03OverPt );
112     fTMVAReader[i]->AddVariable( "NeutralIso03OverPt", &fMVAVar_MuNeutralIso03OverPt );
113     fTMVAReader[i]->AddVariable( "ChargedIso04OverPt", &fMVAVar_MuChargedIso04OverPt );
114     fTMVAReader[i]->AddVariable( "NeutralIso04OverPt", &fMVAVar_MuNeutralIso04OverPt );
115     }
116 sixie 1.1
117     if (i==0) fTMVAReader[i]->BookMVA(fMethodname , Subdet0Pt10To14p5Weights );
118     if (i==1) fTMVAReader[i]->BookMVA(fMethodname , Subdet1Pt10To14p5Weights );
119     if (i==2) fTMVAReader[i]->BookMVA(fMethodname , Subdet0Pt14p5To20Weights );
120     if (i==3) fTMVAReader[i]->BookMVA(fMethodname , Subdet1Pt14p5To20Weights );
121     if (i==4) fTMVAReader[i]->BookMVA(fMethodname , Subdet0Pt20ToInfWeights );
122     if (i==5) fTMVAReader[i]->BookMVA(fMethodname , Subdet1Pt20ToInfWeights );
123    
124     }
125    
126     std::cout << "Muon ID MVA Initialization\n";
127     std::cout << "MethodName : " << fMethodname << " , type == " << type << std::endl;
128     std::cout << "Load weights file : " << Subdet0Pt10To14p5Weights << std::endl;
129     std::cout << "Load weights file : " << Subdet1Pt10To14p5Weights << std::endl;
130     std::cout << "Load weights file : " << Subdet0Pt14p5To20Weights << std::endl;
131     std::cout << "Load weights file : " << Subdet1Pt14p5To20Weights << std::endl;
132     std::cout << "Load weights file : " << Subdet0Pt20ToInfWeights << std::endl;
133     std::cout << "Load weights file : " << Subdet1Pt20ToInfWeights << std::endl;
134    
135     }
136    
137     //--------------------------------------------------------------------------------------------------
138     Double_t MuonIDMVA::MVAValue(Double_t MuPt , Double_t MuEta,
139 sixie 1.2 Double_t MuTkNchi2,
140     Double_t MuGlobalNchi2,
141     Double_t MuNValidHits,
142     Double_t MuNTrackerHits,
143     Double_t MuNPixelHits,
144     Double_t MuNMatches,
145     Double_t MuD0,
146     Double_t MuIP3d,
147     Double_t MuIP3dSig,
148     Double_t MuTrkKink,
149     Double_t MuSegmentCompatibility,
150     Double_t MuCaloCompatibility,
151     Double_t MuHadEnergyOverPt,
152     Double_t MuHoEnergyOverPt,
153     Double_t MuEmEnergyOverPt,
154     Double_t MuHadS9EnergyOverPt,
155     Double_t MuHoS9EnergyOverPt,
156     Double_t MuEmS9EnergyOverPt,
157     Double_t MuChargedIso03OverPt,
158     Double_t MuNeutralIso03OverPt,
159     Double_t MuChargedIso04OverPt,
160     Double_t MuNeutralIso04OverPt
161 sixie 1.1 ) {
162    
163     if (!fIsInitialized) {
164     std::cout << "Error: MuonIDMVA not properly initialized.\n";
165     return -9999;
166     }
167    
168     Int_t subdet = 0;
169     if (fabs(MuEta) < 1.479) subdet = 0;
170     else subdet = 1;
171     Int_t ptBin = 0;
172     if (MuPt > 14.5) ptBin = 1;
173     if (MuPt > 20.0) ptBin = 2;
174    
175    
176     //set all input variables
177     fMVAVar_MuTkNchi2 = MuTkNchi2;
178     fMVAVar_MuGlobalNchi2 = MuGlobalNchi2;
179     fMVAVar_MuNValidHits = MuNValidHits;
180     fMVAVar_MuNTrackerHits = MuNTrackerHits;
181     fMVAVar_MuNPixelHits = MuNPixelHits;
182     fMVAVar_MuNMatches = MuNMatches;
183     fMVAVar_MuD0 = MuD0;
184     fMVAVar_MuIP3d = MuIP3d;
185     fMVAVar_MuIP3dSig = MuIP3dSig;
186     fMVAVar_MuTrkKink = MuTrkKink;
187     fMVAVar_MuSegmentCompatibility = MuSegmentCompatibility;
188     fMVAVar_MuCaloCompatibility = MuCaloCompatibility;
189     fMVAVar_MuHadEnergyOverPt = MuHadEnergyOverPt;
190     fMVAVar_MuHoEnergyOverPt = MuHoEnergyOverPt;
191     fMVAVar_MuEmEnergyOverPt = MuEmEnergyOverPt;
192     fMVAVar_MuHadS9EnergyOverPt = MuHadS9EnergyOverPt;
193     fMVAVar_MuHoS9EnergyOverPt = MuHoS9EnergyOverPt;
194     fMVAVar_MuEmS9EnergyOverPt = MuEmS9EnergyOverPt;
195 sixie 1.2 fMVAVar_MuChargedIso03OverPt = MuChargedIso03OverPt;
196     fMVAVar_MuNeutralIso03OverPt = MuNeutralIso03OverPt;
197     fMVAVar_MuChargedIso04OverPt = MuChargedIso04OverPt;
198     fMVAVar_MuNeutralIso04OverPt = MuNeutralIso04OverPt;
199 sixie 1.1
200     Double_t mva = -9999;
201     TMVA::Reader *reader = 0;
202    
203     Int_t MVABin = -1;
204     if (subdet == 0 && ptBin == 0) MVABin = 0;
205     if (subdet == 1 && ptBin == 0) MVABin = 1;
206     if (subdet == 0 && ptBin == 1) MVABin = 2;
207     if (subdet == 1 && ptBin == 1) MVABin = 3;
208     if (subdet == 0 && ptBin == 2) MVABin = 4;
209     if (subdet == 1 && ptBin == 2) MVABin = 5;
210     assert(MVABin >= 0 && MVABin <= 5);
211     reader = fTMVAReader[MVABin];
212    
213     mva = reader->EvaluateMVA( fMethodname );
214    
215 sixie 1.5 Bool_t printdebug = kFALSE;
216 sixie 1.3 if (printdebug == kTRUE) {
217     std::cout << "Debug Muon MVA: "
218     << MuPt << " " << MuEta << " --> MVABin " << MVABin << " : "
219     << fMVAVar_MuTkNchi2 << " "
220     << fMVAVar_MuGlobalNchi2 << " "
221     << fMVAVar_MuNValidHits << " "
222     << fMVAVar_MuNTrackerHits << " "
223     << fMVAVar_MuNPixelHits << " "
224     << fMVAVar_MuNMatches << " "
225     << fMVAVar_MuD0 << " "
226     << fMVAVar_MuIP3d << " "
227     << fMVAVar_MuIP3dSig << " "
228     << fMVAVar_MuTrkKink << " "
229     << fMVAVar_MuSegmentCompatibility << " "
230     << fMVAVar_MuCaloCompatibility << " "
231     << fMVAVar_MuHadEnergyOverPt << " "
232     << fMVAVar_MuHoEnergyOverPt << " "
233     << fMVAVar_MuEmEnergyOverPt << " "
234     << fMVAVar_MuHadS9EnergyOverPt << " "
235     << fMVAVar_MuHoS9EnergyOverPt << " "
236     << fMVAVar_MuEmS9EnergyOverPt << " "
237     << fMVAVar_MuChargedIso03OverPt << " "
238     << fMVAVar_MuNeutralIso03OverPt << " "
239     << fMVAVar_MuChargedIso04OverPt << " "
240     << fMVAVar_MuNeutralIso04OverPt << " "
241     << " === : === "
242     << mva
243     << std::endl;
244     }
245    
246 sixie 1.1 return mva;
247     }
248    
249    
250    
251 sixie 1.4 //--------------------------------------------------------------------------------------------------
252     Double_t MuonIDMVA::MVAValue(const Muon *mu, const Vertex *vertex, MuonTools *fMuonTools,
253     const PFCandidateCol *PFCands,
254     const PileupEnergyDensityCol *PileupEnergyDensity) {
255 sixie 1.1
256 sixie 1.4 if (!fIsInitialized) {
257     std::cout << "Error: MuonIDMVA not properly initialized.\n";
258     return -9999;
259     }
260    
261     const Track *muTrk=0;
262     if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
263     else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
264 sixie 1.1
265 sixie 1.4 Double_t muNchi2 = 0.0;
266     if(mu->HasGlobalTrk()) { muNchi2 = mu->GlobalTrk()->RChi2(); }
267     else if(mu->HasStandaloneTrk()) { muNchi2 = mu->StandaloneTrk()->RChi2(); }
268     else if(mu->HasTrackerTrk()) { muNchi2 = mu->TrackerTrk()->RChi2(); }
269    
270     Double_t ChargedIso03 = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.1, 99999, 0.3, 0.0, 0.0);
271     Double_t NeutralIso03_05Threshold = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.0, 0.5, 0.3, 0.0, 0.0);
272     Double_t ChargedIso04 = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.1, 99999, 0.4, 0.0, 0.0);
273     Double_t NeutralIso04_05Threshold = IsolationTools::PFMuonIsolation(mu, PFCands, vertex, 0.0, 0.5, 0.4, 0.0, 0.0);
274    
275     Double_t Rho = 0;
276     if (!(TMath::IsNaN(PileupEnergyDensity->At(0)->Rho()) || isinf(PileupEnergyDensity->At(0)->Rho()))) Rho = PileupEnergyDensity->At(0)->Rho();
277    
278     Int_t subdet = 0;
279     if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
280     else subdet = 1;
281     Int_t ptBin = 0;
282     if (muTrk->Pt() > 14.5) ptBin = 1;
283     if (muTrk->Pt() > 20.0) ptBin = 2;
284    
285     //set all input variables
286     fMVAVar_MuTkNchi2 = muTrk->RChi2();
287     fMVAVar_MuGlobalNchi2 = muNchi2;
288     fMVAVar_MuNValidHits = mu->NValidHits();
289     fMVAVar_MuNTrackerHits = muTrk->NHits();
290     fMVAVar_MuNPixelHits = muTrk->NPixelHits();
291     fMVAVar_MuNMatches = mu->NMatches();
292     fMVAVar_MuD0 = muTrk->D0Corrected(*vertex);
293     fMVAVar_MuIP3d = mu->Ip3dPV();
294     fMVAVar_MuIP3dSig = mu->Ip3dPVSignificance();
295     fMVAVar_MuTrkKink = mu->TrkKink();
296     fMVAVar_MuSegmentCompatibility = fMuonTools->GetSegmentCompatability(mu);
297     fMVAVar_MuCaloCompatibility = fMuonTools->GetCaloCompatability(mu, kTRUE, kTRUE);
298     fMVAVar_MuHadEnergyOverPt = (mu->HadEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadEnergy,muTrk->Eta()))/muTrk->Pt();
299     fMVAVar_MuHoEnergyOverPt = (mu->HoEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHoEnergy,muTrk->Eta()))/muTrk->Pt();
300     fMVAVar_MuEmEnergyOverPt = (mu->EmEnergy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEmEnergy,muTrk->Eta()))/muTrk->Pt();
301     fMVAVar_MuHadS9EnergyOverPt = (mu->HadS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadS9Energy,muTrk->Eta()))/muTrk->Pt();
302     fMVAVar_MuHoS9EnergyOverPt = (mu->HoS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuHoS9Energy,muTrk->Eta()))/muTrk->Pt();
303     fMVAVar_MuEmS9EnergyOverPt = (mu->EmS9Energy() - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuEmS9Energy,muTrk->Eta()))/muTrk->Pt();
304     fMVAVar_MuChargedIso03OverPt = (ChargedIso03 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso03,muTrk->Eta()))/muTrk->Pt();
305     fMVAVar_MuNeutralIso03OverPt = (NeutralIso03_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03,muTrk->Eta()))/muTrk->Pt();
306     fMVAVar_MuChargedIso04OverPt = (ChargedIso04 - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuChargedIso04,muTrk->Eta()))/muTrk->Pt();
307     fMVAVar_MuNeutralIso04OverPt = (NeutralIso04_05Threshold - Rho*MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso04,muTrk->Eta()))/muTrk->Pt();
308    
309     Double_t mva = -9999;
310     TMVA::Reader *reader = 0;
311    
312     Int_t MVABin = -1;
313     if (subdet == 0 && ptBin == 0) MVABin = 0;
314     if (subdet == 1 && ptBin == 0) MVABin = 1;
315     if (subdet == 0 && ptBin == 1) MVABin = 2;
316     if (subdet == 1 && ptBin == 1) MVABin = 3;
317     if (subdet == 0 && ptBin == 2) MVABin = 4;
318     if (subdet == 1 && ptBin == 2) MVABin = 5;
319     assert(MVABin >= 0 && MVABin <= 5);
320     reader = fTMVAReader[MVABin];
321 sixie 1.1
322 sixie 1.4 mva = reader->EvaluateMVA( fMethodname );
323    
324 sixie 1.5 Bool_t printdebug = kFALSE;
325 sixie 1.4 if (printdebug == kTRUE) {
326     std::cout << "Debug Muon MVA: "
327     << mu->Pt() << " " << mu->Eta() << " " << mu->Phi() << " : "
328     << muTrk->Pt() << " " << muTrk->Eta() << " --> MVABin " << MVABin << " : "
329     << fMVAVar_MuTkNchi2 << " "
330     << fMVAVar_MuGlobalNchi2 << " "
331     << fMVAVar_MuNValidHits << " "
332     << fMVAVar_MuNTrackerHits << " "
333     << fMVAVar_MuNPixelHits << " "
334     << fMVAVar_MuNMatches << " "
335     << fMVAVar_MuD0 << " "
336     << fMVAVar_MuIP3d << " "
337     << fMVAVar_MuIP3dSig << " "
338     << fMVAVar_MuTrkKink << " "
339     << fMVAVar_MuSegmentCompatibility << " "
340     << fMVAVar_MuCaloCompatibility << " "
341     << fMVAVar_MuHadEnergyOverPt << " "
342     << fMVAVar_MuHoEnergyOverPt << " "
343     << fMVAVar_MuEmEnergyOverPt << " "
344     << fMVAVar_MuHadS9EnergyOverPt << " "
345     << fMVAVar_MuHoS9EnergyOverPt << " "
346     << fMVAVar_MuEmS9EnergyOverPt << " "
347     << fMVAVar_MuChargedIso03OverPt << " "
348     << fMVAVar_MuNeutralIso03OverPt << " "
349     << fMVAVar_MuChargedIso04OverPt << " "
350     << fMVAVar_MuNeutralIso04OverPt << " "
351     << " === : === "
352     << mva
353     << std::endl;
354     }
355 sixie 1.1
356 sixie 1.4 return mva;
357     }