Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
HybridResult.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2
3/*************************************************************************
4 * Project: RooStats *
5 * Package: RooFit/RooStats *
6 * Authors: *
7 * Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke *
8 *************************************************************************
9 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
10 * All rights reserved. *
11 * *
12 * For the licensing terms see $ROOTSYS/LICENSE. *
13 * For the list of contributors see $ROOTSYS/README/CREDITS. *
14 *************************************************************************/
15
16/** \class RooStats::HybridResult
17 \ingroup Roostats
18
19Class encapsulating the result of the HybridCalculatorOriginal.
20This class is a fresh rewrite in RooStats of
21RooStatsCms/LimitResults developed by D. Piparo and G. Schott
22New contributions to this class have been written by Matthias Wolf (error estimation)
23
24The objects of this class store and access with lightweight methods the
25information calculated by LimitResults through a Lent calculation using
26MC toy experiments.
27In some ways can be considered an extended and extensible implementation of
28TConfidenceLevel.
29*/
30
31
32#include "RooDataHist.h"
33#include "RooDataSet.h"
34#include "RooGlobalFunc.h" // for RooFit::Extended()
35#include "RooRealVar.h"
36#include "RooAbsData.h"
37
39#include "RooStats/HybridPlot.h"
40
41#include <TMath.h>
42
43/// ClassImp for building the THtml documentation of the class
44
46
47using namespace RooStats;
48
49////////////////////////////////////////////////////////////////////////////////
50/// Constructor
51
54 fTestStat_data(-999.),
55 fComputationsNulDoneFlag(false),
56 fComputationsAltDoneFlag(false),
57 fSumLargerValues(false)
58{
59 // HybridResult default constructor (with name )
60}
61
62////////////////////////////////////////////////////////////////////////////////
63/// Constructor
64
66 const std::vector<double>& testStat_sb_vals,
67 const std::vector<double>& testStat_b_vals,
68 bool sumLargerValues ) :
70 fTestStat_data(-999.),
71 fComputationsNulDoneFlag(false),
72 fComputationsAltDoneFlag(false),
73 fSumLargerValues(sumLargerValues)
74{
75 // HybridResult constructor (with name, title and vectors of S+B and B values)
76
77 int vector_size_sb = testStat_sb_vals.size();
78 assert(vector_size_sb>0);
79
80 int vector_size_b = testStat_b_vals.size();
81 assert(vector_size_b>0);
82
83 fTestStat_sb.reserve(vector_size_sb);
84 for (int i=0;i<vector_size_sb;++i)
85 fTestStat_sb.push_back(testStat_sb_vals[i]);
86
87 fTestStat_b.reserve(vector_size_b);
88 for (int i=0;i<vector_size_b;++i)
89 fTestStat_b.push_back(testStat_b_vals[i]);
90}
91
92
93////////////////////////////////////////////////////////////////////////////////
94/// HybridResult destructor
95
97{
98
99 fTestStat_sb.clear();
100 fTestStat_b.clear();
101}
102
103////////////////////////////////////////////////////////////////////////////////
104/// set the value of the test statistics on data
105
106void HybridResult::SetDataTestStatistics(double testStat_data_val)
107{
110 fTestStat_data = testStat_data_val;
111 return;
112}
113
114////////////////////////////////////////////////////////////////////////////////
115/// Returns \f$1 - CL_{b}\f$ : the B p-value
116
118{
119 if (fComputationsNulDoneFlag==false) {
120 int nToys = fTestStat_b.size();
121 if (nToys==0) {
122 std::cout << "Error: no toy data present. Returning -1.\n";
123 return -1;
124 }
125
126 double larger_than_measured=0;
127 if (fSumLargerValues) {
128 for (int iToy=0;iToy<nToys;++iToy)
129 if ( fTestStat_b[iToy] >= fTestStat_data ) ++larger_than_measured;
130 } else {
131 for (int iToy=0;iToy<nToys;++iToy)
132 if ( fTestStat_b[iToy] <= fTestStat_data ) ++larger_than_measured;
133 }
134
135 if (larger_than_measured==0) std::cout << "Warning: CLb = 0 ... maybe more toys are needed!\n";
136
138 fNullPValue = 1-larger_than_measured/nToys;
139 }
140
141 return fNullPValue;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145/// Returns \f$CL_{s+b}\f$ : the S+B p-value
146
148{
149 if (fComputationsAltDoneFlag==false) {
150 int nToys = fTestStat_b.size();
151 if (nToys==0) {
152 std::cout << "Error: no toy data present. Returning -1.\n";
153 return -1;
154 }
155
156 double larger_than_measured=0;
157 if (fSumLargerValues) {
158 for (int iToy=0;iToy<nToys;++iToy)
159 if ( fTestStat_sb[iToy] >= fTestStat_data ) ++larger_than_measured;
160 } else {
161 for (int iToy=0;iToy<nToys;++iToy)
162 if ( fTestStat_sb[iToy] <= fTestStat_data ) ++larger_than_measured;
163 }
164
165 if (larger_than_measured==0) std::cout << "Warning: CLsb = 0 ... maybe more toys are needed!\n";
166
168 fAlternatePValue = larger_than_measured/nToys;
169 }
170
171 return fAlternatePValue;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Returns an estimate of the error on \f$CL_{b}\f$ assuming a binomial
176/// error on \f$CL_{b}\f$:
177/// \f[
178/// \sigma_{CL_{b}} = \sqrt{CL_{b} \left( 1 - CL_{b} \right) / n_{toys}}
179/// \f]
180
182{
183 unsigned const int n = fTestStat_b.size();
184 return TMath::Sqrt(CLb() * (1. - CLb()) / n);
185}
186
187////////////////////////////////////////////////////////////////////////////////
188/// Returns an estimate of the error on \f$CL_{s+b}\f$ assuming a binomial
189/// error on \f$CL_{s+b}\f$:
190/// \f[
191/// \sigma_{CL_{s+b}} = \sqrt{CL_{s+b} \left( 1 - CL_{s+b} \right) / n_{toys}}
192/// \f]
193
195{
196 unsigned const int n = fTestStat_sb.size();
197 return TMath::Sqrt(CLsplusb() * (1. - CLsplusb()) / n);
198}
199
200////////////////////////////////////////////////////////////////////////////////
201/// Returns an estimate of the error on \f$CL_{s}\f$ through combination
202/// of the errors on \f$CL_{b}\f$ and \f$CL_{s+b}\f$:
203/// \f[
204/// \sigma_{CL_s} = CL_s \sqrt{\left( \frac{\sigma_{CL_{s+b}}}{CL_{s+b}} \right)^2 + \left( \frac{\sigma_{CL_{b}}}{CL_{b}} \right)^2}
205/// \f]
206
208{
209 unsigned const int n_b = fTestStat_b.size();
210 unsigned const int n_sb = fTestStat_sb.size();
211
212 if (CLb() == 0 || CLsplusb() == 0)
213 return 0;
214
215 double cl_b_err = (1. - CLb()) / (n_b * CLb());
216 double cl_sb_err = (1. - CLsplusb()) / (n_sb * CLsplusb());
217
218 return CLs() * TMath::Sqrt(cl_b_err + cl_sb_err);
219}
220
221////////////////////////////////////////////////////////////////////////////////
222/// add additional toy-MC experiments to the current results
223/// use the data test statistics of the added object if none is already present
224/// (otherwise, ignore the new one)
225
227{
228
229 int other_size_sb = other->GetTestStat_sb().size();
230 for (int i=0;i<other_size_sb;++i)
231 fTestStat_sb.push_back(other->GetTestStat_sb()[i]);
232
233 int other_size_b = other->GetTestStat_b().size();
234 for (int i=0;i<other_size_b;++i)
235 fTestStat_b.push_back(other->GetTestStat_b()[i]);
236
237 // if no data is present use the other's HybridResult's data
238 if (fTestStat_data==-999.)
240
243
244 return;
245}
246
247////////////////////////////////////////////////////////////////////////////////
248/// prepare a plot showing a result and return a pointer to a HybridPlot object
249/// the needed arguments are: an object name, a title and the number of bins in the plot
250
251HybridPlot* HybridResult::GetPlot(const char* name,const char* title, int n_bins)
252{
253 // default plot name
254 TString plot_name;
255 if ( TString(name)=="" ) {
256 plot_name += GetName();
257 plot_name += "_plot";
258 } else plot_name = name;
259
260 // default plot title
261 TString plot_title;
262 if ( TString(title)=="" ) {
263 plot_title += GetTitle();
264 plot_title += "_plot (";
265 plot_title += fTestStat_b.size();
266 plot_title += " toys)";
267 } else plot_title = title;
268
269 HybridPlot* plot = new HybridPlot( plot_name.Data(),
270 plot_title.Data(),
274 n_bins,
275 true );
276 return plot;
277}
278
279////////////////////////////////////////////////////////////////////////////////
280/// Print out some information about the results
281
282void HybridResult::PrintMore(const char* /*options */)
283{
284 std::cout << "\nResults " << GetName() << ":\n"
285 << " - Number of S+B toys: " << fTestStat_b.size() << std::endl
286 << " - Number of B toys: " << fTestStat_sb.size() << std::endl
287 << " - test statistics evaluated on data: " << fTestStat_data << std::endl
288 << " - CL_b " << CLb() << std::endl
289 << " - CL_s+b " << CLsplusb() << std::endl
290 << " - CL_s " << CLs() << std::endl;
291
292 return;
293}
294
#define ClassImp(name)
Definition Rtypes.h:377
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
char name[80]
Definition TGX11.cxx:110
This class provides the plots for the result of a study performed with the HybridCalculatorOriginal c...
Definition HybridPlot.h:36
Class encapsulating the result of the HybridCalculatorOriginal.
std::vector< double > GetTestStat_sb()
Get test statistics values for the sb model.
std::vector< double > fTestStat_b
~HybridResult() override
Destructor of HybridResult.
double AlternatePValue() const override
Returns : the S+B p-value.
double GetTestStat_data()
Get test statistics value for data.
std::vector< double > GetTestStat_b()
Get test statistics values for the b model.
void PrintMore(const char *options)
Print out some information about the results.
HybridPlot * GetPlot(const char *name, const char *title, int n_bins)
prepare a plot showing a result and return a pointer to a HybridPlot object the needed arguments are:...
void Add(HybridResult *other)
add additional toy-MC experiments to the current results use the data test statistics of the added ob...
std::vector< double > fTestStat_sb
double CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
double CLbError() const
The error on the "confidence level" of the null hypothesis.
void SetDataTestStatistics(double testStat_data_val)
set the value of the test statistics on data
HybridResult(const char *name=nullptr)
Default constructor.
double CLsError() const
The error on the ratio .
double NullPValue() const override
Returns : the B p-value.
HypoTestResult is a base class for results from hypothesis tests.
double fNullPValue
p-value for the null hypothesis (small number means disfavoured)
virtual double CLsplusb() const
Convert AlternatePValue into a "confidence level".
virtual double CLs() const
is simply (not a method, but a quantity)
double fAlternatePValue
p-value for the alternate hypothesis (small number means disfavoured)
virtual double CLb() const
Convert NullPValue into a "confidence level".
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:48
Basic string class.
Definition TString.h:139
const char * Data() const
Definition TString.h:376
const Int_t n
Definition legend1.C:16
Namespace for the RooStats classes.
Definition Asimov.h:19
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662