ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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 /**
12  HypoTestInverterPlot class
13 **/
14 
15 // include header file of this class
17 
18 // include other header files
19 #include "RooStats/HybridResult.h"
21 #include "RooStats/HypoTestPlot.h"
23 
24 #include "TGraphErrors.h"
25 #include "TGraphAsymmErrors.h"
26 #include "TMultiGraph.h"
27 #include "TROOT.h"
28 #include "TLine.h"
29 #include "TAxis.h"
30 #include "TLegend.h"
31 #include "TH1.h"
32 #include "TPad.h"
33 #include "Math/DistFuncMathCore.h"
34 
35 #include <cmath>
36 
37 using namespace std;
38 
40 
41 using namespace RooStats;
42 
43 
45  TNamed( results->GetName(), results->GetTitle() ),
46  fResults(results)
47 {
48  // constructor from a HypoTestInverterResult class
49  // name and title are taken from the result class
50 }
51 
52 
53 HypoTestInverterPlot::HypoTestInverterPlot( const char* name,
54  const char* title,
55  HypoTestInverterResult* results ) :
56  TNamed( TString(name), TString(title) ),
57  fResults(results)
58 {
59  // constructor with name and title from a HypoTestInverterResult class
60 }
61 
62 
64 {
65  // Make the plot of the result of the scan
66  // using the observed data
67  // By default plot CLs or CLsb depending if the flag UseCLs is set
68  //
69  // If Option = "CLb" return CLb plot
70  // = "CLs+b" return CLs+b plot independently of the flag
71  // = "CLs" return CLs plot independently of the flag
72 
73  TString option(opt);
74  option.ToUpper();
75  int type = 0; // use default
76  if (option.Contains("CLB")) type = 1; // CLb
77  else if (option.Contains("CLS+B") || option.Contains("CLSPLUSB")) type = 2; // CLs+b
78  else if (option.Contains("CLS" )) type = 3; // CLs
79 
80  const int nEntries = fResults->ArraySize();
81 
82  // sort the arrays based on the x values
83  std::vector<unsigned int> index(nEntries);
84  TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
85 
86  // copy result in sorted arrays
87  std::vector<Double_t> xArray(nEntries);
88  std::vector<Double_t> yArray(nEntries);
89  std::vector<Double_t> yErrArray(nEntries);
90  for (int i=0; i<nEntries; i++) {
91  xArray[i] = fResults->GetXValue(index[i]);
92  if (type == 0) {
93  yArray[i] = fResults->GetYValue(index[i]);
94  yErrArray[i] = fResults->GetYError(index[i]);
95  } else if (type == 1) {
96  yArray[i] = fResults->CLb(index[i]);
97  yErrArray[i] = fResults->CLbError(index[i]);
98  } else if (type == 2) {
99  yArray[i] = fResults->CLsplusb(index[i]);
100  yErrArray[i] = fResults->CLsplusbError(index[i]);
101  } else if (type == 3) {
102  yArray[i] = fResults->CLs(index[i]);
103  yErrArray[i] = fResults->CLsError(index[i]);
104  }
105  }
106 
107  TGraphErrors* graph = new TGraphErrors(nEntries,&xArray.front(),&yArray.front(),0,&yErrArray.front());
108  TString pValueName = "CLs";
109  if (type == 1 ) pValueName = "CLb";
110  if (type == 2 || (type == 0 && !fResults->fUseCLs) ) pValueName = "CLs+b";
111  TString name = pValueName + TString("_observed");
112  TString title = TString("Observed ") + pValueName;
113  graph->SetName(name);
114  graph->SetTitle(title);
115  graph->SetMarkerStyle(20);
116  graph->SetLineWidth(2);
117  return graph;
118 }
119 
121 {
122  // Make the expected plot and the bands
123  // nsig1 and nsig2 indicates the n-sigma value for the bands
124  // if nsig1 = 0 no band is drawn (only expected value)
125  // if nsig2 > nsig1 (default is nsig1=1 and nsig2=2) the second band is also drawn
126  // The first band is drawn in green while the second in yellow
127  // THe return result is a TMultiGraph object
128 
129 
130 
131  const int nEntries = fResults->ArraySize();
132  bool doFirstBand = (nsig1 > 0);
133  bool doSecondBand = (nsig2 > nsig1);
134 
135  nsig1 = std::abs(nsig1);
136  nsig2 = std::abs(nsig2);
137 
138  // sort the arrays based on the x values
139  std::vector<unsigned int> index(nEntries);
140  TMath::SortItr(fResults->fXValues.begin(), fResults->fXValues.end(), index.begin(), false);
141 
142  // create the graphs
143  TGraph * g0 = new TGraph;
144  TString pValueName = "CLs";
145  if (!fResults->fUseCLs) pValueName = "CLs+b";
146  g0->SetTitle(TString::Format("Expected %s - Median",pValueName.Data()) );
147  TGraphAsymmErrors * g1 = 0;
148  TGraphAsymmErrors * g2 = 0;
149  if (doFirstBand) {
150  g1 = new TGraphAsymmErrors;
151  if (nsig1 - int(nsig1) < 0.01)
152  g1->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig1)) );
153  else
154  g1->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig1) );
155  }
156  if (doSecondBand) {
157  g2 = new TGraphAsymmErrors;
158  if (nsig2 - int(nsig2) < 0.01)
159  g2->SetTitle(TString::Format("Expected %s #pm %d #sigma",pValueName.Data(),int(nsig2)) );
160  else
161  g2->SetTitle(TString::Format("Expected %s #pm %3.1f #sigma",pValueName.Data(),nsig2) );
162  }
163  double p[7];
164  double q[7];
165  p[0] = ROOT::Math::normal_cdf(-nsig2);
166  p[1] = ROOT::Math::normal_cdf(-nsig1);
167  p[2] = 0.5;
168  p[3] = ROOT::Math::normal_cdf(nsig1);
169  p[4] = ROOT::Math::normal_cdf(nsig2);
170 
171  bool resultIsAsymptotic = ( !fResults->GetNullTestStatDist(0) && !fResults->GetAltTestStatDist(0) );
172  int np = 0;
173  for (int j=0; j<nEntries; ++j) {
174  int i = index[j]; // i is the order index
176  if ( !s) continue;
177  const std::vector<double> & values = s->GetSamplingDistribution();
178  // special case for asymptotic results (cannot use TMath::quantile in that case)
179  if (resultIsAsymptotic) {
180  double maxSigma = fResults->fgAsymptoticMaxSigma;
181  double dsig = 2* maxSigma/ (values.size() -1) ;
182  int i0 = (int) TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5);
183  int i1 = (int) TMath::Floor ( (-nsig1 + maxSigma )/dsig + 0.5);
184  int i2 = (int) TMath::Floor ( ( maxSigma)/dsig + 0.5);
185  int i3 = (int) TMath::Floor ( ( nsig1 + maxSigma)/dsig + 0.5);
186  int i4 = (int) TMath::Floor ( ( nsig2 + maxSigma)/dsig + 0.5);
187  q[0] = values[i0];
188  q[1] = values[i1];
189  q[2] = values[i2];
190  q[3] = values[i3];
191  q[4] = values[i4];
192  }
193  else {
194  double * x = const_cast<double *>(&values[0]); // need to change TMath::Quantiles
195  TMath::Quantiles(values.size(), 5, x,q,p,false);
196  }
197 
198  g0->SetPoint(np, fResults->GetXValue(i), q[2]);
199  if (g1) {
200  g1->SetPoint(np, fResults->GetXValue(i), q[2]);
201  g1->SetPointEYlow(np, q[2] - q[1]); // -1 sigma errorr
202  g1->SetPointEYhigh(np, q[3] - q[2]);//+1 sigma error
203  }
204  if (g2) {
205  g2->SetPoint(np, fResults->GetXValue(i), q[2]);
206 
207  g2->SetPointEYlow(np, q[2]-q[0]); // -2 -- -1 sigma error
208  g2->SetPointEYhigh(np, q[4]-q[2]);
209  }
210  if (s) delete s;
211  np++;
212  }
213 
214 
215 
216  TString name = GetName() + TString("_expected");
217  TString title = TString("Expected ") + GetTitle();
218  TMultiGraph* graph = new TMultiGraph(name,title);
219 
220  // set the graphics options and add in multi graph
221  // orderof adding is drawing order
222  if (g2) {
223  g2->SetFillColor(kYellow);
224  graph->Add(g2,"3");
225  }
226  if (g1) {
227  g1->SetFillColor(kGreen);
228  graph->Add(g1,"3");
229  }
230  g0->SetLineStyle(2);
231  g0->SetLineWidth(2);
232  graph->Add(g0,"L");
233 
234  return graph;
235 }
236 
238 {
239  // destructor
240 }
241 
243  // Draw the result in the current canvas
244  // Possible options:
245  // SAME : draw in the current axis
246  // OBS : draw only the observed plot
247  // EXP : draw only the expected plot
248  //
249  // CLB : draw also the CLB
250  // 2CL : drow both clsplusb and cls
251  //
252  // default draw observed + expected with 1 and 2 sigma bands
253 
254  TString option(opt);
255  option.ToUpper();
256  bool drawAxis = !option.Contains("SAME");
257  bool drawObs = option.Contains("OBS") || !option.Contains("EXP");
258  bool drawExp = option.Contains("EXP") || !option.Contains("OBS");
259  bool drawCLb = option.Contains("CLB");
260  bool draw2CL = option.Contains("2CL");
261 
262  TGraphErrors * gobs = 0;
263  TGraph * gplot = 0;
264  if (drawObs) {
265  gobs = MakePlot();
266  // add object to top-level directory to avoid mem leak
267  if (gROOT) gROOT->Add(gobs);
268  if (drawAxis) {
269  gobs->Draw("APL");
270  gplot = gobs;
271  gplot->GetHistogram()->SetTitle( GetTitle() );
272  }
273  else gobs->Draw("PL");
274 
275  }
276  TMultiGraph * gexp = 0;
277  if (drawExp) {
278  gexp = MakeExpectedPlot();
279  // add object to current directory to avoid mem leak
280  if (gROOT) gROOT->Add(gexp);
281  if (drawAxis && !drawObs) {
282  gexp->Draw("A");
283  if (gexp->GetHistogram()) gexp->GetHistogram()->SetTitle( GetTitle() );
284  gplot = (TGraph*) gexp->GetListOfGraphs()->First();
285  }
286  else
287  gexp->Draw();
288 
289  }
290 
291  // draw also an horizontal line at the desired conf level
292  if (gplot) {
293  double alpha = 1.-fResults->ConfidenceLevel();
294  double x1 = gplot->GetXaxis()->GetXmin();
295  double x2 = gplot->GetXaxis()->GetXmax();
296  TLine * line = new TLine(x1, alpha, x2,alpha);
297  line->SetLineColor(kRed);
298  line->Draw();
299  // put axis labels
300  RooAbsArg * arg = fResults->fParameters.first();
301  if (arg) gplot->GetXaxis()->SetTitle(arg->GetName());
302  gplot->GetYaxis()->SetTitle("p value");
303  }
304 
305 
306  TGraph *gclb = 0;
307  if (drawCLb) {
308  gclb = MakePlot("CLb");
309  if (gROOT) gROOT->Add(gclb);
310  gclb->SetMarkerColor(kBlue+4);
311  gclb->Draw("PL");
312  // draw in red observed cls or clsb
313  if (gobs) gobs->SetMarkerColor(kRed);
314  }
315  TGraph * gclsb = 0;
316  TGraph * gcls = 0;
317  if (draw2CL) {
318  if (fResults->fUseCLs) {
319  gclsb = MakePlot("CLs+b");
320  if (gROOT) gROOT->Add(gclsb);
321  gclsb->SetMarkerColor(kBlue);
322  gclsb->Draw("PL");
323  gclsb->SetLineStyle(3);
324  }
325  else {
326  gcls = MakePlot("CLs");
327  if (gROOT) gROOT->Add(gcls);
328  gcls->SetMarkerColor(kBlue);
329  gcls->Draw("PL");
330  gcls->SetLineStyle(3);
331  }
332  }
333  // draw again observed values otherwise will be covered by the bands
334  if (gobs) {
335  gobs->Draw("PL");
336  }
337 
338 
339  double y0 = 0.6;
340  double verticalSize = (gexp || draw2CL || drawCLb ) ? 0.3 : 0.15;
341  double y1 = y0 + verticalSize;
342  TLegend * l = new TLegend(0.6,y0,0.9,y1);
343  if (gobs) l->AddEntry(gobs,"","PEL");
344  if (gclsb) l->AddEntry(gclsb,"","PEL");
345  if (gcls) l->AddEntry(gcls,"","PEL");
346  if (gclb) l->AddEntry(gclb,"","PEL");
347  if (gexp) {
348  // loop in reverse order (opposite to drawing one)
349  int ngraphs = gexp->GetListOfGraphs()->GetSize();
350  for (int i = ngraphs-1; i>=0; --i) {
351  TObject * obj = gexp->GetListOfGraphs()->At(i);
352  TString lopt = "F";
353  if (i == ngraphs-1) lopt = "L";
354  if (obj) l->AddEntry(obj,"",lopt);
355  }
356  }
357  l->Draw();
358  // redraw the axis
359  if (gPad) gPad->RedrawAxis();
360 
361 }
362 
364  // plot the test statistic distributions
365  // type =0 null and alt
366  // type = 1 only null (S+B)
367  // type = 2 only alt (B)
368  SamplingDistPlot * pl = 0;
369  if (type == 0) {
371  if (result)
372  pl = new HypoTestPlot(*result, nbins );
373  return pl;
374  }
375  if (type == 1) {
377  if (sbDist) {
378  pl = new SamplingDistPlot( nbins);
379  pl->AddSamplingDistribution(sbDist);
380  return pl;
381  }
382  }
383  if (type == 2) {
385  if (bDist) {
386  pl = new SamplingDistPlot( nbins);
387  pl->AddSamplingDistribution(bDist);
388  return pl;
389  }
390  }
391  return 0;
392 }
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
Double_t Floor(Double_t x)
Definition: TMath.h:473
virtual Double_t ConfidenceLevel() const
return confidence level
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:35
TLine * line
double CLbError(int index) const
return the observed CLb value for the i-th entry
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:71
tuple g2
Definition: multifit.py:39
const char Option_t
Definition: RtypesCore.h:62
int ArraySize() const
number of entries in the results array
Definition: Rtypes.h:61
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:373
std::vector< double > values
Definition: TwoHistoFit2D.C:32
HypoTestResult is a base class for results from hypothesis tests.
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:37
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1088
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:344
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
Basic string class.
Definition: TString.h:137
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:37
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:2153
int nbins[3]
Definition: Rtypes.h:61
Definition: Rtypes.h:61
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...
SamplingDistribution * GetBackgroundTestStatDist(int index) const
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:311
tuple g1
Definition: multifit.py:38
const char * Data() const
Definition: TString.h:349
HypoTestInverterResult * fResults
TGraph with assymetric error bars.
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:2321
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:33
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:1173
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.
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
tuple np
Definition: multifit.py:30
virtual void Draw(Option_t *chopt="")
Draw this multigraph with its current attributes.
double CLsplusb(int index) const
return the observed CLsplusb value for the i-th entry
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...
ClassImp(RooStats::HypoTestInverterPlot) using namespace RooStats
TGraphErrors * MakePlot(Option_t *opt="")
return a TGraphErrors with the obtained observed p-values resultinf from the scan By default (Option ...
SamplingDistribution * GetExpectedPValueDist(int index) const
return expected distribution of p-values (Cls or Clsplusb)
TPaveLabel title(3, 27.1, 15, 28.7,"ROOT Environment and Tools")
virtual void SetFillColor(Color_t fcolor)
Definition: TAttFill.h:50
TLine * l
Definition: textangle.C:4
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
tuple pl
Definition: first.py:10
SamplingDistribution * GetAltTestStatDist(int index) const
This class simply holds a sampling distribution of some test statistic.
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
virtual Int_t GetSize() const
Definition: TCollection.h:95
static const double x1[5]
double CLs(int index) const
return the observed CLb value for the i-th entry
HypoTestInverterResult class: holds the array of hypothesis test results and compute a confidence int...
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:280
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 ...
double CLsplusbError(int index) const
return the observed CLsplusb value for the i-th entry
TH1F * GetHistogram() const
Returns a pointer to the histogram used to draw the axis.
#define name(a, b)
Definition: linkTestLib0.cpp:5
This class provides simple and straightforward utilities to plot SamplingDistribution objects...
Mother of all ROOT objects.
Definition: TObject.h:58
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:557
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2127
void Draw(Option_t *opt="")
Draw the scan result in the current canvas Possible options: "" (default): draw observed + expected w...
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
#define gPad
Definition: TVirtualPad.h:288
double result[121]
Definition: Rtypes.h:61
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
float * q
Definition: THbookFile.cxx:87
double CLb(int index) const
return the observed CLb value for the i-th entry
virtual void SetTitle(const char *title="")
Change (i.e. set) the title of the TNamed.
Definition: TNamed.cxx:152
TObject * obj
virtual void SetPointEYhigh(Int_t i, Double_t eyh)
Set EYhigh for point i.
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
Definition: TMath.h:977
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
std::vector< double > fXValues
max sigma value used to scan asymptotic expected p values