Logo ROOT   6.08/07
Reference Guide
HybridPlot.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 
3 ////////////////////////////////////////////////////////////////////////////////
4 
5 
6 
7 #include "assert.h"
8 #include <cmath>
9 #include <iostream>
10 #include <map>
11 
12 #include "RooStats/HybridPlot.h"
13 #include "TStyle.h"
14 #include "TF1.h"
15 #include "TAxis.h"
16 #include "TH1.h"
17 #include "TLine.h"
18 #include "TLegend.h"
19 #include "TFile.h"
20 #include "TVirtualPad.h"
21 
22 #include <algorithm>
23 
24 /// To build the THtml documentation
25 using namespace std;
26 
28 
29 using namespace RooStats;
30 
31 /*----------------------------------------------------------------------------*/
32 
33 HybridPlot::HybridPlot(const char* name,
34  const char* title,
35  const std::vector<double> & sb_vals,
36  const std::vector<double> & b_vals,
37  double testStat_data,
38  int n_bins,
39  bool verbosity):
40  TNamed(name,title),
41  fSb_histo(NULL),
42  fSb_histo_shaded(NULL),
43  fB_histo(NULL),
44  fB_histo_shaded(NULL),
45  fData_testStat_line(0),
46  fLegend(0),
47  fPad(0),
48  fVerbose(verbosity)
49 {
50  // HybridPlot constructor
51 
52  int nToysSB = sb_vals.size();
53  int nToysB = sb_vals.size();
54  assert (nToysSB >0);
55  assert (nToysB >0);
56 
57  // Get the max and the min of the plots
58  double min = *std::min_element(sb_vals.begin(), sb_vals.end());
59  double max = *std::max_element(sb_vals.begin(), sb_vals.end());
60 
61  double min_b = *std::min_element(b_vals.begin(), b_vals.end());
62  double max_b = *std::max_element(b_vals.begin(), b_vals.end());
63 
64 
65  if ( min_b < min) min = min_b;
66  if ( max_b > max) max = max_b;
67 
68  if (testStat_data<min) min = testStat_data;
69  if (testStat_data>max) max = testStat_data;
70 
71  min *= 1.1;
72  max *= 1.1;
73 
74  // Build the histos
75 
76  fSb_histo = new TH1F ("SB_model",title,n_bins,min,max);
79  fSb_histo->GetXaxis()->SetTitle("test statistics");
81 
82  fB_histo = new TH1F ("B_model",title,n_bins,min,max);
85  fB_histo->GetXaxis()->SetTitle("test statistics");
87 
88  for (int i=0;i<nToysSB;++i) fSb_histo->Fill(sb_vals[i]);
89  for (int i=0;i<nToysB;++i) fB_histo->Fill(b_vals[i]);
90 
91  double histos_max_y = fSb_histo->GetMaximum();
92  double line_hight = histos_max_y/nToysSB;
93  if (histos_max_y<fB_histo->GetMaximum()) histos_max_y = fB_histo->GetMaximum()/nToysB;
94 
95  // Build the line of the measured -2lnQ
96  fData_testStat_line = new TLine(testStat_data,0,testStat_data,line_hight);
99 
100  // The legend
101  double golden_section = (std::sqrt(5.)-1)/2;
102 
103  fLegend = new TLegend(0.75,0.95-0.2*golden_section,0.95,0.95);
104  TString title_leg="test statistics distributions ";
105  title_leg+=sb_vals.size();
106  title_leg+=" toys";
107  fLegend->SetName(title_leg.Data());
108  fLegend->AddEntry(fSb_histo,"SB toy datasets");
109  fLegend->AddEntry(fB_histo,"B toy datasets");
110  fLegend->AddEntry((TLine*)fData_testStat_line,"test statistics on data","L");
111  fLegend->SetFillColor(0);
112 
113 }
114 
115 /*----------------------------------------------------------------------------*/
116 
118 {
119  // destructor
120 
121  if (fSb_histo) delete fSb_histo;
122  if (fB_histo) delete fB_histo;
124  if (fB_histo_shaded) delete fB_histo_shaded;
126  if (fLegend) delete fLegend;
127 }
128 
129 /*----------------------------------------------------------------------------*/
130 
131 void HybridPlot::Draw(const char* )
132 {
133  // draw the S+B and B histogram in the current canvas
134 
135 
136  // We don't want the statistics of the histos
137  gStyle->SetOptStat(0);
138 
139  // The histos
140 
143  fB_histo->DrawNormalized("same");
144  }
145  else{
147  fSb_histo->DrawNormalized("same");
148  }
149 
150  // Shaded
151  fB_histo_shaded = (TH1F*)fB_histo->Clone("b_shaded");
154 
155  fSb_histo_shaded = (TH1F*)fSb_histo->Clone("sb_shaded");
158 
159  // Empty the bins according to the data -2lnQ
160  double data_m2lnq= fData_testStat_line->GetX1();
161  for (int i=1;i<=fSb_histo->GetNbinsX();++i){
162  if (fSb_histo->GetBinCenter(i)<data_m2lnq){
165  }
166  else{
169  }
170  }
171 
172  // Draw the shaded histos
173  fB_histo_shaded->Draw("same");
174  fSb_histo_shaded->Draw("same");
175 
176  // The line
177  fData_testStat_line->Draw("same");
178 
179  // The legend
180  fLegend->Draw("same");
181 
182  if (gPad) {
183  gPad->SetName(GetName());
184  gPad->SetTitle(GetTitle() );
185  }
186 
187  fPad = gPad;
188 
189 }
190 
191 /*----------------------------------------------------------------------------*/
192 
193 void HybridPlot::DumpToFile (const char* RootFileName, const char* options)
194 {
195  // All the objects are written to rootfile
196 
197  TFile ofile(RootFileName,options);
198  ofile.cd();
199 
200  // The histos
201  fSb_histo->Write();
202  fB_histo->Write();
203 
204  // The shaded histos
208  }
209 
210  // The line
211  fData_testStat_line->Write("Measured test statistics line tag");
212 
213  // The legend
214  fLegend->Write();
215 
216  ofile.Close();
217 
218 }
219 
220 void HybridPlot::DumpToImage(const char * filename) {
221  if (!fPad) {
222  Error("HybridPlot","Hybrid plot has not yet been drawn ");
223  return;
224  }
225  fPad->Print(filename);
226 }
227 
228 /*----------------------------------------------------------------------------*/
229 
230 // from Rsc.cxx
231 
232 
233 /**
234  Perform 2 times a gaussian fit to fetch the center of the histo.
235  To get the second fit range get an interval that tries to keep into account
236  the skewness of the distribution.
237 **/
238 double HybridPlot::GetHistoCenter(TH1* histo_orig, double n_rms, bool display_result){
239  // Get the center of the histo
240 
241  TString optfit = "Q0";
242  if (display_result) optfit = "Q";
243 
244  TH1F* histo = (TH1F*)histo_orig->Clone();
245 
246  // get the histo x extremes
247  double x_min = histo->GetXaxis()->GetXmin();
248  double x_max = histo->GetXaxis()->GetXmax();
249 
250  // First fit!
251 
252  TF1* gaus = new TF1("mygaus", "gaus", x_min, x_max);
253 
254  gaus->SetParameter("Constant",histo->GetEntries());
255  gaus->SetParameter("Mean",histo->GetMean());
256  gaus->SetParameter("Sigma",histo->GetRMS());
257 
258  histo->Fit(gaus,optfit);
259 
260  // Second fit!
261  double sigma = gaus->GetParameter("Sigma");
262  double mean = gaus->GetParameter("Mean");
263 
264  delete gaus;
265 
266  std::cout << "Center is 1st pass = " << mean << std::endl;
267 
268  double skewness = histo->GetSkewness();
269 
270  x_min = mean - n_rms*sigma - sigma*skewness/2;
271  x_max = mean + n_rms*sigma - sigma*skewness/2;;
272 
273  TF1* gaus2 = new TF1("mygaus2", "gaus", x_min, x_max);
274  gaus2->SetParameter("Mean",mean);
275 
276  // second fit : likelihood fit
277  optfit += "L";
278  histo->Fit(gaus2,optfit,"", x_min, x_max);
279 
280 
281  double center = gaus2->GetParameter("Mean");
282 
283  if (display_result) {
284  histo->Draw();
285  gaus2->Draw("same");
286  }
287  else {
288  delete histo;
289  }
290  delete gaus2;
291 
292  return center;
293 
294 
295 }
296 
297 /**
298  We let an orizzontal bar go down and we stop when we have the integral
299  equal to the desired one.
300 **/
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 /**
351  Get the median of an histogram.
352 **/
353 double HybridPlot::GetMedian(TH1* histo){
354 
355  //int xbin_median;
356  Double_t* integral = histo->GetIntegral();
357  int median_i = 0;
358  for (int j=0;j<histo->GetNbinsX(); j++)
359  if (integral[j]<0.5)
360  median_i = j;
361 
362  double median_x =
363  histo->GetBinCenter(median_i)+
364  (histo->GetBinCenter(median_i+1)-
365  histo->GetBinCenter(median_i))*
366  (0.5-integral[median_i])/(integral[median_i+1]-integral[median_i]);
367 
368  return median_x;
369 }
370 
371 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
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:830
Double_t GetX1() const
Definition: TLine.h:59
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:49
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3125
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:7664
virtual void SetName(const char *name="")
Definition: TPave.h:76
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8251
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:27
Definition: Rtypes.h:61
double * GetHistoPvals(TH1 *histo, double percentage)
Get the "effective sigmas" of the histo.
Definition: HybridPlot.cxx:302
~HybridPlot()
Destructor.
Definition: HybridPlot.cxx:117
R__EXTERN TStyle * gStyle
Definition: TStyle.h:418
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:373
THist< 1, float, THistStatContent, THistStatUncertainty > TH1F
Definition: THist.hxx:302
Definition: Rtypes.h:60
virtual TH1 * DrawNormalized(Option_t *option="", Double_t norm=1) const
Draw a normalized copy of this histogram.
Definition: TH1.cxx:2925
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7100
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4638
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:6760
tomato 1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
void Draw(const char *options="")
Draw on current pad.
Definition: HybridPlot.cxx:131
TArc * a
Definition: textangle.C:12
Double_t GetRMS(Int_t axis=1) const
Definition: TH1.h:320
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1086
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:255
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:44
void DumpToImage(const char *filename)
Write an image on disk.
Definition: HybridPlot.cxx:220
STL namespace.
double sqrt(double)
Double_t GetXmin() const
Definition: TAxis.h:139
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
Double_t x_min[]
Definition: SparseFit4.cxx:32
const Double_t sigma
This class provides the plots for the result of a study performed with the HybridCalculatorOriginal c...
Definition: HybridPlot.h:53
virtual Double_t GetSkewness(Int_t axis=1) const
Definition: TH1.cxx:6865
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:46
TVirtualPad * fPad
Definition: HybridPlot.h:133
virtual Double_t * GetIntegral()
Return a pointer to the array of bins integral.
Definition: TH1.cxx:2403
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:2851
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:42
void DumpToFile(const char *RootFileName, const char *options)
All the objects are written to rootfile.
Definition: HybridPlot.cxx:193
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:8323
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
A simple line.
Definition: TLine.h:33
TLine * fData_testStat_line
Definition: HybridPlot.h:131
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:279
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:280
double GetMedian(TH1 *histo)
Get the median of an histogram.
Definition: HybridPlot.cxx:353
The TH1 histogram class.
Definition: TH1.h:80
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4053
double GetHistoCenter(TH1 *histo, double n_rms=1, bool display_result=false)
Get the center of the histo.
Definition: HybridPlot.cxx:238
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:359
1-Dim function class
Definition: TF1.h:149
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2544
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:1257
#define NULL
Definition: Rtypes.h:82
#define gPad
Definition: TVirtualPad.h:289
Definition: Rtypes.h:61
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:431
virtual void SetTitle(const char *title)
Change (i.e.
Definition: TH1.cxx:6027
virtual Int_t GetNbinsX() const
Definition: TH1.h:301
Double_t x_max[]
Definition: SparseFit4.cxx:33
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
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:3563
Double_t GetXmax() const
Definition: TAxis.h:140
char name[80]
Definition: TGX11.cxx:109
TAxis * GetXaxis()
Definition: TH1.h:324
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
virtual void Print(const char *filename="") const =0
Print function.