Logo ROOT   6.12/07
Reference Guide
HybridPlot.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 
3 /** \class RooStats::HybridPlot
4  \ingroup Roostats
5 
6 This class provides the plots for the result of a study performed with the
7 HybridCalculatorOriginal class.
8 
9 Authors: D. Piparo, G. Schott - Universitaet Karlsruhe
10 
11 An example plot is available here:
12 http://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
33 using namespace std;
34 
36 
37 using namespace RooStats;
38 
39 ////////////////////////////////////////////////////////////////////////////////
40 /// HybridPlot constructor
41 
42 HybridPlot::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");
118  fLegend->SetFillColor(0);
119 }
120 
121 ////////////////////////////////////////////////////////////////////////////////
122 /// destructor
123 
125 {
126 
127  if (fSb_histo) delete fSb_histo;
128  if (fB_histo) delete fB_histo;
130  if (fB_histo_shaded) delete fB_histo_shaded;
132  if (fLegend) delete fLegend;
133 }
134 
135 ////////////////////////////////////////////////////////////////////////////////
136 /// draw the S+B and B histogram in the current canvas
137 
138 void 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 
198 void 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 
226 void 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 
239 double 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 
302 double* 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 
352 double HybridPlot::GetMedian(TH1* histo){
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 }
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:785
Double_t GetX1() const
Definition: TLine.h:49
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3244
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:7807
virtual void SetName(const char *name="")
Definition: TPave.h:72
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8397
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
Definition: Rtypes.h:59
double * GetHistoPvals(TH1 *histo, double percentage)
Get the "effective sigmas" of the histo.
Definition: HybridPlot.cxx:302
~HybridPlot()
Destructor.
Definition: HybridPlot.cxx:124
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:452
THist< 1, float, THistStatContent, THistStatUncertainty > TH1F
Definition: THist.hxx:285
Definition: Rtypes.h:58
virtual TH1 * DrawNormalized(Option_t *option="", Double_t norm=1) const
Draw a normalized copy of this histogram.
Definition: TH1.cxx:3044
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7232
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4763
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:6892
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:567
void Draw(const char *options="")
Draw on current pad.
Definition: HybridPlot.cxx:138
Double_t GetRMS(Int_t axis=1) const
Definition: TH1.h:310
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1222
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
void DumpToImage(const char *filename)
Write an image on disk.
Definition: HybridPlot.cxx:226
STL namespace.
double sqrt(double)
Double_t GetXmin() const
Definition: TAxis.h:133
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
const Double_t sigma
This class provides the plots for the result of a study performed with the HybridCalculatorOriginal c...
Definition: HybridPlot.h:36
virtual Double_t GetSkewness(Int_t axis=1) const
Definition: TH1.cxx:6997
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
TVirtualPad * fPad
Definition: HybridPlot.h:116
virtual Double_t * GetIntegral()
Return a pointer to the array of bins integral.
Definition: TH1.cxx:2520
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
auto * a
Definition: textangle.C:12
void DumpToFile(const char *RootFileName, const char *options)
All the objects are written to rootfile.
Definition: HybridPlot.cxx:198
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:8477
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
A simple line.
Definition: TLine.h:23
TLine * fData_testStat_line
Definition: HybridPlot.h:114
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:359
double GetMedian(TH1 *histo)
Get the median of an histogram.
Definition: HybridPlot.cxx:352
The TH1 histogram class.
Definition: TH1.h:56
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4178
double GetHistoCenter(TH1 *histo, double n_rms=1, bool display_result=false)
Get the center of the histo.
Definition: HybridPlot.cxx:239
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:496
1-Dim function class
Definition: TF1.h:211
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2662
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
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:1266
#define gPad
Definition: TVirtualPad.h:285
Definition: Rtypes.h:59
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:618
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6154
virtual Int_t GetNbinsX() const
Definition: TH1.h:291
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
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:3688
Double_t GetXmax() const
Definition: TAxis.h:134
char name[80]
Definition: TGX11.cxx:109
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:315
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual void Print(const char *filename="") const =0
Print function.