ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/LJMet/MultivariateAnalysis/macros/TMVApplication.C
Revision: 1.2
Committed: Wed Jan 21 18:36:07 2009 UTC (16 years, 3 months ago) by kukartse
Content type: text/plain
Branch: MAIN
CVS Tags: gak022309, gak021209, gak040209, gak012809, V00-01-09, V00-01-08, V00-01-07, V00-01-06, V00-01-05, V00-01-04
Changes since 1.1: +65 -16 lines
Log Message:
after JTerm status update

File Contents

# User Rev Content
1 kukartse 1.1 /**********************************************************************************
2     * Project : TMVA - a Root-integrated toolkit for multivariate data analysis *
3     * Package : TMVA *
4     * Exectuable: TMVApplication *
5     * *
6     * This macro provides a simple example on how to use the trained classifiers *
7     * within an analysis module *
8     * *
9     * ------------------------------------------------------------------------------ *
10     * see also the alternative (slightly faster) way to retrieve the MVA values in *
11     * examples/TMVApplicationAlternative.cxx *
12     * ------------------------------------------------------------------------------ *
13     **********************************************************************************/
14    
15 kukartse 1.2 #include <iostream>
16    
17     #include "TCut.h"
18     #include "TFile.h"
19     #include "TSystem.h"
20     #include "TTree.h"
21     #include "TStopwatch.h"
22     // requires links
23     #include "TMVA/Factory.h"
24     #include "TMVA/Tools.h"
25     #include "TMVA/Config.h"
26     #include "TMVA/Reader.h"
27    
28 kukartse 1.1 // ---------------------------------------------------------------
29     // choose MVA methods to be applied
30     Bool_t Use_Cuts = 0;
31     Bool_t Use_CutsD = 0;
32     Bool_t Use_CutsGA = 1;
33     Bool_t Use_Likelihood = 1;
34     Bool_t Use_LikelihoodD = 0; // the "D" extension indicates decorrelated input variables (see option strings)
35     Bool_t Use_LikelihoodPCA = 1; // the "PCA" extension indicates PCA-transformed input variables (see option strings)
36     Bool_t Use_PDERS = 1;
37     Bool_t Use_PDERSD = 0;
38     Bool_t Use_PDERSPCA = 0;
39     Bool_t Use_KNN = 1;
40     Bool_t Use_HMatrix = 1;
41     Bool_t Use_Fisher = 1;
42     Bool_t Use_FDA_GA = 0;
43     Bool_t Use_FDA_MT = 1;
44     Bool_t Use_MLP = 1; // this is the recommended ANN
45     Bool_t Use_CFMlpANN = 0;
46     Bool_t Use_TMlpANN = 0;
47     Bool_t Use_SVM_Gauss = 1;
48     Bool_t Use_SVM_Poly = 0;
49     Bool_t Use_SVM_Lin = 0;
50     Bool_t Use_BDT = 1;
51     Bool_t Use_BDTD = 0;
52     Bool_t Use_RuleFit = 1;
53     // ---------------------------------------------------------------
54    
55     void TMVApplication( TString myMethodList = "" )
56     {
57     cout << endl;
58     cout << "==> start TMVApplication" << endl;
59    
60     if (myMethodList != "") {
61     Use_CutsGA = Use_CutsD = Use_Cuts
62     = Use_LikelihoodPCA = Use_LikelihoodD = Use_Likelihood
63     = Use_PDERSPCA = Use_PDERSD = Use_PDERS
64     = Use_KNN
65     = Use_MLP = Use_CFMlpANN = Use_TMlpANN
66     = Use_HMatrix = Use_Fisher = Use_BDTD = Use_BDT = Use_RuleFit
67     = Use_SVM_Gauss = Use_SVM_Poly = Use_SVM_Lin
68     = Use_FDA_GA = Use_FDA_MT
69     = 0;
70    
71     TList* mlist = TMVA::Tools::ParseFormatLine( myMethodList, " :," );
72    
73     if (mlist->FindObject( "Cuts" ) != 0) Use_Cuts = 1;
74     if (mlist->FindObject( "CutsD" ) != 0) Use_CutsD = 1;
75     if (mlist->FindObject( "CutsGA" ) != 0) Use_CutsGA = 1;
76     if (mlist->FindObject( "Likelihood" ) != 0) Use_Likelihood = 1;
77     if (mlist->FindObject( "LikelihoodD" ) != 0) Use_LikelihoodD = 1;
78     if (mlist->FindObject( "LikelihoodPCA" ) != 0) Use_LikelihoodPCA = 1;
79     if (mlist->FindObject( "PDERS" ) != 0) Use_PDERS = 1;
80     if (mlist->FindObject( "PDERSD" ) != 0) Use_PDERSD = 1;
81     if (mlist->FindObject( "PDERSPCA" ) != 0) Use_PDERSPCA = 1;
82     if (mlist->FindObject( "KNN" ) != 0) Use_KNN = 1;
83     if (mlist->FindObject( "HMatrix" ) != 0) Use_HMatrix = 1;
84     if (mlist->FindObject( "Fisher" ) != 0) Use_Fisher = 1;
85     if (mlist->FindObject( "MLP" ) != 0) Use_MLP = 1;
86     if (mlist->FindObject( "CFMlpANN" ) != 0) Use_CFMlpANN = 1;
87     if (mlist->FindObject( "TMlpANN" ) != 0) Use_TMlpANN = 1;
88     if (mlist->FindObject( "BDTD" ) != 0) Use_BDTD = 1;
89     if (mlist->FindObject( "BDT" ) != 0) Use_BDT = 1;
90     if (mlist->FindObject( "RuleFit" ) != 0) Use_RuleFit = 1;
91     if (mlist->FindObject( "SVM_Gauss" ) != 0) Use_SVM_Gauss = 1;
92     if (mlist->FindObject( "SVM_Poly" ) != 0) Use_SVM_Poly = 1;
93     if (mlist->FindObject( "SVM_Lin" ) != 0) Use_SVM_Lin = 1;
94     if (mlist->FindObject( "FDA_MT" ) != 0) Use_FDA_MT = 1;
95     if (mlist->FindObject( "FDA_GA" ) != 0) Use_FDA_GA = 1;
96    
97     delete mlist;
98     }
99    
100     //
101     // create the Reader object
102     //
103     TMVA::Reader *reader = new TMVA::Reader("!Color");
104    
105     // create a set of variables and declare them to the reader
106     // - the variable names must corresponds in name and type to
107     // those given in the weight file(s) that you use
108 kukartse 1.2 Float_t _aplanarity;
109     Float_t _getHt3;
110     Float_t _ktMinPrime;
111     Float_t _DphiJMET;
112     Float_t _dPhiLMet;
113     Float_t _getMt;
114     reader->AddVariable("aplanarity", &_aplanarity);
115     reader->AddVariable("getHt3", &_getHt3);
116     reader->AddVariable("ktMinPrime", &_ktMinPrime);
117     reader->AddVariable("DphiJMET", &_DphiJMET);
118     //reader->AddVariable("dPhiLMet", &_dPhiLMet);
119     reader->AddVariable("getMt", &_getMt);
120 kukartse 1.1
121     //
122     // book the MVA methods
123     //
124     string dir = "weights/";
125     string prefix = "TMVAnalysis";
126    
127     if (Use_Cuts) reader->BookMVA( "Cuts method", dir + prefix + "_Cuts.weights.txt" );
128     if (Use_CutsD) reader->BookMVA( "CutsD method", dir + prefix + "_CutsD.weights.txt" );
129     if (Use_CutsGA) reader->BookMVA( "CutsGA method", dir + prefix + "_CutsGA.weights.txt" );
130     if (Use_Likelihood) reader->BookMVA( "Likelihood method", dir + prefix + "_Likelihood.weights.txt" );
131     if (Use_LikelihoodD) reader->BookMVA( "LikelihoodD method", dir + prefix + "_LikelihoodD.weights.txt" );
132     if (Use_LikelihoodPCA) reader->BookMVA( "LikelihoodPCA method", dir + prefix + "_LikelihoodPCA.weights.txt" );
133     if (Use_PDERS) reader->BookMVA( "PDERS method", dir + prefix + "_PDERS.weights.txt" );
134     if (Use_PDERSD) reader->BookMVA( "PDERSD method", dir + prefix + "_PDERSD.weights.txt" );
135     if (Use_PDERSPCA) reader->BookMVA( "PDERSPCA method", dir + prefix + "_PDERSPCA.weights.txt" );
136     if (Use_KNN) reader->BookMVA( "KNN method", dir + prefix + "_KNN.weights.txt" );
137     if (Use_HMatrix) reader->BookMVA( "HMatrix method", dir + prefix + "_HMatrix.weights.txt" );
138     if (Use_Fisher) reader->BookMVA( "Fisher method", dir + prefix + "_Fisher.weights.txt" );
139     if (Use_MLP) reader->BookMVA( "MLP method", dir + prefix + "_MLP.weights.txt" );
140     if (Use_CFMlpANN) reader->BookMVA( "CFMlpANN method", dir + prefix + "_CFMlpANN.weights.txt" );
141     if (Use_TMlpANN) reader->BookMVA( "TMlpANN method", dir + prefix + "_TMlpANN.weights.txt" );
142     if (Use_BDT) reader->BookMVA( "BDT method", dir + prefix + "_BDT.weights.txt" );
143     if (Use_BDTD) reader->BookMVA( "BDTD method", dir + prefix + "_BDTD.weights.txt" );
144     if (Use_RuleFit) reader->BookMVA( "RuleFit method", dir + prefix + "_RuleFitTMVA.weights.txt" );
145     if (Use_SVM_Gauss) reader->BookMVA( "SVM_Gauss method", dir + prefix + "_SVM_Gauss.weights.txt" );
146     if (Use_SVM_Poly) reader->BookMVA( "SVM_Poly method", dir + prefix + "_SVM_Poly.weights.txt" );
147     if (Use_SVM_Lin) reader->BookMVA( "SVM_Lin method", dir + prefix + "_SVM_Lin.weights.txt" );
148     if (Use_FDA_MT) reader->BookMVA( "FDA_MT method", dir + prefix + "_FDA_MT.weights.txt" );
149     if (Use_FDA_GA) reader->BookMVA( "FDA_GA method", dir + prefix + "_FDA_GA.weights.txt" );
150    
151     // book output histograms
152     UInt_t nbin = 100;
153 kukartse 1.2 TFile *target = new TFile( "TMVApp.root","RECREATE" );
154 kukartse 1.1 TH1F *histLk, *histLkD, *histLkPCA, *histPD, *histPDD, *histPDPCA, *histKNN, *histHm, *histFi;
155     TH1F *histNn, *histNnC, *histNnT, *histBdt, *histBdtD, *histRf;
156     TH1F *histSVMG, *histSVMP, *histSVML;
157     TH1F *histFDAMT, *histFDAGA;
158    
159     if (Use_Likelihood) histLk = new TH1F( "MVA_Likelihood", "MVA_Likelihood", nbin, 0, 1 );
160 kukartse 1.2 TTree * _tree = new TTree("classifier","classifier");
161     Double_t _ll = -1.0;
162     TBranch * b_ll = _tree->Branch("MVA_Likelihood", &_ll, "classifier/D");
163 kukartse 1.1 if (Use_LikelihoodD) histLkD = new TH1F( "MVA_LikelihoodD", "MVA_LikelihoodD", nbin, 0.000001, 0.9999 );
164     if (Use_LikelihoodPCA) histLkPCA = new TH1F( "MVA_LikelihoodPCA", "MVA_LikelihoodPCA", nbin, 0, 1 );
165     if (Use_PDERS) histPD = new TH1F( "MVA_PDERS", "MVA_PDERS", nbin, 0, 1 );
166     if (Use_PDERSD) histPDD = new TH1F( "MVA_PDERSD", "MVA_PDERSD", nbin, 0, 1 );
167     if (Use_PDERSPCA) histPDPCA = new TH1F( "MVA_PDERSPCA", "MVA_PDERSPCA", nbin, 0, 1 );
168     if (Use_KNN) histKNN = new TH1F( "MVA_KNN", "MVA_KNN", nbin, 0, 1 );
169     if (Use_HMatrix) histHm = new TH1F( "MVA_HMatrix", "MVA_HMatrix", nbin, -0.95, 1.55 );
170     if (Use_Fisher) histFi = new TH1F( "MVA_Fisher", "MVA_Fisher", nbin, -4, 4 );
171     if (Use_MLP) histNn = new TH1F( "MVA_MLP", "MVA_MLP", nbin, -0.25, 1.5 );
172     if (Use_CFMlpANN) histNnC = new TH1F( "MVA_CFMlpANN", "MVA_CFMlpANN", nbin, 0, 1 );
173     if (Use_TMlpANN) histNnT = new TH1F( "MVA_TMlpANN", "MVA_TMlpANN", nbin, -1.3, 1.3 );
174     if (Use_BDT) histBdt = new TH1F( "MVA_BDT", "MVA_BDT", nbin, -0.8, 0.8 );
175     if (Use_BDTD) histBdtD = new TH1F( "MVA_BDTD", "MVA_BDTD", nbin, -0.4, 0.6 );
176     if (Use_RuleFit) histRf = new TH1F( "MVA_RuleFitTMVA", "MVA_RuleFitTMVA", nbin, -2.0, 2.0 );
177     if (Use_SVM_Gauss) histSVMG = new TH1F( "MVA_SVM_Gauss", "MVA_SVM_Gauss", nbin, 0.0, 1.0 );
178     if (Use_SVM_Poly) histSVMP = new TH1F( "MVA_SVM_Poly", "MVA_SVM_Poly", nbin, 0.0, 1.0 );
179     if (Use_SVM_Lin) histSVML = new TH1F( "MVA_SVM_Lin", "MVA_SVM_Lin", nbin, 0.0, 1.0 );
180     if (Use_FDA_MT) histFDAMT = new TH1F( "MVA_FDA_MT", "MVA_FDA_MT", nbin, -2.0, 3.0 );
181     if (Use_FDA_GA) histFDAGA = new TH1F( "MVA_FDA_GA", "MVA_FDA_GA", nbin, -2.0, 3.0 );
182    
183     // book examsple histogram for probability (the other methods are done similarly)
184     TH1F *probHistFi, *rarityHistFi;
185     if (Use_Fisher) {
186     probHistFi = new TH1F( "PROBA_MVA_Fisher", "PROBA_MVA_Fisher", nbin, 0, 1 );
187     rarityHistFi = new TH1F( "RARITY_MVA_Fisher", "RARITY_MVA_Fisher", nbin, 0, 1 );
188     }
189    
190     // Prepare input tree (this must be replaced by your data source)
191     // in this example, there is a toy tree with signal and one with background events
192     // we'll later on use only the "signal" events for the test in this example.
193     //
194     TFile *input(0);
195 kukartse 1.2 //TString fname = "/uscms_data/d1/lpcljm/MVA/Summer08/fake_data/tmva_fake_data_100pb-summer08-09jan2009.root";
196     TString fname = "/uscms_data/d1/lpcljm/MVA/Summer08/fake_data/tmva_fake_data_100pb-summer08-09jan2009.root";
197 kukartse 1.1 if (!gSystem->AccessPathName( fname )) {
198     // first we try to find tmva_example.root in the local directory
199     cout << "--- accessing data file: " << fname << endl;
200     input = TFile::Open( fname );
201     }
202     else {
203     // second we try accessing the file via the web from
204     // http://root.cern.ch/files/tmva_example.root
205     cout << "--- accessing tmva_example.root file from http://root.cern.ch/files" << endl;
206     cout << "--- for faster startup you may consider downloading it into you local directory" << endl;
207     input = TFile::Open("http://root.cern.ch/files/tmva_example.root");
208     }
209    
210     if (!input) {
211     std::cout << "ERROR: could not open data file: " << fname << std::endl;
212     exit(1);
213     }
214    
215     //
216     // prepare the tree
217     // - here the variable names have to corresponds to your tree
218     // - you can use the same variables as above which is slightly faster,
219     // but of course you can use different ones and copy the values inside the event loop
220     //
221 kukartse 1.2 TTree* theTree = (TTree*)input->Get("data");
222 kukartse 1.1 cout << "--- select signal sample" << endl;
223 kukartse 1.2 //Float_t userVar1, userVar2;
224     Double_t aplanarity;
225     Double_t getHt3;
226     Double_t ktMinPrime;
227     Double_t DphiJMET;
228     Double_t dPhiLMet;
229     Double_t getMt;
230     theTree->SetBranchAddress( "aplanarity", &aplanarity );
231     theTree->SetBranchAddress( "getHt3", &getHt3 );
232     theTree->SetBranchAddress( "ktMinPrime", &ktMinPrime );
233     theTree->SetBranchAddress( "DphiJMET", &DphiJMET );
234     //theTree->SetBranchAddress( "dPhiLMet", &dPhiLMet );
235     theTree->SetBranchAddress( "getMt", &getMt );
236 kukartse 1.1
237     // efficiency calculator for cut method
238     Int_t nSelCuts = 0, nSelCutsD = 0, nSelCutsGA = 0;
239     Double_t effS = 0.7;
240    
241     cout << "--- processing: " << theTree->GetEntries() << " events" << endl;
242     TStopwatch sw;
243     sw.Start();
244     for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
245    
246     if (ievt%1000 == 0)
247     cout << "--- ... processing event: " << ievt << endl;
248    
249     theTree->GetEntry(ievt);
250    
251 kukartse 1.2 //var1 = userVar1 + userVar2;
252     //var2 = userVar1 - userVar2;
253     _aplanarity = (Float_t)aplanarity;
254     _getHt3 = (Float_t)getHt3;
255     _ktMinPrime = (Float_t)ktMinPrime;
256     _DphiJMET = (Float_t)DphiJMET;
257     //_dPhiLMet = (Float_t)dPhiLMet;
258     _getMt = (Float_t)getMt;
259 kukartse 1.1
260     //
261     // return the MVAs and fill to histograms
262     //
263     if (Use_Cuts) {
264     // Cuts is a special case: give the desired signal efficienciy
265     Bool_t passed = reader->EvaluateMVA( "Cuts method", effS );
266     if (passed) nSelCuts++;
267     }
268     if (Use_CutsD) {
269     // Cuts is a special case: give the desired signal efficienciy
270     Bool_t passed = reader->EvaluateMVA( "CutsD method", effS );
271     if (passed) nSelCutsD++;
272     }
273     if (Use_CutsGA) {
274     // Cuts is a special case: give the desired signal efficienciy
275     Bool_t passed = reader->EvaluateMVA( "CutsGA method", effS );
276     if (passed) nSelCutsGA++;
277     }
278    
279 kukartse 1.2 if (Use_Likelihood ){
280     _ll = (Double_t)(reader->EvaluateMVA("Likelihood method") );
281     histLk ->Fill( reader->EvaluateMVA( "Likelihood method" ) );
282     _tree ->Fill();
283     }
284 kukartse 1.1 if (Use_LikelihoodD ) histLkD ->Fill( reader->EvaluateMVA( "LikelihoodD method" ) );
285     if (Use_LikelihoodPCA) histLkPCA ->Fill( reader->EvaluateMVA( "LikelihoodPCA method" ) );
286     if (Use_PDERS ) histPD ->Fill( reader->EvaluateMVA( "PDERS method" ) );
287     if (Use_PDERSD ) histPDD ->Fill( reader->EvaluateMVA( "PDERSD method" ) );
288     if (Use_PDERSPCA ) histPDPCA ->Fill( reader->EvaluateMVA( "PDERSPCA method" ) );
289     if (Use_KNN ) histKNN ->Fill( reader->EvaluateMVA( "KNN method" ) );
290     if (Use_HMatrix ) histHm ->Fill( reader->EvaluateMVA( "HMatrix method" ) );
291     if (Use_Fisher ) histFi ->Fill( reader->EvaluateMVA( "Fisher method" ) );
292     if (Use_MLP ) histNn ->Fill( reader->EvaluateMVA( "MLP method" ) );
293     if (Use_CFMlpANN ) histNnC ->Fill( reader->EvaluateMVA( "CFMlpANN method" ) );
294     if (Use_TMlpANN ) histNnT ->Fill( reader->EvaluateMVA( "TMlpANN method" ) );
295     if (Use_BDT ) histBdt ->Fill( reader->EvaluateMVA( "BDT method" ) );
296     if (Use_BDTD ) histBdtD ->Fill( reader->EvaluateMVA( "BDTD method" ) );
297     if (Use_RuleFit ) histRf ->Fill( reader->EvaluateMVA( "RuleFit method" ) );
298     if (Use_SVM_Gauss ) histSVMG ->Fill( reader->EvaluateMVA( "SVM_Gauss method" ) );
299     if (Use_SVM_Poly ) histSVMP ->Fill( reader->EvaluateMVA( "SVM_Poly method" ) );
300     if (Use_SVM_Lin ) histSVML ->Fill( reader->EvaluateMVA( "SVM_Lin method" ) );
301     if (Use_FDA_MT ) histFDAMT ->Fill( reader->EvaluateMVA( "FDA_MT method" ) );
302     if (Use_FDA_GA ) histFDAGA ->Fill( reader->EvaluateMVA( "FDA_GA method" ) );
303    
304     // retrieve probability instead of MVA output
305     if (Use_Fisher ) {
306     probHistFi ->Fill( reader->GetProba ( "Fisher method" ) );
307     rarityHistFi->Fill( reader->GetRarity( "Fisher method" ) );
308     }
309     }
310     sw.Stop();
311     cout << "--- end of event loop: "; sw.Print();
312     // get elapsed time
313     if (Use_Cuts) cout << "--- efficiency for Cuts method : " << double(nSelCuts)/theTree->GetEntries()
314     << " (for a required signal efficiency of " << effS << ")" << endl;
315     if (Use_CutsD) cout << "--- efficiency for CutsD method : " << double(nSelCutsD)/theTree->GetEntries()
316     << " (for a required signal efficiency of " << effS << ")" << endl;
317     if (Use_CutsGA) cout << "--- efficiency for CutsGA method: " << double(nSelCutsGA)/theTree->GetEntries()
318     << " (for a required signal efficiency of " << effS << ")" << endl;
319    
320 kukartse 1.2 /*************************
321 kukartse 1.1 //
322     // write histograms
323     //
324     TFile *target = new TFile( "TMVApp.root","RECREATE" );
325     if (Use_Likelihood ) histLk ->Write();
326     if (Use_LikelihoodD ) histLkD ->Write();
327     if (Use_LikelihoodPCA) histLkPCA ->Write();
328     if (Use_PDERS ) histPD ->Write();
329     if (Use_PDERSD ) histPDD ->Write();
330     if (Use_PDERSPCA ) histPDPCA ->Write();
331     if (Use_KNN ) histKNN ->Write();
332     if (Use_HMatrix ) histHm ->Write();
333     if (Use_Fisher ) histFi ->Write();
334     if (Use_MLP ) histNn ->Write();
335     if (Use_CFMlpANN ) histNnC ->Write();
336     if (Use_TMlpANN ) histNnT ->Write();
337     if (Use_BDT ) histBdt ->Write();
338     if (Use_BDTD ) histBdtD ->Write();
339     if (Use_RuleFit ) histRf ->Write();
340     if (Use_SVM_Gauss ) histSVMG ->Write();
341     if (Use_SVM_Poly ) histSVMP ->Write();
342     if (Use_SVM_Lin ) histSVML ->Write();
343     if (Use_FDA_MT ) histFDAMT ->Write();
344     if (Use_FDA_GA ) histFDAGA ->Write();
345    
346     if (Use_Fisher ) { probHistFi->Write(); rarityHistFi->Write(); }
347 kukartse 1.2
348 kukartse 1.1 target->Close();
349 kukartse 1.2 *******************************/
350    
351     target->cd();
352     target->Write();
353     delete target;
354 kukartse 1.1
355     cout << "--- created root file: \"TMVApp.root\" containing the MVA output histograms" << endl;
356    
357     delete reader;
358    
359     cout << "==> TMVApplication is done!" << endl << endl;
360     }