Logo ROOT   6.10/09
Reference Guide
HypoTestInverterPlot.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3 /*************************************************************************
4  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5  * All rights reserved. *
6  * *
7  * For the licensing terms see $ROOTSYS/LICENSE. *
8  * For the list of contributors see $ROOTSYS/README/CREDITS. *
9  *************************************************************************/
10 
11 /** \class RooStats::HypoTestInverterPlot
12  \ingroup Roostats
13 
14 Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator
15 
16 It can be used to plot the obtained p-values ( CLb, CLs+b or CLs) for each scanned point and
17 the test statistic distributions (when a calculator based on pseudo-experiments is used) for the two
18 hypotheses.
19 
20 */
21 
22 // include header file of this class
24 
25 // include other header files
26 #include "RooStats/HybridResult.h"
28 #include "RooStats/HypoTestPlot.h"
30 
31 #include "TGraphErrors.h"
32 #include "TGraphAsymmErrors.h"
33 #include "TMultiGraph.h"
34 #include "TROOT.h"
35 #include "TLine.h"
36 #include "TAxis.h"
37 #include "TLegend.h"
38 #include "TH1.h"
39 #include "TPad.h"
40 #include "Math/DistFuncMathCore.h"
41 
42 #include <cmath>
43 
44 using namespace std;
45 
47 
48 using namespace RooStats;
49 
50 ////////////////////////////////////////////////////////////////////////////////
51 /// constructor from a HypoTestInverterResult class
52 /// name and title are taken from the result class
53 
54 HypoTestInverterPlot::HypoTestInverterPlot(HypoTestInverterResult* results ) :
55  TNamed( results->GetName(), results->GetTitle() ),
56  fResults(results)
57 {
58 }
59 
60 ////////////////////////////////////////////////////////////////////////////////
61 /// constructor with name and title from a HypoTestInverterResult class
62 
64  const char* title,
65  HypoTestInverterResult* results ) :
66  TNamed( TString(name), TString(title) ),
67  fResults(results)
68 {
69 }
70 
71 ////////////////////////////////////////////////////////////////////////////////
72 /// Make the plot of the result of the scan
73 /// using the observed data
74 /// By default plot CLs or CLsb depending if the flag UseCLs is set
75 ///
76 /// - If Option = "CLb" return CLb plot
77 /// - = "CLs+b" return CLs+b plot independently of the flag
78 /// - = "CLs" return CLs plot independently of the flag
79 
81 {
82  TString option(opt);
83  option.ToUpper();
84  int type = 0; // use default
85  if (option.Contains("CLB")) type = 1; // CLb
86  else if (option.Contains("CLS+B") || option.Contains("CLSPLUSB")) type = 2; // CLs+b
87  else if (option.Contains("CLS" )) type = 3; // CLs
88 
89  const int nEntries = fResults->ArraySize();
90 
91  // sort the arrays based on the x values
92  std::vector<unsigned int> index(nEntries);
93  TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
94 
95  // copy result in sorted arrays
96  std::vector<Double_t> xArray(nEntries);
97  std::vector<Double_t> yArray(nEntries);
98  std::vector<Double_t> yErrArray(nEntries);
99  for (int i=0; i<nEntries; i++) {
100  xArray[i] = fResults->GetXValue(index[i]);
101  if (type == 0) {
102  yArray[i] = fResults->GetYValue(index[i]);
103  yErrArray[i] = fResults->GetYError(index[i]);
104  } else if (type == 1) {
105  yArray[i] = fResults->CLb(index[i]);
106  yErrArray[i] = fResults->CLbError(index[i]);
107  } else if (type == 2) {
108  yArray[i] = fResults->CLsplusb(index[i]);
109  yErrArray[i] = fResults->CLsplusbError(index[i]);
110  } else if (type == 3) {
111  yArray[i] = fResults->CLs(index[i]);
112  yErrArray[i] = fResults->CLsError(index[i]);
113  }
114  }
115 
116  TGraphErrors* graph = new TGraphErrors(nEntries,&xArray.front(),&yArray.front(),0,&yErrArray.front());
117  TString pValueName = "CLs";
118  if (type == 1 ) pValueName = "CLb";
119  if (type == 2 || (type == 0 && !fResults->fUseCLs) ) pValueName = "CLs+b";
120  TString name = pValueName + TString("_observed");
121  TString title = TString("Observed ") + pValueName;
122  graph->SetName(name);
123  graph->SetTitle(title);
124  graph->SetMarkerStyle(20);
125  graph->SetLineWidth(2);
126  return graph;
127 }
128 
129 ////////////////////////////////////////////////////////////////////////////////
130 /// Make the expected plot and the bands
131 /// nsig1 and nsig2 indicates the n-sigma value for the bands
132 /// if nsig1 = 0 no band is drawn (only expected value)
133 /// if nsig2 > nsig1 (default is nsig1=1 and nsig2=2) the second band is also drawn
134 /// The first band is drawn in green while the second in yellow
135 /// THe return result is a TMultiGraph object
136 
138 {
139  const int nEntries = fResults->ArraySize();
140  bool doFirstBand = (nsig1 > 0);
141  bool doSecondBand = (nsig2 > nsig1);
142 
143  nsig1 = std::abs(nsig1);
144  nsig2 = std::abs(nsig2);
145 
146  // sort the arrays based on the x values
147  std::vector<unsigned int> index(nEntries);
148  TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
149 
150  // create the graphs
151  TGraph * g0 = new TGraph;
152  TString pValueName = "CLs";
153  if (!fResults->fUseCLs) pValueName = "CLs+b";
154  g0->SetTitle(TString::Format("Expected %s - Median",pValueName.Data()) );
155  TGraphAsymmErrors * g1 = 0;
156  TGraphAsymmErrors * g2 = 0;
157  if (doFirstBand) {
158  g1 = new TGraphAsymmErrors;
159  if (nsig1 - int(nsig1) < 0.01)
160  g1->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig1)) );
161  else
162  g1->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig1) );
163  }
164  if (doSecondBand) {
165  g2 = new TGraphAsymmErrors;
166  if (nsig2 - int(nsig2) < 0.01)
167  g2->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig2)) );
168  else
169  g2->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig2) );
170  }
171  double p[7];
172  double q[7];
173  p[0] = ROOT::Math::normal_cdf(-nsig2);
174  p[1] = ROOT::Math::normal_cdf(-nsig1);
175  p[2] = 0.5;
176  p[3] = ROOT::Math::normal_cdf(nsig1);
177  p[4] = ROOT::Math::normal_cdf(nsig2);
178 
179  bool resultIsAsymptotic = ( !fResults->GetNullTestStatDist(0) && !fResults->GetAltTestStatDist(0) );
180  int np = 0;
181  for (int j=0; j<nEntries; ++j) {
182  int i = index[j]; // i is the order index
184  if ( !s) continue;
185  const std::vector<double> & values = s->GetSamplingDistribution();
186  // special case for asymptotic results (cannot use TMath::quantile in that case)
187  if (resultIsAsymptotic) {
188  double maxSigma = fResults->fgAsymptoticMaxSigma;
189  double dsig = 2* maxSigma/ (values.size() -1) ;
190  int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5);
191  int i1 = (int) TMath::Floor ( (-nsig1 + maxSigma )/dsig + 0.5);
192  int i2 = (int) TMath::Floor ( ( maxSigma)/dsig + 0.5);
193  int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma)/dsig + 0.5);
194  int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma)/dsig + 0.5);
195  q[0] = values[i0];
196  q[1] = values[i1];
197  q[2] = values[i2];
198  q[3] = values[i3];
199  q[4] = values[i4];
200  }
201  else {
202  double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
203  TMath::Quantiles(values.size(), 5, x,q,p,false);
204  }
205 
206  g0->SetPoint(np, fResults->GetXValue(i), q[2]);
207  if (g1) {
208  g1->SetPoint(np, fResults->GetXValue(i), q[2]);
209  g1->SetPointEYlow(np, q[2] - q[1]); // -1 sigma errorr
210  g1->SetPointEYhigh(np, q[3] - q[2]);//+1 sigma error
211  }
212  if (g2) {
213  g2->SetPoint(np, fResults->GetXValue(i), q[2]);
214 
215  g2->SetPointEYlow(np, q[2]-q[0]); // -2 -- -1 sigma error
216  g2->SetPointEYhigh(np, q[4]-q[2]);
217  }
218  if (s) delete s;
219  np++;
220  }
221 
222  TString name = GetName() + TString("_expected");
223  TString title = TString("Expected ") + GetTitle();
224  TMultiGraph* graph = new TMultiGraph(name,title);
225 
226  // set the graphics options and add in multi graph
227  // orderof adding is drawing order
228  if (g2) {
229  g2->SetFillColor(kYellow);
230  graph->Add(g2,"3");
231  }
232  if (g1) {
233  g1->SetFillColor(kGreen);
234  graph->Add(g1,"3");
235  }
236  g0->SetLineStyle(2);
237  g0->SetLineWidth(2);
238  graph->Add(g0,"L");
239 
240  return graph;
241 }
242 
243 ////////////////////////////////////////////////////////////////////////////////
244 /// destructor
245 
247 {
248 }
249 
250 ////////////////////////////////////////////////////////////////////////////////
251 /// Draw the result in the current canvas
252 /// Possible options:
253 /// - SAME : draw in the current axis
254 /// - OBS : draw only the observed plot
255 /// - EXP : draw only the expected plot
256 ///
257 /// - CLB : draw also the CLB
258 /// - 2CL : draw both clsplusb and cls
259 ///
260 /// default draw observed + expected with 1 and 2 sigma bands
261 
263 
264  TString option(opt);
265  option.ToUpper();
266  bool drawAxis = !option.Contains("SAME");
267  bool drawObs = option.Contains("OBS") || !option.Contains("EXP");
268  bool drawExp = option.Contains("EXP") || !option.Contains("OBS");
269  bool drawCLb = option.Contains("CLB");
270  bool draw2CL = option.Contains("2CL");
271 
272  TGraphErrors * gobs = 0;
273  TGraph * gplot = 0;
274  if (drawObs) {
275  gobs = MakePlot();
276  // add object to top-level directory to avoid mem leak
277  if (gROOT) gROOT->Add(gobs);
278  if (drawAxis) {
279  gobs->Draw("APL");
280  gplot = gobs;
281  gplot->GetHistogram()->SetTitle( GetTitle() );
282  }
283  else gobs->Draw("PL");
284 
285  }
286  TMultiGraph * gexp = 0;
287  if (drawExp) {
288  gexp = MakeExpectedPlot();
289  // add object to current directory to avoid mem leak
290  if (gROOT) gROOT->Add(gexp);
291  if (drawAxis && !drawObs) {
292  gexp->Draw("A");
293  if (gexp->GetHistogram()) gexp->GetHistogram()->SetTitle( GetTitle() );
294  gplot = (TGraph*) gexp->GetListOfGraphs()->First();
295  }
296  else
297  gexp->Draw();
298 
299  }
300 
301  // draw also an horizontal line at the desired conf level
302  if (gplot) {
303  double alpha = 1.-fResults->ConfidenceLevel();
304  double x1 = gplot->GetXaxis()->GetXmin();
305  double x2 = gplot->GetXaxis()->GetXmax();
306  TLine * line = new TLine(x1, alpha, x2,alpha);
307  line->SetLineColor(kRed);
308  line->Draw();
309  // put axis labels
310  RooAbsArg * arg = fResults->fParameters.first();
311  if (arg) gplot->GetXaxis()->SetTitle(arg->GetName());
312  gplot->GetYaxis()->SetTitle("p value");
313  }
314 
315 
316  TGraph *gclb = 0;
317  if (drawCLb) {
318  gclb = MakePlot("CLb");
319  if (gROOT) gROOT->Add(gclb);
320  gclb->SetMarkerColor(kBlue+4);
321  gclb->Draw("PL");
322  // draw in red observed cls or clsb
323  if (gobs) gobs->SetMarkerColor(kRed);
324  }
325  TGraph * gclsb = 0;
326  TGraph * gcls = 0;
327  if (draw2CL) {
328  if (fResults->fUseCLs) {
329  gclsb = MakePlot("CLs+b");
330  if (gROOT) gROOT->Add(gclsb);
331  gclsb->SetMarkerColor(kBlue);
332  gclsb->Draw("PL");
333  gclsb->SetLineStyle(3);
334  }
335  else {
336  gcls = MakePlot("CLs");
337  if (gROOT) gROOT->Add(gcls);
338  gcls->SetMarkerColor(kBlue);
339  gcls->Draw("PL");
340  gcls->SetLineStyle(3);
341  }
342  }
343  // draw again observed values otherwise will be covered by the bands
344  if (gobs) {
345  gobs->Draw("PL");
346  }
347 
348 
349  double y0 = 0.6;
350  double verticalSize = (gexp || draw2CL || drawCLb ) ? 0.3 : 0.15;
351  double y1 = y0 + verticalSize;
352  TLegend * l = new TLegend(0.6,y0,0.9,y1);
353  if (gobs) l->AddEntry(gobs,"","PEL");
354  if (gclsb) l->AddEntry(gclsb,"","PEL");
355  if (gcls) l->AddEntry(gcls,"","PEL");
356  if (gclb) l->AddEntry(gclb,"","PEL");
357  if (gexp) {
358  // loop in reverse order (opposite to drawing one)
359  int ngraphs = gexp->GetListOfGraphs()->GetSize();
360  for (int i = ngraphs-1; i>=0; --i) {
361  TObject * obj = gexp->GetListOfGraphs()->At(i);
362  TString lopt = "F";
363  if (i == ngraphs-1) lopt = "L";
364  if (obj) l->AddEntry(obj,"",lopt);
365  }
366  }
367  l->Draw();
368  // redraw the axis
369  if (gPad) gPad->RedrawAxis();
370 
371 }
372 
373 ////////////////////////////////////////////////////////////////////////////////
374 /// plot the test statistic distributions
375 /// - type =0 null and alt
376 /// - type = 1 only null (S+B)
377 /// - type = 2 only alt (B)
378 
380  SamplingDistPlot * pl = 0;
381  if (type == 0) {
383  if (result)
384  pl = new HypoTestPlot(*result, nbins );
385  return pl;
386  }
387  if (type == 1) {
389  if (sbDist) {
390  pl = new SamplingDistPlot( nbins);
391  pl->AddSamplingDistribution(sbDist);
392  return pl;
393  }
394  }
395  if (type == 2) {
397  if (bDist) {
398  pl = new SamplingDistPlot( nbins);
399  pl->AddSamplingDistribution(bDist);
400  return pl;
401  }
402  }
403  return 0;
404 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:145
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
get the signal and background test statistic distribution
Double_t Floor(Double_t x)
Definition: TMath.h:600
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
SamplingDistribution * GetExpectedPValueDist(int index) const
return expected distribution of p-values (Cls or Clsplusb)
TLine * line
const char Option_t
Definition: RtypesCore.h:62
Definition: Rtypes.h:56
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:131
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:452
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
HypoTestResult is a base class for results from hypothesis tests.
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:35
TAxis * GetYaxis() const
Get y axis of the graph.
Definition: TGraph.cxx:1602
Double_t AddSamplingDistribution(const SamplingDistribution *samplingDist, Option_t *drawOptions="NORMALIZE HIST")
adds the sampling distribution and returns the scale factor
#define gROOT
Definition: TROOT.h:375
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:202
double GetYError(int index) const
function to return the estimated error on the value of the confidence level for the i^th entry in the...
virtual void SetTitle(const char *title="")
Set graph title.
Definition: TGraph.cxx:2180
int nbins[3]
Definition: Rtypes.h:56
Definition: Rtypes.h:56
STL namespace.
TMultiGraph * MakeExpectedPlot(double sig1=1, double sig2=2)
Make the expected plot and the bands nsig1 and nsig2 indicates the n-sigma value for the bands if nsi...
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:745
HypoTestInverterPlot(HypoTestInverterResult *results)
constructor from a HypoTestInverterResult class name and title are taken from the result class ...
HypoTestInverterResult * fResults
TGraph with asymmetric error bars.
Double_t GetXmin() const
Definition: TAxis.h:133
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2345
double CLsError(int index) const
return the observed CLb value for the i-th entry
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities Parameters: x -the data sample n ...
Definition: TMath.cxx:1177
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
virtual void SetPointEYlow(Int_t i, Double_t eyl)
Set EYlow for point i.
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
double CLsplusb(int index) const
return the observed CLsplusb value for the i-th entry
virtual void Draw(Option_t *chopt="")
Draw this multigraph with its current attributes.
double CLs(int index) const
return the observed CLb value for the i-th entry
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
SamplingDistPlot * MakeTestStatPlot(int index, int type=0, int nbins=100)
Plot the test statistic distributions type =0 null and alt type = 1 only null (S+B) type = 2 only alt...
TH1F * GetHistogram() const
Returns a pointer to the histogram used to draw the axis.
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:561
SamplingDistribution * GetBackgroundTestStatDist(int index) const
get the background test statistic distribution
TGraphErrors * MakePlot(Option_t *opt="")
return a TGraphErrors with the obtained observed p-values resultinf from the scan By default (Option ...
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
RooAbsArg * first() const
virtual Double_t ConfidenceLevel() const
return confidence level
double CLbError(int index) const
return the observed CLb value for the i-th entry
double CLb(int index) const
return the observed CLb value for the i-th entry
A simple line.
Definition: TLine.h:23
TLine * l
Definition: textangle.C:4
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:315
Definition: graph.py:1
This class simply holds a sampling distribution of some test statistic.
TAxis * GetXaxis() const
Get x axis of the graph.
Definition: TGraph.cxx:1592
Namespace for the RooStats classes.
Definition: Asimov.h:20
static const double x1[5]
#define ClassImp(name)
Definition: Rtypes.h:336
TH1F * GetHistogram() const
Returns a pointer to the histogram used to draw the axis Takes into account the two following cases...
Definition: TGraph.cxx:1469
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:359
Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator.
int type
Definition: TGX11.cxx:120
double GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results ...
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:42
This class provides simple and straightforward utilities to plot SamplingDistribution objects...
Mother of all ROOT objects.
Definition: TObject.h:37
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2156
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:69
void Draw(Option_t *opt="")
Draw the scan result in the current canvas Possible options: "" (default): draw observed + expected w...
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
#define gPad
Definition: TVirtualPad.h:284
double result[121]
Definition: Rtypes.h:56
virtual void SetTitle(const char *title)
Change (i.e.
Definition: TH1.cxx:6028
virtual void Add(TGraph *graph, Option_t *chopt="")
Add a new graph to the list of graphs.
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
virtual Int_t GetSize() const
Definition: TCollection.h:89
float * q
Definition: THbookFile.cxx:87
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
double CLsplusbError(int index) const
return the observed CLsplusb value for the i-th entry
int ArraySize() const
number of entries in the results array
Double_t GetXmax() const
Definition: TAxis.h:134
virtual void SetPointEYhigh(Int_t i, Double_t eyh)
Set EYhigh for point i.
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
SamplingDistribution * GetAltTestStatDist(int index) const
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Definition: TMath.h:1126
std::vector< double > fXValues
number of points used to build expected p-values