benhoob/HWW/rulevisCorr.C

#include "tmvaglob.C"

// This macro plots the distributions of the different input variables overlaid on
// the sum of importance per bin.
// The scale goes from violett (no importance) to red (high importance).
// Areas where many important rules are active, will thus be very red.
//
// input: - Input file (result from TMVA),
//        - normal/decorrelated/PCA
//        - use of TMVA plotting TStyle
void rulevisCorr( TString fin = "TMVA.root", TMVAGlob::TypeOfPlot type = TMVAGlob::kNorm, bool useTMVAStyle=kTRUE )
{

   // set style and remove existing canvas'
   //   TMVAGlob::Initialize( useTMVAStyle );
   
   // checks if file with name "fin" is already open, and if not opens one
   TFile *file = TMVAGlob::OpenFile( fin );

   // get top dir containing all hists of the variables
//    TDirectory* vardir = (TDirectory*)file->Get( "InputVariables_Id" );
//    TDirectory* vardir = TMVAGlob::GetInputVariablesDir( type );
//    if (vardir==0) return;
   
//    TDirectory* corrdir = TMVAGlob::GetCorrelationPlotsDir( type, vardir );
//    if (corrdir==0) return;

   // get all titles of the method rulefit
   TList titles;
   UInt_t ninst = TMVAGlob::GetListOfTitles("Method_RuleFit",titles);
   if (ninst==0) return;

   TDirectory* dir = (TDirectory*)file->Get( "Method_RuleFit" );

   // loop over rulefit methods 
   TIter next(dir->GetListOfKeys());
   TKey *key(0);   
   while ((key = (TKey*)next())) {

      if (!gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory")) continue;

      TDirectory* rfdir   = (TDirectory*)key->ReadObj();      
      TDirectory* vardir  = rfdir;
      TDirectory* corrdir = rfdir;
      
      // loop over all titles
      TIter keyIter(&titles);
      TDirectory *rfdir;
      TKey *rfkey;
      //      while ((rfkey = TMVAGlob::NextKey(keyIter,"TDirectory"))) {
      //         rfdir = (TDirectory *)rfkey->ReadObj();
         rulevisCorr( rfdir, vardir, corrdir, type );
         //      }
   }
}

void rulevisCorr( TDirectory *rfdir, TDirectory *vardir, TDirectory *corrdir, TMVAGlob::TypeOfPlot type) {
   //
   if (rfdir==0)   return;
   if (vardir==0)  return;
   if (corrdir==0) return;
   //
   const TString rfName    = rfdir->GetName();
   const TString maintitle = rfName + " : Rule Importance, 2D";
   const TString rfNameOpt = "_RF2D_";
   const TString outfname[TMVAGlob::kNumOfMethods] = { "rulevisCorr",
                                                       "rulevisCorr_decorr",
                                                       "rulevisCorr_pca",
                                                       "rulevisCorr_gaussdecorr" };
   const TString outputName = outfname[type]+"_"+rfdir->GetName();
   //
   TIter rfnext(rfdir->GetListOfKeys());
   TKey *rfkey;
   Double_t rfmax = -1;
   Double_t rfmin = -1;
   Bool_t allEmpty=kTRUE;
   Bool_t first=kTRUE;
   TH2F *hrf;
   while ((rfkey = TMVAGlob::NextKey(rfnext,"TH2F"))) {
      hrf = (TH2F *)rfkey->ReadObj();
      TString hname= hrf->GetName();
      if (hname.Contains(rfNameOpt)){ // found a new RF2D plot
         Double_t valmin = hrf->GetMinimum();
         Double_t valmax = hrf->GetMaximum();
         if (first) {
            rfmax = valmax;
            rfmin = valmin;
         } else {
            if (valmax>rfmax) rfmax = valmax;
            if (valmin<rfmin) rfmin = valmin;
         }
         if (hrf->GetEntries()>0) allEmpty=kFALSE;
         first=kFALSE;
      }
   }
   if (first) {
      cout << "ERROR: no RF2D plots found..." << endl;
      return;
   }
   Double_t minrange   = rfmin;
   Double_t maxrange   = rfmax;
   Double_t targetrange = maxrange - minrange;

   const Int_t nContours = 100;
   Double_t contourLevels[nContours];
   Double_t dcl = targetrange/Double_t(nContours-1);
   //
   for (Int_t i=0; i<nContours; i++) {
      contourLevels[i] = rfmin+dcl*Double_t(i);
   }

   ///////////////////////////
   vardir->cd();
 
   // how many plots are in the directory?
   Int_t noVars = ((vardir->GetListOfKeys())->GetEntries()) / 2;
   Int_t noPlots = (noVars*(noVars+1)/2) - noVars;

   // *** CONTINUE HERE *** 
   // define Canvas layout here!
   // default setting
   Int_t xPad;  // no of plots in x
   Int_t yPad;  // no of plots in y
   Int_t width; // size of canvas
   Int_t height;
   switch (noPlots) {
   case 1:
      xPad = 1; yPad = 1; width = 500; height = 0.7*width; break;
   case 2:
      xPad = 2; yPad = 1; width = 600; height = 0.7*width; break;
   case 3:
      xPad = 3; yPad = 1; width = 900; height = 0.4*width; break;
   case 4:
      xPad = 2; yPad = 2; width = 600; height = width; break;
   default:
      xPad = 3; yPad = 2; width = 800; height = 0.7*width; break;
   }
   Int_t noPad = xPad * yPad ;   

   // this defines how many canvases we need
   const Int_t noCanvas = 1 + (Int_t)((noPlots - 0.001)/noPad);
   TCanvas **c = new TCanvas*[noCanvas];
   for (Int_t ic=0; ic<noCanvas; ic++) c[ic] = 0;

   // counter variables
   Int_t countCanvas = 0;
   Int_t countPad    = 1;

   // loop over all objects in directory
   TIter next(corrdir->GetListOfKeys());
   TKey *key;
   TH2F *sigCpy=0;
   TH2F *bgdCpy=0;
   Bool_t first = kTRUE;
   //
   while ((key = (TKey*)next())) {

      // make sure, that we only look at histograms
      TClass *cl = gROOT->GetClass(key->GetClassName());
      if (!cl->InheritsFrom("TH2")) continue;
      sig = (TH2F*)key->ReadObj();
      TString hname= sig->GetName();
      // check for all signal histograms
      if (hname.Contains("_sig_")){ // found a new signal plot

         // create new canvas
         if ((c[countCanvas]==NULL) || (countPad>noPad)) {
            char cn[20];
            sprintf( cn, "rulecorr%d_", countCanvas+1 );
            TString cname(cn);
            cname += rfdir->GetName();
            c[countCanvas] = new TCanvas( cname, maintitle,
                                          countCanvas*50+200, countCanvas*20, width, height ); 
            // style
            c[countCanvas]->Divide(xPad,yPad);
            countPad = 1;
         }       
         // save canvas to file
         TPad *cPad = (TPad *)(c[countCanvas]->GetPad(countPad));
         c[countCanvas]->cd(countPad);
         countPad++;

         // find the corredponding background histo
         TString bgname = hname;
         bgname.ReplaceAll("_sig_","_bgd_");
         hkey = corrdir->GetKey(bgname);
         bgd = (TH2F*)hkey->ReadObj();
         if (bgd == NULL) {
            cout << "ERROR!!! couldn't find backgroung histo for" << hname << endl;
            exit;
         }
         const Int_t rebin=6;
         sig->Rebin2D(rebin,rebin);
         bgd->Rebin2D(rebin,rebin);
         //
         TString rfname = hname;
         rfname.ReplaceAll("_sig_",rfNameOpt);
         TKey *hrfkey = rfdir->GetKey(rfname);
         hrf2 = (TH2F*)hrfkey->ReadObj();
         Double_t wmin = hrf2->GetMinimum();
         Double_t wmax = hrf2->GetMaximum();
         Double_t wmean  = (wmax+wmin)/2.0;
         Double_t wrange = wmax-wmin;
         Double_t wscale = (wrange>0.0 ? targetrange/wrange : 1.0);
         //         if (rfmax>0.0)
         //            hrf2->Scale(1.0/rfmax);
         hrf2->SetMinimum(minrange); // make sure it's zero  -> for palette axis
         hrf2->SetMaximum(maxrange); // make sure max is 1.0 -> idem
         hrf2->SetContour(nContours,&contourLevels[0]);

         // this is set but not stored during plot creation in MVA_Factory
         //         TMVAGlob::SetSignalAndBackgroundStyle( sigK, bgd );
         sig->SetFillColor(1);
         sig->SetLineColor(1);

         bgd->SetFillColor(15);
         bgd->SetLineColor(15);

         // chop off "signal" 
         TString title(hrf2->GetTitle());
         title.ReplaceAll("signal","");
         if (first) {
            hrf2->SetTitle( maintitle );
            first=kFALSE;
         } else {
            hrf2->SetTitle( "" );
         }
         TMVAGlob::SetFrameStyle( hrf2, 1.2 );

         // finally plot and overlay       
         hrf2->Draw("colz ah");
         Float_t sc = 1.1;
         if (countPad==2) sc = 1.3;
         sig->SetMaximum( TMath::Max( sig->GetMaximum(), bgd->GetMaximum() )*sc );
         Double_t smax = sig->GetMaximum();

         sig->Scale(1.0/smax);
         sig->SetContour(5);
         sig->Draw("same cont3");
         TMVAGlob::SetFrameStyle( sig, 1.2 );

         bgd->Scale(1.0/smax);
         bgd->SetContour(5);
         bgd->Draw("same cont3");
         TMVAGlob::SetFrameStyle( bgd, 1.2 );
         //         sig->GetXaxis()->SetTitle( title );
         sig->GetYaxis()->SetTitleOffset( 1.30 );
         //         sig->GetYaxis()->SetTitle("Events");

         // redraw axes
         sig->Draw("sameaxis");

         cPad->SetRightMargin(0.13);
         cPad->Update();

         // Draw legend
         if (countPad==2){
            TLegend *legend= new TLegend( cPad->GetLeftMargin(), 
                                          1-cPad->GetTopMargin()-.18, 
                                          cPad->GetLeftMargin()+.4, 
                                          1-cPad->GetTopMargin() );
            legend->AddEntry(sig,"Signal","F");
            legend->AddEntry(bgd,"Background","F");
            legend->Draw("same");
            legend->SetBorderSize(1);
            legend->SetMargin( 0.3 );
            legend->SetFillColor(19);
            legend->SetFillStyle(3001);
         } 

         // save canvas to file
         if (countPad > noPad) {
            c[countCanvas]->Update();
            TString fname = Form( "plots/%s_c%i", outputName.Data(), countCanvas+1 );
            TMVAGlob::imgconv( c[countCanvas], fname );
            //        TMVAGlob::plot_logo(); // don't understand why this doesn't work ... :-(
            countCanvas++;
         }
      }
   }

   if (countPad <= noPad) {
      c[countCanvas]->Update();
      TString fname = Form( "plots/%s_c%i", outfname[type].Data(), countCanvas+1 );
      TMVAGlob::imgconv( c[countCanvas], fname );
   }
}
Revision:	1.1
Committed:	Mon Feb 14 12:39:14 2011 UTC (14 years, 3 months ago) by benhoob
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	Initial commit
#	User	Rev	Content
1	benhoob	1.1	#include "tmvaglob.C"
2
3			// This macro plots the distributions of the different input variables overlaid on
4			// the sum of importance per bin.
5			// The scale goes from violett (no importance) to red (high importance).
6			// Areas where many important rules are active, will thus be very red.
7			//
8			// input: - Input file (result from TMVA),
9			// - normal/decorrelated/PCA
10			// - use of TMVA plotting TStyle
11			void rulevisCorr( TString fin = "TMVA.root", TMVAGlob::TypeOfPlot type = TMVAGlob::kNorm, bool useTMVAStyle=kTRUE )
12			{
13
14			// set style and remove existing canvas'
15			// TMVAGlob::Initialize( useTMVAStyle );
16
17			// checks if file with name "fin" is already open, and if not opens one
18			TFile *file = TMVAGlob::OpenFile( fin );
19
20			// get top dir containing all hists of the variables
21			// TDirectory* vardir = (TDirectory*)file->Get( "InputVariables_Id" );
22			// TDirectory* vardir = TMVAGlob::GetInputVariablesDir( type );
23			// if (vardir==0) return;
24
25			// TDirectory* corrdir = TMVAGlob::GetCorrelationPlotsDir( type, vardir );
26			// if (corrdir==0) return;
27
28			// get all titles of the method rulefit
29			TList titles;
30			UInt_t ninst = TMVAGlob::GetListOfTitles("Method_RuleFit",titles);
31			if (ninst==0) return;
32
33			TDirectory* dir = (TDirectory*)file->Get( "Method_RuleFit" );
34
35			// loop over rulefit methods
36			TIter next(dir->GetListOfKeys());
37			TKey *key(0);
38			while ((key = (TKey*)next())) {
39
40			if (!gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory")) continue;
41
42			TDirectory* rfdir = (TDirectory*)key->ReadObj();
43			TDirectory* vardir = rfdir;
44			TDirectory* corrdir = rfdir;
45
46			// loop over all titles
47			TIter keyIter(&titles);
48			TDirectory *rfdir;
49			TKey *rfkey;
50			// while ((rfkey = TMVAGlob::NextKey(keyIter,"TDirectory"))) {
51			// rfdir = (TDirectory *)rfkey->ReadObj();
52			rulevisCorr( rfdir, vardir, corrdir, type );
53			// }
54			}
55			}
56
57			void rulevisCorr( TDirectory rfdir, TDirectory vardir, TDirectory *corrdir, TMVAGlob::TypeOfPlot type) {
58			//
59			if (rfdir==0) return;
60			if (vardir==0) return;
61			if (corrdir==0) return;
62			//
63			const TString rfName = rfdir->GetName();
64			const TString maintitle = rfName + " : Rule Importance, 2D";
65			const TString rfNameOpt = "_RF2D_";
66			const TString outfname[TMVAGlob::kNumOfMethods] = { "rulevisCorr",
67			"rulevisCorr_decorr",
68			"rulevisCorr_pca",
69			"rulevisCorr_gaussdecorr" };
70			const TString outputName = outfname[type]+"_"+rfdir->GetName();
71			//
72			TIter rfnext(rfdir->GetListOfKeys());
73			TKey *rfkey;
74			Double_t rfmax = -1;
75			Double_t rfmin = -1;
76			Bool_t allEmpty=kTRUE;
77			Bool_t first=kTRUE;
78			TH2F *hrf;
79			while ((rfkey = TMVAGlob::NextKey(rfnext,"TH2F"))) {
80			hrf = (TH2F *)rfkey->ReadObj();
81			TString hname= hrf->GetName();
82			if (hname.Contains(rfNameOpt)){ // found a new RF2D plot
83			Double_t valmin = hrf->GetMinimum();
84			Double_t valmax = hrf->GetMaximum();
85			if (first) {
86			rfmax = valmax;
87			rfmin = valmin;
88			} else {
89			if (valmax>rfmax) rfmax = valmax;
90			if (valmin<rfmin) rfmin = valmin;
91			}
92			if (hrf->GetEntries()>0) allEmpty=kFALSE;
93			first=kFALSE;
94			}
95			}
96			if (first) {
97			cout << "ERROR: no RF2D plots found..." << endl;
98			return;
99			}
100			Double_t minrange = rfmin;
101			Double_t maxrange = rfmax;
102			Double_t targetrange = maxrange - minrange;
103
104			const Int_t nContours = 100;
105			Double_t contourLevels[nContours];
106			Double_t dcl = targetrange/Double_t(nContours-1);
107			//
108			for (Int_t i=0; i<nContours; i++) {
109			contourLevels[i] = rfmin+dcl*Double_t(i);
110			}
111
112			///////////////////////////
113			vardir->cd();
114
115			// how many plots are in the directory?
116			Int_t noVars = ((vardir->GetListOfKeys())->GetEntries()) / 2;
117			Int_t noPlots = (noVars*(noVars+1)/2) - noVars;
118
119			// * CONTINUE HERE *
120			// define Canvas layout here!
121			// default setting
122			Int_t xPad; // no of plots in x
123			Int_t yPad; // no of plots in y
124			Int_t width; // size of canvas
125			Int_t height;
126			switch (noPlots) {
127			case 1:
128			xPad = 1; yPad = 1; width = 500; height = 0.7*width; break;
129			case 2:
130			xPad = 2; yPad = 1; width = 600; height = 0.7*width; break;
131			case 3:
132			xPad = 3; yPad = 1; width = 900; height = 0.4*width; break;
133			case 4:
134			xPad = 2; yPad = 2; width = 600; height = width; break;
135			default:
136			xPad = 3; yPad = 2; width = 800; height = 0.7*width; break;
137			}
138			Int_t noPad = xPad * yPad ;
139
140			// this defines how many canvases we need
141			const Int_t noCanvas = 1 + (Int_t)((noPlots - 0.001)/noPad);
142			TCanvas *c = new TCanvas[noCanvas];
143			for (Int_t ic=0; ic<noCanvas; ic++) c[ic] = 0;
144
145			// counter variables
146			Int_t countCanvas = 0;
147			Int_t countPad = 1;
148
149			// loop over all objects in directory
150			TIter next(corrdir->GetListOfKeys());
151			TKey *key;
152			TH2F *sigCpy=0;
153			TH2F *bgdCpy=0;
154			Bool_t first = kTRUE;
155			//
156			while ((key = (TKey*)next())) {
157
158			// make sure, that we only look at histograms
159			TClass *cl = gROOT->GetClass(key->GetClassName());
160			if (!cl->InheritsFrom("TH2")) continue;
161			sig = (TH2F*)key->ReadObj();
162			TString hname= sig->GetName();
163			// check for all signal histograms
164			if (hname.Contains("_sig_")){ // found a new signal plot
165
166			// create new canvas
167			if ((c[countCanvas]==NULL) \|\| (countPad>noPad)) {
168			char cn[20];
169			sprintf( cn, "rulecorr%d_", countCanvas+1 );
170			TString cname(cn);
171			cname += rfdir->GetName();
172			c[countCanvas] = new TCanvas( cname, maintitle,
173			countCanvas50+200, countCanvas20, width, height );
174			// style
175			c[countCanvas]->Divide(xPad,yPad);
176			countPad = 1;
177			}
178			// save canvas to file
179			TPad cPad = (TPad )(c[countCanvas]->GetPad(countPad));
180			c[countCanvas]->cd(countPad);
181			countPad++;
182
183			// find the corredponding background histo
184			TString bgname = hname;
185			bgname.ReplaceAll("_sig_","_bgd_");
186			hkey = corrdir->GetKey(bgname);
187			bgd = (TH2F*)hkey->ReadObj();
188			if (bgd == NULL) {
189			cout << "ERROR!!! couldn't find backgroung histo for" << hname << endl;
190			exit;
191			}
192			const Int_t rebin=6;
193			sig->Rebin2D(rebin,rebin);
194			bgd->Rebin2D(rebin,rebin);
195			//
196			TString rfname = hname;
197			rfname.ReplaceAll("_sig_",rfNameOpt);
198			TKey *hrfkey = rfdir->GetKey(rfname);
199			hrf2 = (TH2F*)hrfkey->ReadObj();
200			Double_t wmin = hrf2->GetMinimum();
201			Double_t wmax = hrf2->GetMaximum();
202			Double_t wmean = (wmax+wmin)/2.0;
203			Double_t wrange = wmax-wmin;
204			Double_t wscale = (wrange>0.0 ? targetrange/wrange : 1.0);
205			// if (rfmax>0.0)
206			// hrf2->Scale(1.0/rfmax);
207			hrf2->SetMinimum(minrange); // make sure it's zero -> for palette axis
208			hrf2->SetMaximum(maxrange); // make sure max is 1.0 -> idem
209			hrf2->SetContour(nContours,&contourLevels[0]);
210
211			// this is set but not stored during plot creation in MVA_Factory
212			// TMVAGlob::SetSignalAndBackgroundStyle( sigK, bgd );
213			sig->SetFillColor(1);
214			sig->SetLineColor(1);
215
216			bgd->SetFillColor(15);
217			bgd->SetLineColor(15);
218
219			// chop off "signal"
220			TString title(hrf2->GetTitle());
221			title.ReplaceAll("signal","");
222			if (first) {
223			hrf2->SetTitle( maintitle );
224			first=kFALSE;
225			} else {
226			hrf2->SetTitle( "" );
227			}
228			TMVAGlob::SetFrameStyle( hrf2, 1.2 );
229
230			// finally plot and overlay
231			hrf2->Draw("colz ah");
232			Float_t sc = 1.1;
233			if (countPad==2) sc = 1.3;
234			sig->SetMaximum( TMath::Max( sig->GetMaximum(), bgd->GetMaximum() )*sc );
235			Double_t smax = sig->GetMaximum();
236
237			sig->Scale(1.0/smax);
238			sig->SetContour(5);
239			sig->Draw("same cont3");
240			TMVAGlob::SetFrameStyle( sig, 1.2 );
241
242			bgd->Scale(1.0/smax);
243			bgd->SetContour(5);
244			bgd->Draw("same cont3");
245			TMVAGlob::SetFrameStyle( bgd, 1.2 );
246			// sig->GetXaxis()->SetTitle( title );
247			sig->GetYaxis()->SetTitleOffset( 1.30 );
248			// sig->GetYaxis()->SetTitle("Events");
249
250			// redraw axes
251			sig->Draw("sameaxis");
252
253			cPad->SetRightMargin(0.13);
254			cPad->Update();
255
256			// Draw legend
257			if (countPad==2){
258			TLegend *legend= new TLegend( cPad->GetLeftMargin(),
259			1-cPad->GetTopMargin()-.18,
260			cPad->GetLeftMargin()+.4,
261			1-cPad->GetTopMargin() );
262			legend->AddEntry(sig,"Signal","F");
263			legend->AddEntry(bgd,"Background","F");
264			legend->Draw("same");
265			legend->SetBorderSize(1);
266			legend->SetMargin( 0.3 );
267			legend->SetFillColor(19);
268			legend->SetFillStyle(3001);
269			}
270
271			// save canvas to file
272			if (countPad > noPad) {
273			c[countCanvas]->Update();
274			TString fname = Form( "plots/%s_c%i", outputName.Data(), countCanvas+1 );
275			TMVAGlob::imgconv( c[countCanvas], fname );
276			// TMVAGlob::plot_logo(); // don't understand why this doesn't work ... :-(
277			countCanvas++;
278			}
279			}
280			}
281
282			if (countPad <= noPad) {
283			c[countCanvas]->Update();
284			TString fname = Form( "plots/%s_c%i", outfname[type].Data(), countCanvas+1 );
285			TMVAGlob::imgconv( c[countCanvas], fname );
286			}
287			}