Logo ROOT  
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
14Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator
15
16It can be used to plot the obtained p-values ( CLb, CLs+b or CLs) for each scanned point and
17the test statistic distributions (when a calculator based on pseudo-experiments is used) for the two
18hypotheses.
19
20*/
21
22// include header file of this class
24
25// include other header files
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"
41
42#include <cmath>
43
44using namespace std;
45
47
48using namespace RooStats;
49
50////////////////////////////////////////////////////////////////////////////////
51/// constructor from a HypoTestInverterResult class
52/// name and title are taken from the result class
53
54HypoTestInverterPlot::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) {
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);
308 line->Draw();
309 // put axis labels
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/// \param index Index of the result stored in HypoTestInverterResult
375/// \param type Type of the test (see below)
376/// \param nbins Number of bins
377/// - type =0 null and alt
378/// - type = 1 only null (S+B)
379/// - type = 2 only alt (B)
380
382 SamplingDistPlot * pl = 0;
383 if (type == 0) {
384 HypoTestResult * result = (HypoTestResult*) fResults->fYObjects.At(index);
385 if (result)
386 pl = new HypoTestPlot(*result, nbins );
387 return pl;
388 }
389 if (type == 1) {
391 if (sbDist) {
392 pl = new SamplingDistPlot( nbins);
393 pl->AddSamplingDistribution(sbDist);
394 return pl;
395 }
396 }
397 if (type == 2) {
399 if (bDist) {
400 pl = new SamplingDistPlot( nbins);
401 pl->AddSamplingDistribution(bDist);
402 return pl;
403 }
404 }
405 return 0;
406}
static const double x2[5]
static const double x1[5]
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
@ kRed
Definition: Rtypes.h:64
@ kGreen
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
@ kYellow
Definition: Rtypes.h:64
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
float * q
Definition: THbookFile.cxx:87
#define gROOT
Definition: TROOT.h:415
#define gPad
Definition: TVirtualPad.h:286
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:71
RooAbsArg * first() const
Class to plot an HypoTestInverterResult, result of the HypoTestInverter calculator.
void Draw(Option_t *opt="")
Draw the scan result in the current canvas Possible options: "" (default): draw observed + expected w...
HypoTestInverterResult * fResults
SamplingDistPlot * MakeTestStatPlot(int index, int type=0, int nbins=100)
Plot the test statistic distributions.
HypoTestInverterPlot(HypoTestInverterResult *results)
constructor from a HypoTestInverterResult class name and title are taken from the result class
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...
TGraphErrors * MakePlot(Option_t *opt="")
return a TGraphErrors with the obtained observed p-values resultinf from the scan By default (Option ...
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
double GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
get the signal and background test statistic distribution
int ArraySize() const
number of entries in the results array
std::vector< double > fXValues
number of points used to build expected p-values
double CLsplusbError(int index) const
return the observed CLsplusb value for the i-th entry
double CLsError(int index) const
return the observed CLb value for the i-th entry
double CLbError(int index) const
return the observed CLb value for the i-th entry
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results
double CLs(int index) const
return the observed CLb value for the i-th entry
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...
double CLb(int index) const
return the observed CLb value for the i-th entry
SamplingDistribution * GetExpectedPValueDist(int index) const
return expected distribution of p-values (Cls or Clsplusb)
double CLsplusb(int index) const
return the observed CLsplusb value for the i-th entry
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
SamplingDistribution * GetAltTestStatDist(int index) const
SamplingDistribution * GetBackgroundTestStatDist(int index) const
get the background test statistic distribution
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
HypoTestResult is a base class for results from hypothesis tests.
This class provides simple and straightforward utilities to plot SamplingDistribution objects.
Double_t AddSamplingDistribution(const SamplingDistribution *samplingDist, Option_t *drawOptions="NORMALIZE HIST")
adds the sampling distribution and returns the scale factor
This class simply holds a sampling distribution of some test statistic.
virtual Double_t ConfidenceLevel() const
return confidence level
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:42
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
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
Double_t GetXmax() const
Definition: TAxis.h:134
Double_t GetXmin() const
Definition: TAxis.h:133
virtual Int_t GetSize() const
Return the capacity of the collection, i.e.
Definition: TCollection.h:182
TGraph with asymmetric error bars.
virtual void SetPointEYlow(Int_t i, Double_t eyl)
Set EYlow for point i.
virtual void SetPointEYhigh(Int_t i, Double_t eyh)
Set EYhigh for point i.
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2257
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2312
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:753
TAxis * GetXaxis() const
Get x axis of the graph.
Definition: TGraph.cxx:1619
TAxis * GetYaxis() const
Get y axis of the graph.
Definition: TGraph.cxx:1629
TH1F * GetHistogram() const
Returns a pointer to the histogram used to draw the axis Takes into account the two following cases.
Definition: TGraph.cxx:1475
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6333
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
A simple line.
Definition: TLine.h:23
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:354
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:656
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:35
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:69
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.
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
Mother of all ROOT objects.
Definition: TObject.h:37
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
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
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TLine * line
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Double_t x[n]
Definition: legend1.C:17
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:150
Namespace for the RooStats classes.
Definition: Asimov.h:20
static constexpr double s
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Definition: TMathBase.h:337
Double_t Floor(Double_t x)
Definition: TMath.h:693
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:1190
Definition: graph.py:1
auto * l
Definition: textangle.C:4