ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/benhoob/HWW/Classify_HWW.C
Revision: 1.3
Committed: Mon Feb 14 12:49:15 2011 UTC (14 years, 3 months ago) by benhoob
Content type: text/plain
Branch: MAIN
Changes since 1.2: +36 -19 lines
Log Message:
Updates to selection

File Contents

# User Rev Content
1 benhoob 1.1 /**********************************************************************************
2     * Project : TMVA - a Root-integrated toolkit for multivariate data analysis *
3     * Package : TMVA *
4     * Exectuable: TMVAClassificationApplication *
5     * *
6     * This macro provides a simple example on how to use the trained classifiers *
7     * within an analysis module *
8     **********************************************************************************/
9    
10     #include <cstdlib>
11     #include <vector>
12     #include <iostream>
13     #include <map>
14     #include <string>
15    
16     #include "TFile.h"
17     #include "TTree.h"
18     #include "TString.h"
19     #include "TSystem.h"
20     #include "TROOT.h"
21     #include "TStopwatch.h"
22     #include "TChain.h"
23    
24     #include "TMVAGui.C"
25    
26     #if not defined(__CINT__) || defined(__MAKECINT__)
27     #include "TMVA/Tools.h"
28     #include "TMVA/Reader.h"
29     #include "TMVA/MethodCuts.h"
30     #endif
31    
32     using namespace TMVA;
33    
34     void Classify_HWW( TString myMethodList = "" )
35     {
36     #ifdef __CINT__
37     gROOT->ProcessLine( ".O0" ); // turn off optimization in CINT
38     #endif
39    
40     //---------------------------------------------------------------
41    
42     // This loads the library
43     TMVA::Tools::Instance();
44    
45     // Default MVA methods to be trained + tested
46     std::map<std::string,int> Use;
47    
48     // --- Cut optimisation
49     Use["Cuts"] = 1;
50     Use["CutsD"] = 1;
51     Use["CutsPCA"] = 0;
52     Use["CutsGA"] = 0;
53     Use["CutsSA"] = 0;
54     //
55     // --- 1-dimensional likelihood ("naive Bayes estimator")
56     Use["Likelihood"] = 1;
57     Use["LikelihoodD"] = 0; // the "D" extension indicates decorrelated input variables (see option strings)
58     Use["LikelihoodPCA"] = 1; // the "PCA" extension indicates PCA-transformed input variables (see option strings)
59     Use["LikelihoodKDE"] = 0;
60     Use["LikelihoodMIX"] = 0;
61     //
62     // --- Mutidimensional likelihood and Nearest-Neighbour methods
63     Use["PDERS"] = 1;
64     Use["PDERSD"] = 0;
65     Use["PDERSPCA"] = 0;
66     Use["PDEFoam"] = 1;
67     Use["PDEFoamBoost"] = 0; // uses generalised MVA method boosting
68     Use["KNN"] = 1; // k-nearest neighbour method
69     //
70     // --- Linear Discriminant Analysis
71     Use["LD"] = 1; // Linear Discriminant identical to Fisher
72     Use["Fisher"] = 0;
73     Use["FisherG"] = 0;
74     Use["BoostedFisher"] = 0; // uses generalised MVA method boosting
75     Use["HMatrix"] = 0;
76     //
77     // --- Function Discriminant analysis
78     Use["FDA_GA"] = 1; // minimisation of user-defined function using Genetics Algorithm
79     Use["FDA_SA"] = 0;
80     Use["FDA_MC"] = 0;
81     Use["FDA_MT"] = 0;
82     Use["FDA_GAMT"] = 0;
83     Use["FDA_MCMT"] = 0;
84     //
85     // --- Neural Networks (all are feed-forward Multilayer Perceptrons)
86     Use["MLP"] = 0; // Recommended ANN
87     Use["MLPBFGS"] = 0; // Recommended ANN with optional training method
88     Use["MLPBNN"] = 1; // Recommended ANN with BFGS training method and bayesian regulator
89     Use["CFMlpANN"] = 0; // Depreciated ANN from ALEPH
90     Use["TMlpANN"] = 0; // ROOT's own ANN
91     //
92     // --- Support Vector Machine
93     Use["SVM"] = 1;
94     //
95     // --- Boosted Decision Trees
96     Use["BDT"] = 1; // uses Adaptive Boost
97     Use["BDTG"] = 0; // uses Gradient Boost
98     Use["BDTB"] = 0; // uses Bagging
99     Use["BDTD"] = 0; // decorrelation + Adaptive Boost
100     //
101     // --- Friedman's RuleFit method, ie, an optimised series of cuts ("rules")
102     Use["RuleFit"] = 1;
103     // ---------------------------------------------------------------
104     Use["Plugin"] = 0;
105     Use["Category"] = 0;
106     Use["SVM_Gauss"] = 0;
107     Use["SVM_Poly"] = 0;
108     Use["SVM_Lin"] = 0;
109    
110     std::cout << std::endl;
111     std::cout << "==> Start TMVAClassificationApplication" << std::endl;
112    
113     // Select methods (don't look at this code - not of interest)
114     if (myMethodList != "") {
115     for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;
116    
117     std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' );
118     for (UInt_t i=0; i<mlist.size(); i++) {
119     std::string regMethod(mlist[i]);
120    
121     if (Use.find(regMethod) == Use.end()) {
122     std::cout << "Method \"" << regMethod
123     << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
124     for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
125     std::cout << it->first << " ";
126     }
127     std::cout << std::endl;
128     return;
129     }
130     Use[regMethod] = 1;
131     }
132     }
133    
134     // --------------------------------------------------------------------------------------------------
135    
136     //-----------------------------------
137     //select samples to run over
138     //-----------------------------------
139    
140     vector<char*> samples;
141     samples.push_back("WWTo2L2Nu");
142 benhoob 1.2 samples.push_back("GluGluToWWTo4L");
143     samples.push_back("WZ");
144     samples.push_back("ZZ");
145     samples.push_back("TTJets");
146     samples.push_back("tW");
147     samples.push_back("WJetsToLNu");
148     samples.push_back("DY");
149     samples.push_back("Higgs130");
150 benhoob 1.1 // samples.push_back("Higgs160");
151     // samples.push_back("Higgs200");
152    
153     const unsigned int nsamples = samples.size();
154    
155     for( unsigned int i = 0 ; i < nsamples ; ++i ){
156    
157     // --- Create the Reader object
158    
159     TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );
160    
161     // Create a set of variables and declare them to the reader
162     // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
163     // Float_t var1, var2;
164     // Float_t var3, var4;
165     // reader->AddVariable( "myvar1 := var1+var2", &var1 );
166     // reader->AddVariable( "myvar2 := var1-var2", &var2 );
167     // reader->AddVariable( "var3", &var3 );
168     // reader->AddVariable( "var4", &var4 );
169    
170     Float_t lephard_pt;
171     Float_t lepsoft_pt;
172     Float_t dil_dphi;
173     Float_t dil_mass;
174    
175     reader->AddVariable( "lephard_pt" , &lephard_pt );
176     reader->AddVariable( "lepsoft_pt" , &lepsoft_pt );
177     reader->AddVariable( "dil_dphi" , &dil_dphi );
178     reader->AddVariable( "dil_mass" , &dil_mass );
179    
180     // Spectator variables declared in the training have to be added to the reader, too
181     // Float_t spec1,spec2;
182     // reader->AddSpectator( "spec1 := var1*2", &spec1 );
183     // reader->AddSpectator( "spec2 := var1*3", &spec2 );
184    
185     Float_t Category_cat1, Category_cat2, Category_cat3;
186     if (Use["Category"]){
187     // Add artificial spectators for distinguishing categories
188     // reader->AddSpectator( "Category_cat1 := var3<=0", &Category_cat1 );
189     // reader->AddSpectator( "Category_cat2 := (var3>0)&&(var4<0)", &Category_cat2 );
190     // reader->AddSpectator( "Category_cat3 := (var3>0)&&(var4>=0)", &Category_cat3 );
191     }
192    
193     // --- Book the MVA methods
194    
195 benhoob 1.3 //--------------------------------------------------------------------------------------
196     // tell Classify_HWW where to find the weights dir, which contains the trained MVA's.
197     // In this example, the weights dir is located at [path]/[dir]
198     // and the output root file is written to [path]/[output]
199     //--------------------------------------------------------------------------------------
200    
201 benhoob 1.1 //TString dir = "weights/";
202 benhoob 1.2 //TString path = "Trainings/H130_WWTo2L2Nu/";
203     //TString path = "Trainings/H130_WWTo2L2Nu_WJetsToLNu/";
204     TString path = "Trainings/H130_allbkg/";
205 benhoob 1.1 TString dir = path + "weights/";
206     TString outdir = path + "output/";
207    
208     TString prefix = "TMVAClassification";
209    
210     // Book method(s)
211     for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
212     if (it->second) {
213     TString methodName = TString(it->first) + TString(" method");
214     TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
215     reader->BookMVA( methodName, weightfile );
216     }
217     }
218    
219     // Book output histograms
220 benhoob 1.2 UInt_t nbin = 1000;
221 benhoob 1.1 TH1F *histLk(0), *histLkD(0), *histLkPCA(0), *histLkKDE(0), *histLkMIX(0), *histPD(0), *histPDD(0);
222     TH1F *histPDPCA(0), *histPDEFoam(0), *histPDEFoamErr(0), *histPDEFoamSig(0), *histKNN(0), *histHm(0);
223     TH1F *histFi(0), *histFiG(0), *histFiB(0), *histLD(0), *histNn(0),*histNnbfgs(0),*histNnbnn(0);
224     TH1F *histNnC(0), *histNnT(0), *histBdt(0), *histBdtG(0), *histBdtD(0), *histRf(0), *histSVMG(0);
225     TH1F *histSVMP(0), *histSVML(0), *histFDAMT(0), *histFDAGA(0), *histCat(0), *histPBdt(0);
226    
227     if (Use["Likelihood"]) histLk = new TH1F( "MVA_Likelihood", "MVA_Likelihood", nbin, -1, 1 );
228     if (Use["LikelihoodD"]) histLkD = new TH1F( "MVA_LikelihoodD", "MVA_LikelihoodD", nbin, -1, 0.9999 );
229     if (Use["LikelihoodPCA"]) histLkPCA = new TH1F( "MVA_LikelihoodPCA", "MVA_LikelihoodPCA", nbin, -1, 1 );
230     if (Use["LikelihoodKDE"]) histLkKDE = new TH1F( "MVA_LikelihoodKDE", "MVA_LikelihoodKDE", nbin, -0.00001, 0.99999 );
231     if (Use["LikelihoodMIX"]) histLkMIX = new TH1F( "MVA_LikelihoodMIX", "MVA_LikelihoodMIX", nbin, 0, 1 );
232     if (Use["PDERS"]) histPD = new TH1F( "MVA_PDERS", "MVA_PDERS", nbin, 0, 1 );
233     if (Use["PDERSD"]) histPDD = new TH1F( "MVA_PDERSD", "MVA_PDERSD", nbin, 0, 1 );
234     if (Use["PDERSPCA"]) histPDPCA = new TH1F( "MVA_PDERSPCA", "MVA_PDERSPCA", nbin, 0, 1 );
235     if (Use["KNN"]) histKNN = new TH1F( "MVA_KNN", "MVA_KNN", nbin, 0, 1 );
236     if (Use["HMatrix"]) histHm = new TH1F( "MVA_HMatrix", "MVA_HMatrix", nbin, -0.95, 1.55 );
237     if (Use["Fisher"]) histFi = new TH1F( "MVA_Fisher", "MVA_Fisher", nbin, -4, 4 );
238     if (Use["FisherG"]) histFiG = new TH1F( "MVA_FisherG", "MVA_FisherG", nbin, -1, 1 );
239     if (Use["BoostedFisher"]) histFiB = new TH1F( "MVA_BoostedFisher", "MVA_BoostedFisher", nbin, -2, 2 );
240     if (Use["LD"]) histLD = new TH1F( "MVA_LD", "MVA_LD", nbin, -2, 2 );
241     if (Use["MLP"]) histNn = new TH1F( "MVA_MLP", "MVA_MLP", nbin, -1.25, 1.5 );
242     if (Use["MLPBFGS"]) histNnbfgs = new TH1F( "MVA_MLPBFGS", "MVA_MLPBFGS", nbin, -1.25, 1.5 );
243     if (Use["MLPBNN"]) histNnbnn = new TH1F( "MVA_MLPBNN", "MVA_MLPBNN", nbin, -1.25, 1.5 );
244     if (Use["CFMlpANN"]) histNnC = new TH1F( "MVA_CFMlpANN", "MVA_CFMlpANN", nbin, 0, 1 );
245     if (Use["TMlpANN"]) histNnT = new TH1F( "MVA_TMlpANN", "MVA_TMlpANN", nbin, -1.3, 1.3 );
246     if (Use["BDT"]) histBdt = new TH1F( "MVA_BDT", "MVA_BDT", nbin, -0.8, 0.8 );
247     if (Use["BDTD"]) histBdtD = new TH1F( "MVA_BDTD", "MVA_BDTD", nbin, -0.8, 0.8 );
248     if (Use["BDTG"]) histBdtG = new TH1F( "MVA_BDTG", "MVA_BDTG", nbin, -1.0, 1.0 );
249     if (Use["RuleFit"]) histRf = new TH1F( "MVA_RuleFit", "MVA_RuleFit", nbin, -2.0, 2.0 );
250     if (Use["SVM_Gauss"]) histSVMG = new TH1F( "MVA_SVM_Gauss", "MVA_SVM_Gauss", nbin, 0.0, 1.0 );
251     if (Use["SVM_Poly"]) histSVMP = new TH1F( "MVA_SVM_Poly", "MVA_SVM_Poly", nbin, 0.0, 1.0 );
252     if (Use["SVM_Lin"]) histSVML = new TH1F( "MVA_SVM_Lin", "MVA_SVM_Lin", nbin, 0.0, 1.0 );
253     if (Use["FDA_MT"]) histFDAMT = new TH1F( "MVA_FDA_MT", "MVA_FDA_MT", nbin, -2.0, 3.0 );
254     if (Use["FDA_GA"]) histFDAGA = new TH1F( "MVA_FDA_GA", "MVA_FDA_GA", nbin, -2.0, 3.0 );
255     if (Use["Category"]) histCat = new TH1F( "MVA_Category", "MVA_Category", nbin, -2., 2. );
256     if (Use["Plugin"]) histPBdt = new TH1F( "MVA_PBDT", "MVA_BDT", nbin, -0.8, 0.8 );
257    
258     // PDEFoam also returns per-event error, fill in histogram, and also fill significance
259     if (Use["PDEFoam"]) {
260     histPDEFoam = new TH1F( "MVA_PDEFoam", "MVA_PDEFoam", nbin, 0, 1 );
261     histPDEFoamErr = new TH1F( "MVA_PDEFoamErr", "MVA_PDEFoam error", nbin, 0, 1 );
262     histPDEFoamSig = new TH1F( "MVA_PDEFoamSig", "MVA_PDEFoam significance", nbin, 0, 10 );
263     }
264    
265     // Book example histogram for probability (the other methods are done similarly)
266     TH1F *probHistFi(0), *rarityHistFi(0);
267     if (Use["Fisher"]) {
268     probHistFi = new TH1F( "MVA_Fisher_Proba", "MVA_Fisher_Proba", nbin, 0, 1 );
269     rarityHistFi = new TH1F( "MVA_Fisher_Rarity", "MVA_Fisher_Rarity", nbin, 0, 1 );
270     }
271    
272     // Prepare input tree (this must be replaced by your data source)
273     // in this example, there is a toy tree with signal and one with background events
274     // we'll later on use only the "signal" events for the test in this example.
275     //
276    
277 benhoob 1.3 char* prefix = "babies/v3";
278 benhoob 1.1 TChain *ch = new TChain("Events");
279    
280     if( strcmp( samples.at(i) , "DY" ) == 0 ){
281 benhoob 1.3 ch -> Add( Form("%s/DYToMuMuM20_PU_testFinal_baby.root",prefix) );
282     ch -> Add( Form("%s/DYToMuMuM10To20_PU_testFinal_baby.root",prefix) );
283     ch -> Add( Form("%s/DYToEEM20_PU_testFinal_baby.root",prefix) );
284     ch -> Add( Form("%s/DYToEEM10To20_PU_testFinal_baby.root",prefix) );
285     ch -> Add( Form("%s/DYToTauTauM20_PU_testFinal_baby.root",prefix) );
286     ch -> Add( Form("%s/DYToTauTauM10To20_PU_testFinal_baby.root",prefix) );
287 benhoob 1.1 }
288     else if( strcmp( samples.at(i) , "Higgs130" ) == 0 ){
289 benhoob 1.3 ch -> Add( Form("%s/HToWWTo2L2NuM130_PU_testFinal_baby.root",prefix) );
290     ch -> Add( Form("%s/HToWWToLNuTauNuM130_PU_testFinal_baby.root",prefix) );
291     ch -> Add( Form("%s/HToWWTo2Tau2NuM130_PU_testFinal_baby.root",prefix) );
292 benhoob 1.1 }
293     else if( strcmp( samples.at(i) , "Higgs160" ) == 0 ){
294 benhoob 1.3 ch -> Add( Form("%s/HToWWTo2L2NuM160_PU_testFinal_baby.root",prefix) );
295     ch -> Add( Form("%s/HToWWToLNuTauNuM160_PU_testFinal_baby.root",prefix) );
296     ch -> Add( Form("%s/HToWWTo2Tau2NuM160_PU_testFinal_baby.root",prefix) );
297 benhoob 1.1 }
298     else if( strcmp( samples.at(i) , "Higgs200" ) == 0 ){
299 benhoob 1.3 ch -> Add( Form("%s/HToWWTo2L2NuM200_PU_testFinal_baby.root",prefix) );
300     ch -> Add( Form("%s/HToWWToLNuTauNuM200_PU_testFinal_baby.root",prefix) );
301     ch -> Add( Form("%s/HToWWTo2Tau2NuM200_PU_testFinal_baby.root",prefix) );
302 benhoob 1.1 }
303     else{
304 benhoob 1.3 ch -> Add( Form("%s/%s_PU_testFinal_baby.root",prefix,samples.at(i)) );
305 benhoob 1.1 }
306    
307    
308    
309     // --- Event loop
310    
311     // Prepare the event tree
312     // - here the variable names have to corresponds to your tree
313     // - you can use the same variables as above which is slightly faster,
314     // but of course you can use different ones and copy the values inside the event loop
315     //
316    
317     TTree *theTree = (TTree*) ch;
318    
319     std::cout << "--- Using input files: -------------------" << std::endl;
320    
321     TObjArray *listOfFiles = ch->GetListOfFiles();
322     TIter fileIter(listOfFiles);
323     TChainElement* currentFile = 0;
324    
325     while((currentFile = (TChainElement*)fileIter.Next())) {
326     std::cout << currentFile->GetTitle() << std::endl;
327     }
328    
329     Float_t lephard_pt_;
330     Float_t lepsoft_pt_;
331     Float_t dil_dphi_;
332     Float_t dil_mass_;
333     Int_t event_type_;
334     Float_t met_projpt_;
335     Int_t jets_num_;
336     Int_t extralep_num_;
337     Int_t lowptbtags_num_;
338     Int_t softmu_num_;
339     Float_t event_scale1fb_;
340    
341     theTree->SetBranchAddress( "lephard_pt_" , &lephard_pt_ );
342     theTree->SetBranchAddress( "lepsoft_pt_" , &lepsoft_pt_ );
343     theTree->SetBranchAddress( "dil_dphi_" , &dil_dphi_ );
344     theTree->SetBranchAddress( "dil_mass_" , &dil_mass_ );
345     theTree->SetBranchAddress( "event_type_" , &event_type_ );
346     theTree->SetBranchAddress( "met_projpt_" , &met_projpt_ );
347     theTree->SetBranchAddress( "jets_num_" , &jets_num_ );
348     theTree->SetBranchAddress( "extralep_num_" , &extralep_num_ );
349     theTree->SetBranchAddress( "lowptbtags_num_" , &lowptbtags_num_ );
350     theTree->SetBranchAddress( "softmu_num_" , &softmu_num_ );
351     theTree->SetBranchAddress( "event_scale1fb_" , &event_scale1fb_ );
352    
353     // Efficiency calculator for cut method
354     Int_t nSelCutsGA = 0;
355     Double_t effS = 0.7;
356    
357     std::vector<Float_t> vecVar(4); // vector for EvaluateMVA tests
358    
359     std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
360     TStopwatch sw;
361     sw.Start();
362    
363     int npass = 0;
364    
365     for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
366    
367     if (ievt%1000 == 0) std::cout << "--- ... Processing event: " << ievt << std::endl;
368    
369     theTree->GetEntry(ievt);
370    
371 benhoob 1.3 //-------------------------------------------------------
372     // event selection
373     //-------------------------------------------------------
374    
375     //em
376     if( event_type_ == 1 || event_type == 2 ){
377     if( met_projpt_ < 20. ) continue;
378     }
379     //ee/mm
380     if( event_type_ == 0 || event_type == 3 ){
381     if( met_projpt_ < 35. ) continue;
382     }
383    
384 benhoob 1.1 if( lephard_pt_ < 20. ) continue;
385 benhoob 1.2 //if( lepsoft_pt_ < 10. ) continue;
386     if( lepsoft_pt_ < 20. ) continue;
387 benhoob 1.1 if( jets_num_ > 0 ) continue;
388     if( extralep_num_ > 0 ) continue;
389     if( lowptbtags_num_ > 0 ) continue;
390     if( softmu_num_ > 0 ) continue;
391     if( dil_mass_ < 12. ) continue;
392    
393     float weight = event_scale1fb_ * 0.0355;
394    
395     lephard_pt = lephard_pt_;
396     lepsoft_pt = lepsoft_pt_;
397     dil_mass = dil_mass_;
398     dil_dphi = dil_dphi_;
399    
400     npass++;
401     // var1 = userVar1 + userVar2;
402     // var2 = userVar1 - userVar2;
403    
404     // --- Return the MVA outputs and fill into histograms
405    
406     if (Use["CutsGA"]) {
407     // Cuts is a special case: give the desired signal efficienciy
408     Bool_t passed = reader->EvaluateMVA( "CutsGA method", effS );
409     if (passed) nSelCutsGA++;
410     }
411    
412     if (Use["Likelihood" ]) histLk ->Fill( reader->EvaluateMVA( "Likelihood method" ) , weight);
413     if (Use["LikelihoodD" ]) histLkD ->Fill( reader->EvaluateMVA( "LikelihoodD method" ) , weight);
414     if (Use["LikelihoodPCA"]) histLkPCA ->Fill( reader->EvaluateMVA( "LikelihoodPCA method" ) , weight);
415     if (Use["LikelihoodKDE"]) histLkKDE ->Fill( reader->EvaluateMVA( "LikelihoodKDE method" ) , weight);
416     if (Use["LikelihoodMIX"]) histLkMIX ->Fill( reader->EvaluateMVA( "LikelihoodMIX method" ) , weight);
417     if (Use["PDERS" ]) histPD ->Fill( reader->EvaluateMVA( "PDERS method" ) , weight);
418     if (Use["PDERSD" ]) histPDD ->Fill( reader->EvaluateMVA( "PDERSD method" ) , weight);
419     if (Use["PDERSPCA" ]) histPDPCA ->Fill( reader->EvaluateMVA( "PDERSPCA method" ) , weight);
420     if (Use["KNN" ]) histKNN ->Fill( reader->EvaluateMVA( "KNN method" ) , weight);
421     if (Use["HMatrix" ]) histHm ->Fill( reader->EvaluateMVA( "HMatrix method" ) , weight);
422     if (Use["Fisher" ]) histFi ->Fill( reader->EvaluateMVA( "Fisher method" ) , weight);
423     if (Use["FisherG" ]) histFiG ->Fill( reader->EvaluateMVA( "FisherG method" ) , weight);
424     if (Use["BoostedFisher"]) histFiB ->Fill( reader->EvaluateMVA( "BoostedFisher method" ) , weight);
425     if (Use["LD" ]) histLD ->Fill( reader->EvaluateMVA( "LD method" ) , weight);
426     if (Use["MLP" ]) histNn ->Fill( reader->EvaluateMVA( "MLP method" ) , weight);
427     if (Use["MLPBFGS" ]) histNnbfgs ->Fill( reader->EvaluateMVA( "MLPBFGS method" ) , weight);
428     if (Use["MLPBNN" ]) histNnbnn ->Fill( reader->EvaluateMVA( "MLPBNN method" ) , weight);
429     if (Use["CFMlpANN" ]) histNnC ->Fill( reader->EvaluateMVA( "CFMlpANN method" ) , weight);
430     if (Use["TMlpANN" ]) histNnT ->Fill( reader->EvaluateMVA( "TMlpANN method" ) , weight);
431     if (Use["BDT" ]) histBdt ->Fill( reader->EvaluateMVA( "BDT method" ) , weight);
432     if (Use["BDTD" ]) histBdtD ->Fill( reader->EvaluateMVA( "BDTD method" ) , weight);
433     if (Use["BDTG" ]) histBdtG ->Fill( reader->EvaluateMVA( "BDTG method" ) , weight);
434     if (Use["RuleFit" ]) histRf ->Fill( reader->EvaluateMVA( "RuleFit method" ) , weight);
435     if (Use["SVM_Gauss" ]) histSVMG ->Fill( reader->EvaluateMVA( "SVM_Gauss method" ) , weight);
436     if (Use["SVM_Poly" ]) histSVMP ->Fill( reader->EvaluateMVA( "SVM_Poly method" ) , weight);
437     if (Use["SVM_Lin" ]) histSVML ->Fill( reader->EvaluateMVA( "SVM_Lin method" ) , weight);
438     if (Use["FDA_MT" ]) histFDAMT ->Fill( reader->EvaluateMVA( "FDA_MT method" ) , weight);
439     if (Use["FDA_GA" ]) histFDAGA ->Fill( reader->EvaluateMVA( "FDA_GA method" ) , weight);
440     if (Use["Category" ]) histCat ->Fill( reader->EvaluateMVA( "Category method" ) , weight);
441     if (Use["Plugin" ]) histPBdt ->Fill( reader->EvaluateMVA( "P_BDT method" ) , weight);
442    
443     // Retrieve also per-event error
444     if (Use["PDEFoam"]) {
445     Double_t val = reader->EvaluateMVA( "PDEFoam method" );
446     Double_t err = reader->GetMVAError();
447     histPDEFoam ->Fill( val );
448     histPDEFoamErr->Fill( err );
449     if (err>1.e-50) histPDEFoamSig->Fill( val/err , weight);
450     }
451    
452     // Retrieve probability instead of MVA output
453     if (Use["Fisher"]) {
454     probHistFi ->Fill( reader->GetProba ( "Fisher method" ) , weight);
455     rarityHistFi->Fill( reader->GetRarity( "Fisher method" ) , weight);
456     }
457     }
458    
459     std::cout << npass << " events passing selection" << std::endl;
460    
461     // Get elapsed time
462     sw.Stop();
463     std::cout << "--- End of event loop: "; sw.Print();
464    
465     // Get efficiency for cuts classifier
466     if (Use["CutsGA"]) std::cout << "--- Efficiency for CutsGA method: " << double(nSelCutsGA)/theTree->GetEntries()
467     << " (for a required signal efficiency of " << effS << ")" << std::endl;
468    
469     if (Use["CutsGA"]) {
470    
471     // test: retrieve cuts for particular signal efficiency
472     // CINT ignores dynamic_casts so we have to use a cuts-secific Reader function to acces the pointer
473     TMVA::MethodCuts* mcuts = reader->FindCutsMVA( "CutsGA method" ) ;
474    
475     if (mcuts) {
476     std::vector<Double_t> cutsMin;
477     std::vector<Double_t> cutsMax;
478     mcuts->GetCuts( 0.7, cutsMin, cutsMax );
479     std::cout << "--- -------------------------------------------------------------" << std::endl;
480     std::cout << "--- Retrieve cut values for signal efficiency of 0.7 from Reader" << std::endl;
481     for (UInt_t ivar=0; ivar<cutsMin.size(); ivar++) {
482     std::cout << "... Cut: "
483     << cutsMin[ivar]
484     << " < \""
485     << mcuts->GetInputVar(ivar)
486     << "\" <= "
487     << cutsMax[ivar] << std::endl;
488     }
489     std::cout << "--- -------------------------------------------------------------" << std::endl;
490     }
491     }
492    
493     // --- Write histograms
494     cout << "dir " << dir << endl;
495     char* mydir = outdir;
496     TFile *target = new TFile( Form("%s/%s.root",mydir,samples.at(i) ) ,"RECREATE" );
497     cout << "Writing to file " << Form("%s/%s.root",mydir,samples.at(i) ) << endl;
498    
499     if (Use["Likelihood" ]) histLk ->Write();
500     if (Use["LikelihoodD" ]) histLkD ->Write();
501     if (Use["LikelihoodPCA"]) histLkPCA ->Write();
502     if (Use["LikelihoodKDE"]) histLkKDE ->Write();
503     if (Use["LikelihoodMIX"]) histLkMIX ->Write();
504     if (Use["PDERS" ]) histPD ->Write();
505     if (Use["PDERSD" ]) histPDD ->Write();
506     if (Use["PDERSPCA" ]) histPDPCA ->Write();
507     if (Use["KNN" ]) histKNN ->Write();
508     if (Use["HMatrix" ]) histHm ->Write();
509     if (Use["Fisher" ]) histFi ->Write();
510     if (Use["FisherG" ]) histFiG ->Write();
511     if (Use["BoostedFisher"]) histFiB ->Write();
512     if (Use["LD" ]) histLD ->Write();
513     if (Use["MLP" ]) histNn ->Write();
514     if (Use["MLPBFGS" ]) histNnbfgs ->Write();
515     if (Use["MLPBNN" ]) histNnbnn ->Write();
516     if (Use["CFMlpANN" ]) histNnC ->Write();
517     if (Use["TMlpANN" ]) histNnT ->Write();
518     if (Use["BDT" ]) histBdt ->Write();
519     if (Use["BDTD" ]) histBdtD ->Write();
520     if (Use["BDTG" ]) histBdtG ->Write();
521     if (Use["RuleFit" ]) histRf ->Write();
522     if (Use["SVM_Gauss" ]) histSVMG ->Write();
523     if (Use["SVM_Poly" ]) histSVMP ->Write();
524     if (Use["SVM_Lin" ]) histSVML ->Write();
525     if (Use["FDA_MT" ]) histFDAMT ->Write();
526     if (Use["FDA_GA" ]) histFDAGA ->Write();
527     if (Use["Category" ]) histCat ->Write();
528     if (Use["Plugin" ]) histPBdt ->Write();
529    
530     // Write also error and significance histos
531     if (Use["PDEFoam"]) { histPDEFoam->Write(); histPDEFoamErr->Write(); histPDEFoamSig->Write(); }
532    
533     // Write also probability hists
534     if (Use["Fisher"]) { if (probHistFi != 0) probHistFi->Write(); if (rarityHistFi != 0) rarityHistFi->Write(); }
535     target->Close();
536    
537     delete reader;
538    
539     std::cout << "==> TMVAClassificationApplication is done with sample " << samples.at(i) << endl << std::endl;
540     }
541     }