ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/benhoob/HWW/Classify_HWW.C
Revision: 1.6
Committed: Wed Feb 23 15:14:14 2011 UTC (14 years, 2 months ago) by benhoob
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.5: +83 -71 lines
Log Message:
Minor updates

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