Logo ROOT  
Reference Guide
HybridPlot.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2
3/** \class RooStats::HybridPlot
4 \ingroup Roostats
5
6This class provides the plots for the result of a study performed with the
7HybridCalculatorOriginal class.
8
9Authors: D. Piparo, G. Schott - Universitaet Karlsruhe
10
11An example plot is available here:
12http://www-ekp.physik.uni-karlsruhe.de/~schott/roostats/hybridplot_example.png
13*/
14
15#include "assert.h"
16#include <cmath>
17#include <iostream>
18#include <map>
19
20#include "RooStats/HybridPlot.h"
21#include "TStyle.h"
22#include "TF1.h"
23#include "TAxis.h"
24#include "TH1.h"
25#include "TLine.h"
26#include "TLegend.h"
27#include "TFile.h"
28#include "TVirtualPad.h"
29
30#include <algorithm>
31
32/// To build the THtml documentation
33using namespace std;
34
36
37using namespace RooStats;
38
39////////////////////////////////////////////////////////////////////////////////
40/// HybridPlot constructor
41
42HybridPlot::HybridPlot(const char* name,
43 const char* title,
44 const std::vector<double> & sb_vals,
45 const std::vector<double> & b_vals,
46 double testStat_data,
47 int n_bins,
48 bool verbosity):
49 TNamed(name,title),
50 fSb_histo(NULL),
51 fSb_histo_shaded(NULL),
52 fB_histo(NULL),
53 fB_histo_shaded(NULL),
54 fData_testStat_line(0),
55 fLegend(0),
56 fPad(0),
57 fVerbose(verbosity)
58{
59 int nToysSB = sb_vals.size();
60 int nToysB = sb_vals.size();
61 assert (nToysSB >0);
62 assert (nToysB >0);
63
64 // Get the max and the min of the plots
65 double min = *std::min_element(sb_vals.begin(), sb_vals.end());
66 double max = *std::max_element(sb_vals.begin(), sb_vals.end());
67
68 double min_b = *std::min_element(b_vals.begin(), b_vals.end());
69 double max_b = *std::max_element(b_vals.begin(), b_vals.end());
70
71
72 if ( min_b < min) min = min_b;
73 if ( max_b > max) max = max_b;
74
75 if (testStat_data<min) min = testStat_data;
76 if (testStat_data>max) max = testStat_data;
77
78 min *= 1.1;
79 max *= 1.1;
80
81 // Build the histos
82
83 fSb_histo = new TH1F ("SB_model",title,n_bins,min,max);
86 fSb_histo->GetXaxis()->SetTitle("test statistics");
88
89 fB_histo = new TH1F ("B_model",title,n_bins,min,max);
92 fB_histo->GetXaxis()->SetTitle("test statistics");
94
95 for (int i=0;i<nToysSB;++i) fSb_histo->Fill(sb_vals[i]);
96 for (int i=0;i<nToysB;++i) fB_histo->Fill(b_vals[i]);
97
98 double histos_max_y = fSb_histo->GetMaximum();
99 double line_hight = histos_max_y/nToysSB;
100 if (histos_max_y<fB_histo->GetMaximum()) histos_max_y = fB_histo->GetMaximum()/nToysB;
101
102 // Build the line of the measured -2lnQ
103 fData_testStat_line = new TLine(testStat_data,0,testStat_data,line_hight);
106
107 // The legend
108 double golden_section = (std::sqrt(5.)-1)/2;
109
110 fLegend = new TLegend(0.75,0.95-0.2*golden_section,0.95,0.95);
111 TString title_leg="test statistics distributions ";
112 title_leg+=sb_vals.size();
113 title_leg+=" toys";
114 fLegend->SetName(title_leg.Data());
115 fLegend->AddEntry(fSb_histo,"SB toy datasets");
116 fLegend->AddEntry(fB_histo,"B toy datasets");
117 fLegend->AddEntry((TLine*)fData_testStat_line,"test statistics on data","L");
119}
120
121////////////////////////////////////////////////////////////////////////////////
122/// destructor
123
125{
126
127 if (fSb_histo) delete fSb_histo;
128 if (fB_histo) delete fB_histo;
132 if (fLegend) delete fLegend;
133}
134
135////////////////////////////////////////////////////////////////////////////////
136/// draw the S+B and B histogram in the current canvas
137
138void HybridPlot::Draw(const char* )
139{
140 // We don't want the statistics of the histos
141 gStyle->SetOptStat(0);
142
143 // The histos
144
147 fB_histo->DrawNormalized("same");
148 }
149 else{
151 fSb_histo->DrawNormalized("same");
152 }
153
154 // Shaded
155 fB_histo_shaded = (TH1F*)fB_histo->Clone("b_shaded");
158
159 fSb_histo_shaded = (TH1F*)fSb_histo->Clone("sb_shaded");
162
163 // Empty the bins according to the data -2lnQ
164 double data_m2lnq= fData_testStat_line->GetX1();
165 for (int i=1;i<=fSb_histo->GetNbinsX();++i){
166 if (fSb_histo->GetBinCenter(i)<data_m2lnq){
169 }
170 else{
173 }
174 }
175
176 // Draw the shaded histos
177 fB_histo_shaded->Draw("same");
178 fSb_histo_shaded->Draw("same");
179
180 // The line
181 fData_testStat_line->Draw("same");
182
183 // The legend
184 fLegend->Draw("same");
185
186 if (gPad) {
187 gPad->SetName(GetName());
188 gPad->SetTitle(GetTitle() );
189 }
190
191 fPad = gPad;
192
193}
194
195////////////////////////////////////////////////////////////////////////////////
196/// All the objects are written to rootfile
197
198void HybridPlot::DumpToFile (const char* RootFileName, const char* options)
199{
200
201 TFile ofile(RootFileName,options);
202 ofile.cd();
203
204 // The histos
205 fSb_histo->Write();
206 fB_histo->Write();
207
208 // The shaded histos
209 if (fB_histo_shaded!=NULL && fSb_histo_shaded!=NULL){
212 }
213
214 // The line
215 fData_testStat_line->Write("Measured test statistics line tag");
216
217 // The legend
218 fLegend->Write();
219
220 ofile.Close();
221
222}
223
224////////////////////////////////////////////////////////////////////////////////
225
226void HybridPlot::DumpToImage(const char * filename) {
227 if (!fPad) {
228 Error("HybridPlot","Hybrid plot has not yet been drawn ");
229 return;
230 }
231 fPad->Print(filename);
232}
233
234////////////////////////////////////////////////////////////////////////////////
235/// Perform 2 times a gaussian fit to fetch the center of the histo.
236/// To get the second fit range get an interval that tries to keep into account
237/// the skewness of the distribution.
238
239double HybridPlot::GetHistoCenter(TH1* histo_orig, double n_rms, bool display_result){
240 // Get the center of the histo
241
242 TString optfit = "Q0";
243 if (display_result) optfit = "Q";
244
245 TH1F* histo = (TH1F*)histo_orig->Clone();
246
247 // get the histo x extremes
248 double x_min = histo->GetXaxis()->GetXmin();
249 double x_max = histo->GetXaxis()->GetXmax();
250
251 // First fit!
252
253 TF1* gaus = new TF1("mygaus", "gaus", x_min, x_max);
254
255 gaus->SetParameter("Constant",histo->GetEntries());
256 gaus->SetParameter("Mean",histo->GetMean());
257 gaus->SetParameter("Sigma",histo->GetRMS());
258
259 histo->Fit(gaus,optfit);
260
261 // Second fit!
262 double sigma = gaus->GetParameter("Sigma");
263 double mean = gaus->GetParameter("Mean");
264
265 delete gaus;
266
267 std::cout << "Center is 1st pass = " << mean << std::endl;
268
269 double skewness = histo->GetSkewness();
270
271 x_min = mean - n_rms*sigma - sigma*skewness/2;
272 x_max = mean + n_rms*sigma - sigma*skewness/2;;
273
274 TF1* gaus2 = new TF1("mygaus2", "gaus", x_min, x_max);
275 gaus2->SetParameter("Mean",mean);
276
277 // second fit : likelihood fit
278 optfit += "L";
279 histo->Fit(gaus2,optfit,"", x_min, x_max);
280
281
282 double center = gaus2->GetParameter("Mean");
283
284 if (display_result) {
285 histo->Draw();
286 gaus2->Draw("same");
287 }
288 else {
289 delete histo;
290 }
291 delete gaus2;
292
293 return center;
294
295
296}
297
298////////////////////////////////////////////////////////////////////////////////
299/// We let an horizontal bar go down and we stop when we have the integral
300/// equal to the desired one.
301
302double* HybridPlot::GetHistoPvals (TH1* histo, double percentage){
303
304 if (percentage>1){
305 std::cerr << "Percentage greater or equal to 1!\n";
306 return NULL;
307 }
308
309 // Get the integral of the histo
310 double* h_integral=histo->GetIntegral();
311
312 // Start the iteration
313 std::map<int,int> extremes_map;
314
315 for (int i=0;i<histo->GetNbinsX();++i){
316 for (int j=0;j<histo->GetNbinsX();++j){
317 double integral = h_integral[j]-h_integral[i];
318 if (integral>percentage){
319 extremes_map[i]=j;
320 break;
321 }
322 }
323 }
324
325 // Now select the couple of extremes which have the lower bin content diff
326 std::map<int,int>::iterator it;
327 int a,b;
328 double left_bin_center(0.),right_bin_center(0.);
329 double diff=10e40;
330 double current_diff;
331 for (it = extremes_map.begin();it != extremes_map.end();++it){
332 a=it->first;
333 b=it->second;
334 current_diff=std::fabs(histo->GetBinContent(a)-histo->GetBinContent(b));
335 if (current_diff<diff){
336 //std::cout << "a=" << a << " b=" << b << std::endl;
337 diff=current_diff;
338 left_bin_center=histo->GetBinCenter(a);
339 right_bin_center=histo->GetBinCenter(b);
340 }
341 }
342
343 double* d = new double[2];
344 d[0]=left_bin_center;
345 d[1]=right_bin_center;
346 return d;
347}
348
349////////////////////////////////////////////////////////////////////////////////
350/// Get the median of an histogram.
351
353
354 //int xbin_median;
355 Double_t* integral = histo->GetIntegral();
356 int median_i = 0;
357 for (int j=0;j<histo->GetNbinsX(); j++)
358 if (integral[j]<0.5)
359 median_i = j;
360
361 double median_x =
362 histo->GetBinCenter(median_i)+
363 (histo->GetBinCenter(median_i+1)-
364 histo->GetBinCenter(median_i))*
365 (0.5-integral[median_i])/(integral[median_i+1]-integral[median_i]);
366
367 return median_x;
368}
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define ClassImp(name)
Definition: Rtypes.h:361
@ kRed
Definition: Rtypes.h:64
@ kBlack
Definition: Rtypes.h:63
@ kBlue
Definition: Rtypes.h:64
char name[80]
Definition: TGX11.cxx:109
double sqrt(double)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
#define gPad
Definition: TVirtualPad.h:287
This class provides the plots for the result of a study performed with the HybridCalculatorOriginal c...
Definition: HybridPlot.h:36
double GetMedian(TH1 *histo)
Get the median of an histogram.
Definition: HybridPlot.cxx:352
TVirtualPad * fPad
Definition: HybridPlot.h:118
void DumpToImage(const char *filename)
Write an image on disk.
Definition: HybridPlot.cxx:226
TLine * fData_testStat_line
Definition: HybridPlot.h:116
void Draw(const char *options="")
Draw on current pad.
Definition: HybridPlot.cxx:138
~HybridPlot()
Destructor.
Definition: HybridPlot.cxx:124
double * GetHistoPvals(TH1 *histo, double percentage)
Get the "effective sigmas" of the histo, call delete [] res to release memory.
Definition: HybridPlot.cxx:302
double GetHistoCenter(TH1 *histo, double n_rms=1, bool display_result=false)
Get the center of the histo.
Definition: HybridPlot.cxx:239
void DumpToFile(const char *RootFileName, const char *options)
All the objects are written to rootfile.
Definition: HybridPlot.cxx:198
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
Double_t GetXmax() const
Definition: TAxis.h:134
Double_t GetXmin() const
Definition: TAxis.h:133
Bool_t cd(const char *path=nullptr) override
Change current directory to "this" directory.
1-Dim function class
Definition: TF1.h:210
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1320
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:628
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:506
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
void Close(Option_t *option="") override
Close a file.
Definition: TFile.cxx:873
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
The TH1 histogram class.
Definition: TH1.h:56
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6345
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8597
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:7086
virtual Double_t GetSkewness(Int_t axis=1) const
Definition: TH1.cxx:7192
virtual TH1 * DrawNormalized(Option_t *option="", Double_t norm=1) const
Draw a normalized copy of this histogram.
Definition: TH1.cxx:3075
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2665
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3808
virtual Double_t GetMaximum(Double_t maxval=FLT_MAX) const
Return maximum value smaller than maxval of bins in the range, unless the value has been overridden b...
Definition: TH1.cxx:8006
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
Double_t GetRMS(Int_t axis=1) const
Definition: TH1.h:311
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8678
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4302
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t * GetIntegral()
Return a pointer to the array of bins integral.
Definition: TH1.cxx:2523
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7426
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:330
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:423
A simple line.
Definition: TLine.h:23
Double_t GetX1() const
Definition: TLine.h:51
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:796
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
virtual void SetName(const char *name="")
Definition: TPave.h:75
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1590
virtual void Print(const char *filename="") const =0
Print function.
const Double_t sigma
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Namespace for the RooStats classes.
Definition: Asimov.h:19
auto * a
Definition: textangle.C:12