Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 a HypoTestInverterResult, the output of the HypoTestInverter calculator.
15
16It can be used to plot the obtained p-values ( CLb, CLs+b or CLs) for each scanned point, as well as
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 "TVirtualPad.h"
39#include "TError.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
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 using the observed data.
73/// By default plot CLs or CLsb depending if the flag UseCLs is set for the results
74/// that are passed to this instance.
75///
76/// \param opt Options according to following list:
77/// - Empty: Return CLs or CLs+b depending on the value of UseCLs.ƒ
78/// - "CLB": return CLb plot
79/// - "CLS+B" / "CLSPLUSB": return CLs+b plot independently of the flag
80/// - "CLS": return CLs plot independently of the flag
81
83{
84 TString option(opt);
85 option.ToUpper();
86 enum Mode_t { Default, CLb, CLsPlusb, CLs };
87 Mode_t type = Default;
88 if (option.Contains("CLB")) type = CLb;
89 else if (option.Contains("CLS+B") || option.Contains("CLSPLUSB")) type = CLsPlusb;
90 else if (option.Contains("CLS" )) type = CLs;
91
92 const int nEntries = fResults->ArraySize();
93
94 // sort the arrays based on the x values
95 std::vector<unsigned int> index(nEntries);
96 TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
97
98 // copy result in sorted arrays
99 std::vector<Double_t> xArray;
100 std::vector<Double_t> yArray;
101 std::vector<Double_t> yErrArray;
102
103 for (int i=0; i<nEntries; i++) {
104 double CLVal = 0., CLErr = 0.;
105 if (type == Default) {
106 CLVal = fResults->GetYValue(index[i]);
107 CLErr = fResults->GetYError(index[i]);
108 } else if (type == CLb) {
109 CLVal = fResults->CLb(index[i]);
110 CLErr = fResults->CLbError(index[i]);
111 } else if (type == CLsPlusb) {
112 CLVal = fResults->CLsplusb(index[i]);
113 CLErr = fResults->CLsplusbError(index[i]);
114 } else if (type == CLs) {
115 CLVal = fResults->CLs(index[i]);
116 CLErr = fResults->CLsError(index[i]);
117 }
118
119 if (CLVal < 0. || !std::isfinite(CLVal)) {
120 Warning("HypoTestInverterPlot::MakePlot", "Got a confidence level of %f at x=%f (failed fit?). Skipping this point.", CLVal, fResults->GetXValue(index[i]));
121 continue;
122 }
123
124 yArray.push_back(CLVal);
125 yErrArray.push_back(CLErr);
126 xArray.push_back(fResults->GetXValue(index[i]));
127 }
128
129 TGraphErrors* graph = new TGraphErrors(static_cast<Int_t>(xArray.size()),&xArray.front(),&yArray.front(),0,&yErrArray.front());
130 TString pValueName = "CLs";
131 if (type == CLb ) pValueName = "CLb";
132 if (type == CLsPlusb || (type == Default && !fResults->fUseCLs) ) pValueName = "CLs+b";
133 TString name = pValueName + TString("_observed");
134 TString title = TString("Observed ") + pValueName;
135 graph->SetName(name);
136 graph->SetTitle(title);
137 graph->SetMarkerStyle(20);
138 graph->SetLineWidth(2);
139 return graph;
140}
141
142////////////////////////////////////////////////////////////////////////////////
143/// Make the expected plot and the bands
144/// nsig1 and nsig2 indicates the n-sigma value for the bands
145/// if nsig1 = 0 no band is drawn (only expected value)
146/// if nsig2 > nsig1 (default is nsig1=1 and nsig2=2) the second band is also drawn
147/// The first band is drawn in green while the second in yellow
148/// THe return result is a TMultiGraph object
149
151{
152 const int nEntries = fResults->ArraySize();
153 bool doFirstBand = (nsig1 > 0);
154 bool doSecondBand = (nsig2 > nsig1);
155
156 nsig1 = std::abs(nsig1);
157 nsig2 = std::abs(nsig2);
158
159 // sort the arrays based on the x values
160 std::vector<unsigned int> index(nEntries);
161 TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
162
163 // create the graphs
164 TGraph * g0 = new TGraph;
165 TString pValueName = "CLs";
166 if (!fResults->fUseCLs) pValueName = "CLs+b";
167 g0->SetTitle(TString::Format("Expected %s - Median",pValueName.Data()) );
168 TGraphAsymmErrors * g1 = 0;
169 TGraphAsymmErrors * g2 = 0;
170 if (doFirstBand) {
171 g1 = new TGraphAsymmErrors;
172 if (nsig1 - int(nsig1) < 0.01)
173 g1->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig1)) );
174 else
175 g1->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig1) );
176 }
177 if (doSecondBand) {
178 g2 = new TGraphAsymmErrors;
179 if (nsig2 - int(nsig2) < 0.01)
180 g2->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig2)) );
181 else
182 g2->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig2) );
183 }
184 double p[7];
185 double q[7];
186 p[0] = ROOT::Math::normal_cdf(-nsig2);
187 p[1] = ROOT::Math::normal_cdf(-nsig1);
188 p[2] = 0.5;
189 p[3] = ROOT::Math::normal_cdf(nsig1);
190 p[4] = ROOT::Math::normal_cdf(nsig2);
191
192 bool resultIsAsymptotic = ( !fResults->GetNullTestStatDist(0) && !fResults->GetAltTestStatDist(0) );
193 int np = 0;
194 for (int j=0; j<nEntries; ++j) {
195 int i = index[j]; // i is the order index
197 if ( !s) continue;
198 const std::vector<double> & values = s->GetSamplingDistribution();
199 // special case for asymptotic results (cannot use TMath::quantile in that case)
200 if (resultIsAsymptotic) {
201 double maxSigma = fResults->fgAsymptoticMaxSigma;
202 double dsig = 2* maxSigma/ (values.size() -1) ;
203 int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5);
204 int i1 = (int) TMath::Floor ( (-nsig1 + maxSigma )/dsig + 0.5);
205 int i2 = (int) TMath::Floor ( ( maxSigma)/dsig + 0.5);
206 int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma)/dsig + 0.5);
207 int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma)/dsig + 0.5);
208 q[0] = values[i0];
209 q[1] = values[i1];
210 q[2] = values[i2];
211 q[3] = values[i3];
212 q[4] = values[i4];
213 }
214 else {
215 double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
216 TMath::Quantiles(values.size(), 5, x,q,p,false);
217 }
218
219 g0->SetPoint(np, fResults->GetXValue(i), q[2]);
220 if (g1) {
221 g1->SetPoint(np, fResults->GetXValue(i), q[2]);
222 g1->SetPointEYlow(np, q[2] - q[1]); // -1 sigma errorr
223 g1->SetPointEYhigh(np, q[3] - q[2]);//+1 sigma error
224 }
225 if (g2) {
226 g2->SetPoint(np, fResults->GetXValue(i), q[2]);
227
228 g2->SetPointEYlow(np, q[2]-q[0]); // -2 -- -1 sigma error
229 g2->SetPointEYhigh(np, q[4]-q[2]);
230 }
231 if (s) delete s;
232 np++;
233 }
234
235 TString name = GetName() + TString("_expected");
236 TString title = TString("Expected ") + GetTitle();
237 TMultiGraph* graph = new TMultiGraph(name,title);
238
239 // set the graphics options and add in multi graph
240 // orderof adding is drawing order
241 if (g2) {
243 graph->Add(g2,"3");
244 }
245 if (g1) {
246 g1->SetFillColor(kGreen);
247 graph->Add(g1,"3");
248 }
249 g0->SetLineStyle(2);
250 g0->SetLineWidth(2);
251 graph->Add(g0,"L");
252
253 return graph;
254}
255
256////////////////////////////////////////////////////////////////////////////////
257/// destructor
258
260{
261}
262
263////////////////////////////////////////////////////////////////////////////////
264/// Draw the result in the current canvas
265/// Possible options:
266/// - SAME : draw in the current axis
267/// - OBS : draw only the observed plot
268/// - EXP : draw only the expected plot
269///
270/// - CLB : draw also the CLB
271/// - 2CL : draw both clsplusb and cls
272///
273/// default draw observed + expected with 1 and 2 sigma bands
274
276
277 TString option(opt);
278 option.ToUpper();
279 bool drawAxis = !option.Contains("SAME");
280 bool drawObs = option.Contains("OBS") || !option.Contains("EXP");
281 bool drawExp = option.Contains("EXP") || !option.Contains("OBS");
282 bool drawCLb = option.Contains("CLB");
283 bool draw2CL = option.Contains("2CL");
284
285 TGraphErrors * gobs = 0;
286 TGraph * gplot = 0;
287 if (drawObs) {
288 gobs = MakePlot();
289 // add object to top-level directory to avoid mem leak
290 if (gROOT) gROOT->Add(gobs);
291 if (drawAxis) {
292 gobs->Draw("APL");
293 gplot = gobs;
294 gplot->GetHistogram()->SetTitle( GetTitle() );
295 }
296 else gobs->Draw("PL");
297
298 }
299 TMultiGraph * gexp = 0;
300 if (drawExp) {
301 gexp = MakeExpectedPlot();
302 // add object to current directory to avoid mem leak
303 if (gROOT) gROOT->Add(gexp);
304 if (drawAxis && !drawObs) {
305 gexp->Draw("A");
306 if (gexp->GetHistogram()) gexp->GetHistogram()->SetTitle( GetTitle() );
307 gplot = (TGraph*) gexp->GetListOfGraphs()->First();
308 }
309 else
310 gexp->Draw();
311
312 }
313
314 // draw also an horizontal line at the desired conf level
315 if (gplot) {
316 double alpha = 1.-fResults->ConfidenceLevel();
317 double x1 = gplot->GetXaxis()->GetXmin();
318 double x2 = gplot->GetXaxis()->GetXmax();
319 TLine * line = new TLine(x1, alpha, x2,alpha);
321 line->Draw();
322 // put axis labels
324 if (arg) gplot->GetXaxis()->SetTitle(arg->GetName());
325 gplot->GetYaxis()->SetTitle("p value");
326 }
327
328
329 TGraph *gclb = 0;
330 if (drawCLb) {
331 gclb = MakePlot("CLb");
332 if (gROOT) gROOT->Add(gclb);
333 gclb->SetMarkerColor(kBlue+4);
334 gclb->Draw("PL");
335 // draw in red observed cls or clsb
336 if (gobs) gobs->SetMarkerColor(kRed);
337 }
338 TGraph * gclsb = 0;
339 TGraph * gcls = 0;
340 if (draw2CL) {
341 if (fResults->fUseCLs) {
342 gclsb = MakePlot("CLs+b");
343 if (gROOT) gROOT->Add(gclsb);
344 gclsb->SetMarkerColor(kBlue);
345 gclsb->Draw("PL");
346 gclsb->SetLineStyle(3);
347 }
348 else {
349 gcls = MakePlot("CLs");
350 if (gROOT) gROOT->Add(gcls);
351 gcls->SetMarkerColor(kBlue);
352 gcls->Draw("PL");
353 gcls->SetLineStyle(3);
354 }
355 }
356 // draw again observed values otherwise will be covered by the bands
357 if (gobs) {
358 gobs->Draw("PL");
359 }
360
361
362 double y0 = 0.6;
363 double verticalSize = (gexp || draw2CL || drawCLb ) ? 0.3 : 0.15;
364 double y1 = y0 + verticalSize;
365 TLegend * l = new TLegend(0.6,y0,0.9,y1);
366 if (gobs) l->AddEntry(gobs,"","PEL");
367 if (gclsb) l->AddEntry(gclsb,"","PEL");
368 if (gcls) l->AddEntry(gcls,"","PEL");
369 if (gclb) l->AddEntry(gclb,"","PEL");
370 if (gexp) {
371 // loop in reverse order (opposite to drawing one)
372 int ngraphs = gexp->GetListOfGraphs()->GetSize();
373 for (int i = ngraphs-1; i>=0; --i) {
374 TObject * obj = gexp->GetListOfGraphs()->At(i);
375 TString lopt = "F";
376 if (i == ngraphs-1) lopt = "L";
377 if (obj) l->AddEntry(obj,"",lopt);
378 }
379 }
380 l->Draw();
381 // redraw the axis
382 if (gPad) gPad->RedrawAxis();
383
384}
385
386////////////////////////////////////////////////////////////////////////////////
387/// \param index Index of the result stored in HypoTestInverterResult
388/// \param type Type of the test (see below)
389/// \param nbins Number of bins
390/// - type =0 null and alt
391/// - type = 1 only null (S+B)
392/// - type = 2 only alt (B)
393
395 SamplingDistPlot * pl = 0;
396 if (type == 0) {
397 HypoTestResult * result = (HypoTestResult*) fResults->fYObjects.At(index);
398 if (result)
399 pl = new HypoTestPlot(*result, nbins );
400 return pl;
401 }
402 if (type == 1) {
404 if (sbDist) {
405 pl = new SamplingDistPlot( nbins);
406 pl->AddSamplingDistribution(sbDist);
407 return pl;
408 }
409 }
410 if (type == 2) {
412 if (bDist) {
413 pl = new SamplingDistPlot( nbins);
414 pl->AddSamplingDistribution(bDist);
415 return pl;
416 }
417 }
418 return 0;
419}
static const double x2[5]
static const double x1[5]
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
@ kRed
Definition Rtypes.h:66
@ kGreen
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
@ kYellow
Definition Rtypes.h:66
char name[80]
Definition TGX11.cxx:110
int type
Definition TGX11.cxx:121
float * q
#define gROOT
Definition TROOT.h:406
#define gPad
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:72
RooAbsArg * first() const
Class to plot a HypoTestInverterResult, the output 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...
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.
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
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.
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.
A TGraph is an 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:2284
virtual void SetTitle(const char *title="")
Change (i.e.
Definition TGraph.cxx:2339
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition TGraph.cxx:769
TAxis * GetXaxis() const
Get x axis of the graph.
Definition TGraph.cxx:1640
TAxis * GetYaxis() const
Get y axis of the graph.
Definition TGraph.cxx:1650
TH1F * GetHistogram() const
Returns a pointer to the histogram used to draw the axis Takes into account the two following cases.
Definition TGraph.cxx:1491
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition TH1.cxx:6678
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
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:357
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition TList.cxx:659
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition TMultiGraph.h:36
TList * GetListOfGraphs() const
Definition TMultiGraph.h:70
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 Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:879
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:197
Basic string class.
Definition TString.h:136
const char * Data() const
Definition TString.h:369
void ToUpper()
Change string to upper case.
Definition TString.cxx:1158
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:2331
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
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
Namespace for the RooStats classes.
Definition Asimov.h:19
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:703
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
Definition graph.py:1
auto * l
Definition textangle.C:4