ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/benhoob/HWW/Classify_HWW.C
Revision: 1.5
Committed: Mon Feb 21 14:30:47 2011 UTC (14 years, 2 months ago) by benhoob
Content type: text/plain
Branch: MAIN
Changes since 1.4: +185 -64 lines
Log Message:
Make code user-friendly

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