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
44using namespace std;
45
47
48using namespace RooStats;
49
50////////////////////////////////////////////////////////////////////////////////
51/// Constructor
52
55 fTestStat_data(-999.),
56 fComputationsNulDoneFlag(false),
57 fComputationsAltDoneFlag(false),
58 fSumLargerValues(false)
59{
60 // HybridResult default constructor (with name )
61}
62
63////////////////////////////////////////////////////////////////////////////////
64/// Constructor
65
67 const std::vector<double>& testStat_sb_vals,
68 const std::vector<double>& testStat_b_vals,
69 bool sumLargerValues ) :
71 fTestStat_data(-999.),
72 fComputationsNulDoneFlag(false),
73 fComputationsAltDoneFlag(false),
74 fSumLargerValues(sumLargerValues)
75{
76 // HybridResult constructor (with name, title and vectors of S+B and B values)
77
78 int vector_size_sb = testStat_sb_vals.size();
79 assert(vector_size_sb>0);
80
81 int vector_size_b = testStat_b_vals.size();
82 assert(vector_size_b>0);
83
84 fTestStat_sb.reserve(vector_size_sb);
85 for (int i=0;i<vector_size_sb;++i)
86 fTestStat_sb.push_back(testStat_sb_vals[i]);
87
88 fTestStat_b.reserve(vector_size_b);
89 for (int i=0;i<vector_size_b;++i)
90 fTestStat_b.push_back(testStat_b_vals[i]);
91}
92
93
94////////////////////////////////////////////////////////////////////////////////
95/// HybridResult destructor
96
98{
99
100 fTestStat_sb.clear();
101 fTestStat_b.clear();
102}
103
104////////////////////////////////////////////////////////////////////////////////
105/// set the value of the test statistics on data
106
107void HybridResult::SetDataTestStatistics(double testStat_data_val)
108{
111 fTestStat_data = testStat_data_val;
112 return;
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// Returns \f$1 - CL_{b}\f$ : the B p-value
117
119{
120 if (fComputationsNulDoneFlag==false) {
121 int nToys = fTestStat_b.size();
122 if (nToys==0) {
123 std::cout << "Error: no toy data present. Returning -1.\n";
124 return -1;
125 }
126
127 double larger_than_measured=0;
128 if (fSumLargerValues) {
129 for (int iToy=0;iToy<nToys;++iToy)
130 if ( fTestStat_b[iToy] >= fTestStat_data ) ++larger_than_measured;
131 } else {
132 for (int iToy=0;iToy<nToys;++iToy)
133 if ( fTestStat_b[iToy] <= fTestStat_data ) ++larger_than_measured;
134 }
135
136 if (larger_than_measured==0) std::cout << "Warning: CLb = 0 ... maybe more toys are needed!\n";
137
139 fNullPValue = 1-larger_than_measured/nToys;
140 }
141
142 return fNullPValue;
143}
144
145////////////////////////////////////////////////////////////////////////////////
146/// Returns \f$CL_{s+b}\f$ : the S+B p-value
147
149{
150 if (fComputationsAltDoneFlag==false) {
151 int nToys = fTestStat_b.size();
152 if (nToys==0) {
153 std::cout << "Error: no toy data present. Returning -1.\n";
154 return -1;
155 }
156
157 double larger_than_measured=0;
158 if (fSumLargerValues) {
159 for (int iToy=0;iToy<nToys;++iToy)
160 if ( fTestStat_sb[iToy] >= fTestStat_data ) ++larger_than_measured;
161 } else {
162 for (int iToy=0;iToy<nToys;++iToy)
163 if ( fTestStat_sb[iToy] <= fTestStat_data ) ++larger_than_measured;
164 }
165
166 if (larger_than_measured==0) std::cout << "Warning: CLsb = 0 ... maybe more toys are needed!\n";
167
169 fAlternatePValue = larger_than_measured/nToys;
170 }
171
172 return fAlternatePValue;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Returns an estimate of the error on \f$CL_{b}\f$ assuming a binomial
177/// error on \f$CL_{b}\f$:
178/// \f[
179/// \sigma_{CL_{b}} = \sqrt{CL_{b} \left( 1 - CL_{b} \right) / n_{toys}}
180/// \f]
181
183{
184 unsigned const int n = fTestStat_b.size();
185 return TMath::Sqrt(CLb() * (1. - CLb()) / n);
186}
187
188////////////////////////////////////////////////////////////////////////////////
189/// Returns an estimate of the error on \f$CL_{s+b}\f$ assuming a binomial
190/// error on \f$CL_{s+b}\f$:
191/// \f[
192/// \sigma_{CL_{s+b}} = \sqrt{CL_{s+b} \left( 1 - CL_{s+b} \right) / n_{toys}}
193/// \f]
194
196{
197 unsigned const int n = fTestStat_sb.size();
198 return TMath::Sqrt(CLsplusb() * (1. - CLsplusb()) / n);
199}
200
201////////////////////////////////////////////////////////////////////////////////
202/// Returns an estimate of the error on \f$CL_{s}\f$ through combination
203/// of the errors on \f$CL_{b}\f$ and \f$CL_{s+b}\f$:
204/// \f[
205/// \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}
206/// \f]
207
209{
210 unsigned const int n_b = fTestStat_b.size();
211 unsigned const int n_sb = fTestStat_sb.size();
212
213 if (CLb() == 0 || CLsplusb() == 0)
214 return 0;
215
216 double cl_b_err = (1. - CLb()) / (n_b * CLb());
217 double cl_sb_err = (1. - CLsplusb()) / (n_sb * CLsplusb());
218
219 return CLs() * TMath::Sqrt(cl_b_err + cl_sb_err);
220}
221
222////////////////////////////////////////////////////////////////////////////////
223/// add additional toy-MC experiments to the current results
224/// use the data test statistics of the added object if none is already present
225/// (otherwise, ignore the new one)
226
228{
229
230 int other_size_sb = other->GetTestStat_sb().size();
231 for (int i=0;i<other_size_sb;++i)
232 fTestStat_sb.push_back(other->GetTestStat_sb()[i]);
233
234 int other_size_b = other->GetTestStat_b().size();
235 for (int i=0;i<other_size_b;++i)
236 fTestStat_b.push_back(other->GetTestStat_b()[i]);
237
238 // if no data is present use the other's HybridResult's data
239 if (fTestStat_data==-999.)
241
244
245 return;
246}
247
248////////////////////////////////////////////////////////////////////////////////
249/// prepare a plot showing a result and return a pointer to a HybridPlot object
250/// the needed arguments are: an object name, a title and the number of bins in the plot
251
252HybridPlot* HybridResult::GetPlot(const char* name,const char* title, int n_bins)
253{
254 // default plot name
255 TString plot_name;
256 if ( TString(name)=="" ) {
257 plot_name += GetName();
258 plot_name += "_plot";
259 } else plot_name = name;
260
261 // default plot title
262 TString plot_title;
263 if ( TString(title)=="" ) {
264 plot_title += GetTitle();
265 plot_title += "_plot (";
266 plot_title += fTestStat_b.size();
267 plot_title += " toys)";
268 } else plot_title = title;
269
270 HybridPlot* plot = new HybridPlot( plot_name.Data(),
271 plot_title.Data(),
275 n_bins,
276 true );
277 return plot;
278}
279
280////////////////////////////////////////////////////////////////////////////////
281/// Print out some information about the results
282
283void HybridResult::PrintMore(const char* /*options */)
284{
285 std::cout << "\nResults " << GetName() << ":\n"
286 << " - Number of S+B toys: " << fTestStat_b.size() << std::endl
287 << " - Number of B toys: " << fTestStat_sb.size() << std::endl
288 << " - test statistics evaluated on data: " << fTestStat_data << std::endl
289 << " - CL_b " << CLb() << std::endl
290 << " - CL_s+b " << CLsplusb() << std::endl
291 << " - CL_s " << CLs() << std::endl;
292
293 return;
294}
295
#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:378
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