Logo ROOT   master
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 "TVirtualPad.h"
39 #include "Math/DistFuncMathCore.h"
40 
41 #include <cmath>
42 
43 using namespace std;
44 
46 
47 using namespace RooStats;
48 
49 ////////////////////////////////////////////////////////////////////////////////
50 /// constructor from a HypoTestInverterResult class
51 /// name and title are taken from the result class
52 
53 HypoTestInverterPlot::HypoTestInverterPlot(HypoTestInverterResult* results ) :
54  TNamed( results->GetName(), results->GetTitle() ),
55  fResults(results)
56 {
57 }
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// constructor with name and title from a HypoTestInverterResult class
61 
63  const char* title,
64  HypoTestInverterResult* results ) :
65  TNamed( TString(name), TString(title) ),
66  fResults(results)
67 {
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// Make the plot of the result of the scan
72 /// using the observed data
73 /// By default plot CLs or CLsb depending if the flag UseCLs is set
74 ///
75 /// - If Option = "CLb" return CLb plot
76 /// - = "CLs+b" return CLs+b plot independently of the flag
77 /// - = "CLs" return CLs plot independently of the flag
78 
80 {
81  TString option(opt);
82  option.ToUpper();
83  int type = 0; // use default
84  if (option.Contains("CLB")) type = 1; // CLb
85  else if (option.Contains("CLS+B") || option.Contains("CLSPLUSB")) type = 2; // CLs+b
86  else if (option.Contains("CLS" )) type = 3; // CLs
87 
88  const int nEntries = fResults->ArraySize();
89 
90  // sort the arrays based on the x values
91  std::vector<unsigned int> index(nEntries);
92  TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
93 
94  // copy result in sorted arrays
95  std::vector<Double_t> xArray(nEntries);
96  std::vector<Double_t> yArray(nEntries);
97  std::vector<Double_t> yErrArray(nEntries);
98  for (int i=0; i<nEntries; i++) {
99  xArray[i] = fResults->GetXValue(index[i]);
100  if (type == 0) {
101  yArray[i] = fResults->GetYValue(index[i]);
102  yErrArray[i] = fResults->GetYError(index[i]);
103  } else if (type == 1) {
104  yArray[i] = fResults->CLb(index[i]);
105  yErrArray[i] = fResults->CLbError(index[i]);
106  } else if (type == 2) {
107  yArray[i] = fResults->CLsplusb(index[i]);
108  yErrArray[i] = fResults->CLsplusbError(index[i]);
109  } else if (type == 3) {
110  yArray[i] = fResults->CLs(index[i]);
111  yErrArray[i] = fResults->CLsError(index[i]);
112  }
113  }
114 
115  TGraphErrors* graph = new TGraphErrors(nEntries,&xArray.front(),&yArray.front(),0,&yErrArray.front());
116  TString pValueName = "CLs";
117  if (type == 1 ) pValueName = "CLb";
118  if (type == 2 || (type == 0 && !fResults->fUseCLs) ) pValueName = "CLs+b";
119  TString name = pValueName + TString("_observed");
120  TString title = TString("Observed ") + pValueName;
121  graph->SetName(name);
122  graph->SetTitle(title);
123  graph->SetMarkerStyle(20);
124  graph->SetLineWidth(2);
125  return graph;
126 }
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 /// Make the expected plot and the bands
130 /// nsig1 and nsig2 indicates the n-sigma value for the bands
131 /// if nsig1 = 0 no band is drawn (only expected value)
132 /// if nsig2 > nsig1 (default is nsig1=1 and nsig2=2) the second band is also drawn
133 /// The first band is drawn in green while the second in yellow
134 /// THe return result is a TMultiGraph object
135 
137 {
138  const int nEntries = fResults->ArraySize();
139  bool doFirstBand = (nsig1 > 0);
140  bool doSecondBand = (nsig2 > nsig1);
141 
142  nsig1 = std::abs(nsig1);
143  nsig2 = std::abs(nsig2);
144 
145  // sort the arrays based on the x values
146  std::vector<unsigned int> index(nEntries);
147  TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
148 
149  // create the graphs
150  TGraph * g0 = new TGraph;
151  TString pValueName = "CLs";
152  if (!fResults->fUseCLs) pValueName = "CLs+b";
153  g0->SetTitle(TString::Format("Expected %s - Median",pValueName.Data()) );
154  TGraphAsymmErrors * g1 = 0;
155  TGraphAsymmErrors * g2 = 0;
156  if (doFirstBand) {
157  g1 = new TGraphAsymmErrors;
158  if (nsig1 - int(nsig1) < 0.01)
159  g1->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig1)) );
160  else
161  g1->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig1) );
162  }
163  if (doSecondBand) {
164  g2 = new TGraphAsymmErrors;
165  if (nsig2 - int(nsig2) < 0.01)
166  g2->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig2)) );
167  else
168  g2->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig2) );
169  }
170  double p[7];
171  double q[7];
172  p[0] = ROOT::Math::normal_cdf(-nsig2);
173  p[1] = ROOT::Math::normal_cdf(-nsig1);
174  p[2] = 0.5;
175  p[3] = ROOT::Math::normal_cdf(nsig1);
176  p[4] = ROOT::Math::normal_cdf(nsig2);
177 
178  bool resultIsAsymptotic = ( !fResults->GetNullTestStatDist(0) && !fResults->GetAltTestStatDist(0) );
179  int np = 0;
180  for (int j=0; j<nEntries; ++j) {
181  int i = index[j]; // i is the order index
183  if ( !s) continue;
184  const std::vector<double> & values = s->GetSamplingDistribution();
185  // special case for asymptotic results (cannot use TMath::quantile in that case)
186  if (resultIsAsymptotic) {
187  double maxSigma = fResults->fgAsymptoticMaxSigma;
188  double dsig = 2* maxSigma/ (values.size() -1) ;
189  int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5);
190  int i1 = (int) TMath::Floor ( (-nsig1 + maxSigma )/dsig + 0.5);
191  int i2 = (int) TMath::Floor ( ( maxSigma)/dsig + 0.5);
192  int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma)/dsig + 0.5);
193  int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma)/dsig + 0.5);
194  q[0] = values[i0];
195  q[1] = values[i1];
196  q[2] = values[i2];
197  q[3] = values[i3];
198  q[4] = values[i4];
199  }
200  else {
201  double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
202  TMath::Quantiles(values.size(), 5, x,q,p,false);
203  }
204 
205  g0->SetPoint(np, fResults->GetXValue(i), q[2]);
206  if (g1) {
207  g1->SetPoint(np, fResults->GetXValue(i), q[2]);
208  g1->SetPointEYlow(np, q[2] - q[1]); // -1 sigma errorr
209  g1->SetPointEYhigh(np, q[3] - q[2]);//+1 sigma error
210  }
211  if (g2) {
212  g2->SetPoint(np, fResults->GetXValue(i), q[2]);
213 
214  g2->SetPointEYlow(np, q[2]-q[0]); // -2 -- -1 sigma error
215  g2->SetPointEYhigh(np, q[4]-q[2]);
216  }
217  if (s) delete s;
218  np++;
219  }
220 
221  TString name = GetName() + TString("_expected");
222  TString title = TString("Expected ") + GetTitle();
223  TMultiGraph* graph = new TMultiGraph(name,title);
224 
225  // set the graphics options and add in multi graph
226  // orderof adding is drawing order
227  if (g2) {
228  g2->SetFillColor(kYellow);
229  graph->Add(g2,"3");
230  }
231  if (g1) {
232  g1->SetFillColor(kGreen);
233  graph->Add(g1,"3");
234  }
235  g0->SetLineStyle(2);
236  g0->SetLineWidth(2);
237  graph->Add(g0,"L");
238 
239  return graph;
240 }
241 
242 ////////////////////////////////////////////////////////////////////////////////
243 /// destructor
244 
246 {
247 }
248 
249 ////////////////////////////////////////////////////////////////////////////////
250 /// Draw the result in the current canvas
251 /// Possible options:
252 /// - SAME : draw in the current axis
253 /// - OBS : draw only the observed plot
254 /// - EXP : draw only the expected plot
255 ///
256 /// - CLB : draw also the CLB
257 /// - 2CL : draw both clsplusb and cls
258 ///
259 /// default draw observed + expected with 1 and 2 sigma bands
260 
262 
263  TString option(opt);
264  option.ToUpper();
265  bool drawAxis = !option.Contains("SAME");
266  bool drawObs = option.Contains("OBS") || !option.Contains("EXP");
267  bool drawExp = option.Contains("EXP") || !option.Contains("OBS");
268  bool drawCLb = option.Contains("CLB");
269  bool draw2CL = option.Contains("2CL");
270 
271  TGraphErrors * gobs = 0;
272  TGraph * gplot = 0;
273  if (drawObs) {
274  gobs = MakePlot();
275  // add object to top-level directory to avoid mem leak
276  if (gROOT) gROOT->Add(gobs);
277  if (drawAxis) {
278  gobs->Draw("APL");
279  gplot = gobs;
280  gplot->GetHistogram()->SetTitle( GetTitle() );
281  }
282  else gobs->Draw("PL");
283 
284  }
285  TMultiGraph * gexp = 0;
286  if (drawExp) {
287  gexp = MakeExpectedPlot();
288  // add object to current directory to avoid mem leak
289  if (gROOT) gROOT->Add(gexp);
290  if (drawAxis && !drawObs) {
291  gexp->Draw("A");
292  if (gexp->GetHistogram()) gexp->GetHistogram()->SetTitle( GetTitle() );
293  gplot = (TGraph*) gexp->GetListOfGraphs()->First();
294  }
295  else
296  gexp->Draw();
297 
298  }
299 
300  // draw also an horizontal line at the desired conf level
301  if (gplot) {
302  double alpha = 1.-fResults->ConfidenceLevel();
303  double x1 = gplot->GetXaxis()->GetXmin();
304  double x2 = gplot->GetXaxis()->GetXmax();
305  TLine * line = new TLine(x1, alpha, x2,alpha);
307  line->Draw();
308  // put axis labels
309  RooAbsArg * arg = fResults->fParameters.first();
310  if (arg) gplot->GetXaxis()->SetTitle(arg->GetName());
311  gplot->GetYaxis()->SetTitle("p value");
312  }
313 
314 
315  TGraph *gclb = 0;
316  if (drawCLb) {
317  gclb = MakePlot("CLb");
318  if (gROOT) gROOT->Add(gclb);
319  gclb->SetMarkerColor(kBlue+4);
320  gclb->Draw("PL");
321  // draw in red observed cls or clsb
322  if (gobs) gobs->SetMarkerColor(kRed);
323  }
324  TGraph * gclsb = 0;
325  TGraph * gcls = 0;
326  if (draw2CL) {
327  if (fResults->fUseCLs) {
328  gclsb = MakePlot("CLs+b");
329  if (gROOT) gROOT->Add(gclsb);
330  gclsb->SetMarkerColor(kBlue);
331  gclsb->Draw("PL");
332  gclsb->SetLineStyle(3);
333  }
334  else {
335  gcls = MakePlot("CLs");
336  if (gROOT) gROOT->Add(gcls);
337  gcls->SetMarkerColor(kBlue);
338  gcls->Draw("PL");
339  gcls->SetLineStyle(3);
340  }
341  }
342  // draw again observed values otherwise will be covered by the bands
343  if (gobs) {
344  gobs->Draw("PL");
345  }
346 
347 
348  double y0 = 0.6;
349  double verticalSize = (gexp || draw2CL || drawCLb ) ? 0.3 : 0.15;
350  double y1 = y0 + verticalSize;
351  TLegend * l = new TLegend(0.6,y0,0.9,y1);
352  if (gobs) l->AddEntry(gobs,"","PEL");
353  if (gclsb) l->AddEntry(gclsb,"","PEL");
354  if (gcls) l->AddEntry(gcls,"","PEL");
355  if (gclb) l->AddEntry(gclb,"","PEL");
356  if (gexp) {
357  // loop in reverse order (opposite to drawing one)
358  int ngraphs = gexp->GetListOfGraphs()->GetSize();
359  for (int i = ngraphs-1; i>=0; --i) {
360  TObject * obj = gexp->GetListOfGraphs()->At(i);
361  TString lopt = "F";
362  if (i == ngraphs-1) lopt = "L";
363  if (obj) l->AddEntry(obj,"",lopt);
364  }
365  }
366  l->Draw();
367  // redraw the axis
368  if (gPad) gPad->RedrawAxis();
369 
370 }
371 
372 ////////////////////////////////////////////////////////////////////////////////
373 /// \param index Index of the result stored in HypoTestInverterResult
374 /// \param type Type of the test (see below)
375 /// \param nbins Number of bins
376 /// - type =0 null and alt
377 /// - type = 1 only null (S+B)
378 /// - type = 2 only alt (B)
379 
381  SamplingDistPlot * pl = 0;
382  if (type == 0) {
383  HypoTestResult * result = (HypoTestResult*) fResults->fYObjects.At(index);
384  if (result)
385  pl = new HypoTestPlot(*result, nbins );
386  return pl;
387  }
388  if (type == 1) {
390  if (sbDist) {
391  pl = new SamplingDistPlot( nbins);
392  pl->AddSamplingDistribution(sbDist);
393  return pl;
394  }
395  }
396  if (type == 2) {
398  if (bDist) {
399  pl = new SamplingDistPlot( nbins);
400  pl->AddSamplingDistribution(bDist);
401  return pl;
402  }
403  }
404  return 0;
405 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:150
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
get the signal and background test statistic distribution
Double_t Floor(Double_t x)
Definition: TMath.h:693
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:64
Definition: Rtypes.h:64
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:36
TAxis * GetYaxis() const
Get y axis of the graph.
Definition: TGraph.cxx:1628
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
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:405
Basic string class.
Definition: TString.h:131
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:195
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="")
Change (i.e.
Definition: TGraph.cxx:2316
Definition: Rtypes.h:64
Definition: Rtypes.h:64
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:752
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:2311
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.
Definition: TMath.cxx:1183
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
static constexpr double s
TH1F * GetHistogram()
Returns a pointer to the histogram used to draw the axis.
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.
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:658
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:22
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:356
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:1618
Namespace for the RooStats classes.
Definition: Asimov.h:19
static const double x1[5]
#define ClassImp(name)
Definition: Rtypes.h:361
TH1F * GetHistogram() const
Returns a pointer to the histogram used to draw the axis Takes into account the two following cases...
Definition: TGraph.cxx:1474
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator.
int type
Definition: TGX11.cxx:120
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
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:2261
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:70
void Draw(Option_t *opt="")
Draw the scan result in the current canvas Possible options: "" (default): draw observed + expected w...
auto * l
Definition: textangle.C:4
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
A TGraph is an 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:287
Definition: Rtypes.h:64
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6345
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
float * q
Definition: THbookFile.cxx:87
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
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
char name[80]
Definition: TGX11.cxx:109
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: TMathBase.h:337
std::vector< double > fXValues
number of points used to build expected p-values
const char * Data() const
Definition: TString.h:364