ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MuonIDMVA.cc
Revision: 1.3
Committed: Wed Jan 4 10:32:05 2012 UTC (13 years, 4 months ago) by sixie
Content type: text/plain
Branch: MAIN
Changes since 1.2: +31 -0 lines
Log Message:
add debug printouts

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.3 Bool_t printdebug = kFALSE;
216     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     // //--------------------------------------------------------------------------------------------------
252     // Double_t MuonIDMVA::MVAValue(const Muon *ele, const Vertex *vertex) {
253    
254     // if (!fIsInitialized) {
255     // std::cout << "Error: MuonIDMVA not properly initialized.\n";
256     // return -9999;
257     // }
258    
259     // Int_t subdet = 0;
260     // if (ele->SCluster()->AbsEta() < 1.0) subdet = 0;
261     // else if (ele->SCluster()->AbsEta() < 1.479) subdet = 1;
262     // else subdet = 2;
263     // Int_t ptBin = 0;
264     // if (ele->Pt() > 20.0) ptBin = 1;
265    
266     // //set all input variables
267     // fMVAVar_EleSigmaIEtaIEta = ele->CoviEtaiEta() ;
268     // fMVAVar_EleDEtaIn = ele->DeltaEtaSuperClusterTrackAtVtx();
269     // fMVAVar_EleDPhiIn = ele->DeltaPhiSuperClusterTrackAtVtx();
270     // fMVAVar_EleHoverE = ele->HadronicOverEm();
271     // fMVAVar_EleD0 = ele->BestTrk()->D0Corrected(*vertex);
272     // fMVAVar_EleDZ = ele->BestTrk()->DzCorrected(*vertex);
273     // fMVAVar_EleFBrem = ele->FBrem();
274     // fMVAVar_EleEOverP = ele->ESuperClusterOverP();
275     // fMVAVar_EleESeedClusterOverPout = ele->ESeedClusterOverPout();
276     // if (!TMath::IsNaN(ele->SCluster()->Seed()->CoviPhiiPhi())) fMVAVar_EleSigmaIPhiIPhi = TMath::Sqrt(ele->SCluster()->Seed()->CoviPhiiPhi());
277     // else fMVAVar_EleSigmaIPhiIPhi = ele->CoviEtaiEta();
278     // fMVAVar_EleNBrem = ele->NumberOfClusters() - 1;
279     // fMVAVar_EleOneOverEMinusOneOverP = (1.0/(ele->SCluster()->Energy())) - 1.0 / ele->BestTrk()->P();
280     // fMVAVar_EleESeedClusterOverPIn = ele->ESeedClusterOverPIn();
281     // fMVAVar_EleIP3d = ele->Ip3dPV();
282     // fMVAVar_EleIP3dSig = ele->Ip3dPVSignificance();
283    
284     // Double_t mva = -9999;
285     // TMVA::Reader *reader = 0;
286     // Int_t MVABin = -1;
287     // if (subdet == 0 && ptBin == 0) MVABin = 0;
288     // if (subdet == 1 && ptBin == 0) MVABin = 1;
289     // if (subdet == 2 && ptBin == 0) MVABin = 2;
290     // if (subdet == 0 && ptBin == 1) MVABin = 3;
291     // if (subdet == 1 && ptBin == 1) MVABin = 4;
292     // if (subdet == 2 && ptBin == 1) MVABin = 5;
293     // assert(MVABin >= 0 && MVABin <= 5);
294     // reader = fTMVAReader[MVABin];
295    
296     // mva = reader->EvaluateMVA( fMethodname );
297    
298     // return mva;
299     // }