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 the
28TConfidenceLevel class (http://root.cern.ch/root/html/TConfidenceLevel.html).
29
30*/
31
32
33#include "RooDataHist.h"
34#include "RooDataSet.h"
35#include "RooGlobalFunc.h" // for RooFit::Extended()
36#include "RooNLLVar.h"
37#include "RooRealVar.h"
38#include "RooAbsData.h"
39
41#include "RooStats/HybridPlot.h"
42
43
44/// ClassImp for building the THtml documentation of the class
45using namespace std;
46
48
49using namespace RooStats;
50
51////////////////////////////////////////////////////////////////////////////////
52/// Construtor
53
56 fTestStat_data(-999.),
57 fComputationsNulDoneFlag(false),
58 fComputationsAltDoneFlag(false),
59 fSumLargerValues(false)
60{
61 // HybridResult default constructor (with name )
62}
63
64////////////////////////////////////////////////////////////////////////////////
65/// Construtor
66
68 const std::vector<double>& testStat_sb_vals,
69 const std::vector<double>& testStat_b_vals,
70 bool sumLargerValues ) :
72 fTestStat_data(-999.),
73 fComputationsNulDoneFlag(false),
74 fComputationsAltDoneFlag(false),
75 fSumLargerValues(sumLargerValues)
76{
77 // HybridResult constructor (with name, title and vectors of S+B and B values)
78
79 int vector_size_sb = testStat_sb_vals.size();
80 assert(vector_size_sb>0);
81
82 int vector_size_b = testStat_b_vals.size();
83 assert(vector_size_b>0);
84
85 fTestStat_sb.reserve(vector_size_sb);
86 for (int i=0;i<vector_size_sb;++i)
87 fTestStat_sb.push_back(testStat_sb_vals[i]);
88
89 fTestStat_b.reserve(vector_size_b);
90 for (int i=0;i<vector_size_b;++i)
91 fTestStat_b.push_back(testStat_b_vals[i]);
92}
93
94
95////////////////////////////////////////////////////////////////////////////////
96/// HybridResult destructor
97
99{
100
101 fTestStat_sb.clear();
102 fTestStat_b.clear();
103}
104
105////////////////////////////////////////////////////////////////////////////////
106/// set the value of the test statistics on data
107
108void HybridResult::SetDataTestStatistics(double testStat_data_val)
109{
112 fTestStat_data = testStat_data_val;
113 return;
114}
115
116////////////////////////////////////////////////////////////////////////////////
117/// Returns \f$1 - CL_{b}\f$ : the B p-value
118
120{
121 if (fComputationsNulDoneFlag==false) {
122 int nToys = fTestStat_b.size();
123 if (nToys==0) {
124 std::cout << "Error: no toy data present. Returning -1.\n";
125 return -1;
126 }
127
128 double larger_than_measured=0;
129 if (fSumLargerValues) {
130 for (int iToy=0;iToy<nToys;++iToy)
131 if ( fTestStat_b[iToy] >= fTestStat_data ) ++larger_than_measured;
132 } else {
133 for (int iToy=0;iToy<nToys;++iToy)
134 if ( fTestStat_b[iToy] <= fTestStat_data ) ++larger_than_measured;
135 }
136
137 if (larger_than_measured==0) std::cout << "Warning: CLb = 0 ... maybe more toys are needed!\n";
138
140 fNullPValue = 1-larger_than_measured/nToys;
141 }
142
143 return fNullPValue;
144}
145
146////////////////////////////////////////////////////////////////////////////////
147/// Returns \f$CL_{s+b}\f$ : the S+B p-value
148
150{
151 if (fComputationsAltDoneFlag==false) {
152 int nToys = fTestStat_b.size();
153 if (nToys==0) {
154 std::cout << "Error: no toy data present. Returning -1.\n";
155 return -1;
156 }
157
158 double larger_than_measured=0;
159 if (fSumLargerValues) {
160 for (int iToy=0;iToy<nToys;++iToy)
161 if ( fTestStat_sb[iToy] >= fTestStat_data ) ++larger_than_measured;
162 } else {
163 for (int iToy=0;iToy<nToys;++iToy)
164 if ( fTestStat_sb[iToy] <= fTestStat_data ) ++larger_than_measured;
165 }
166
167 if (larger_than_measured==0) std::cout << "Warning: CLsb = 0 ... maybe more toys are needed!\n";
168
170 fAlternatePValue = larger_than_measured/nToys;
171 }
172
173 return fAlternatePValue;
174}
175
176////////////////////////////////////////////////////////////////////////////////
177/// Returns an estimate of the error on \f$CL_{b}\f$ assuming a binomial
178/// error on \f$CL_{b}\f$:
179/// \f[
180/// \sigma_{CL_{b}} = \sqrt{CL_{b} \left( 1 - CL_{b} \right) / n_{toys}}
181/// \f]
182
184{
185 unsigned const int n = fTestStat_b.size();
186 return TMath::Sqrt(CLb() * (1. - CLb()) / n);
187}
188
189////////////////////////////////////////////////////////////////////////////////
190/// Returns an estimate of the error on \f$CL_{s+b}\f$ assuming a binomial
191/// error on \f$CL_{s+b}\f$:
192/// \f[
193/// \sigma_{CL_{s+b}} = \sqrt{CL_{s+b} \left( 1 - CL_{s+b} \right) / n_{toys}}
194/// \f]
195
197{
198 unsigned const int n = fTestStat_sb.size();
199 return TMath::Sqrt(CLsplusb() * (1. - CLsplusb()) / n);
200}
201
202////////////////////////////////////////////////////////////////////////////////
203/// Returns an estimate of the error on \f$CL_{s}\f$ through combination
204/// of the errors on \f$CL_{b}\f$ and \f$CL_{s+b}\f$:
205/// \f[
206/// \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}
207/// \f]
208
210{
211 unsigned const int n_b = fTestStat_b.size();
212 unsigned const int n_sb = fTestStat_sb.size();
213
214 if (CLb() == 0 || CLsplusb() == 0)
215 return 0;
216
217 double cl_b_err = (1. - CLb()) / (n_b * CLb());
218 double cl_sb_err = (1. - CLsplusb()) / (n_sb * CLsplusb());
219
220 return CLs() * TMath::Sqrt(cl_b_err + cl_sb_err);
221}
222
223////////////////////////////////////////////////////////////////////////////////
224/// add additional toy-MC experiments to the current results
225/// use the data test statistics of the added object if none is already present
226/// (otherwise, ignore the new one)
227
229{
230
231 int other_size_sb = other->GetTestStat_sb().size();
232 for (int i=0;i<other_size_sb;++i)
233 fTestStat_sb.push_back(other->GetTestStat_sb()[i]);
234
235 int other_size_b = other->GetTestStat_b().size();
236 for (int i=0;i<other_size_b;++i)
237 fTestStat_b.push_back(other->GetTestStat_b()[i]);
238
239 // if no data is present use the other's HybridResult's data
240 if (fTestStat_data==-999.)
242
245
246 return;
247}
248
249////////////////////////////////////////////////////////////////////////////////
250/// prepare a plot showing a result and return a pointer to a HybridPlot object
251/// the needed arguments are: an object name, a title and the number of bins in the plot
252
253HybridPlot* HybridResult::GetPlot(const char* name,const char* title, int n_bins)
254{
255 // default plot name
256 TString plot_name;
257 if ( TString(name)=="" ) {
258 plot_name += GetName();
259 plot_name += "_plot";
260 } else plot_name = name;
261
262 // default plot title
263 TString plot_title;
264 if ( TString(title)=="" ) {
265 plot_title += GetTitle();
266 plot_title += "_plot (";
267 plot_title += fTestStat_b.size();
268 plot_title += " toys)";
269 } else plot_title = title;
270
271 HybridPlot* plot = new HybridPlot( plot_name.Data(),
272 plot_title.Data(),
276 n_bins,
277 true );
278 return plot;
279}
280
281////////////////////////////////////////////////////////////////////////////////
282/// Print out some information about the results
283
284void HybridResult::PrintMore(const char* /*options */)
285{
286 std::cout << "\nResults " << GetName() << ":\n"
287 << " - Number of S+B toys: " << fTestStat_b.size() << std::endl
288 << " - Number of B toys: " << fTestStat_sb.size() << std::endl
289 << " - test statistics evaluated on data: " << fTestStat_data << std::endl
290 << " - CL_b " << CLb() << std::endl
291 << " - CL_s+b " << CLsplusb() << std::endl
292 << " - CL_s " << CLs() << std::endl;
293
294 return;
295}
296
#define ClassImp(name)
Definition Rtypes.h:364
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
Double_t CLsError() const
The error on the ratio .
double GetTestStat_data()
Get test statistics value for data.
Double_t NullPValue() const
Returns : the B p-value.
Double_t CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
std::vector< double > GetTestStat_b()
Get test statistics values for the b model.
virtual ~HybridResult()
Destructor of HybridResult.
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...
Double_t CLbError() const
The error on the "confidence level" of the null hypothesis.
HybridResult(const char *name=0)
Default constructor.
Double_t AlternatePValue() const
Returns : the S+B p-value.
std::vector< double > fTestStat_sb
void SetDataTestStatistics(double testStat_data_val)
set the value of the test statistics on data
HypoTestResult is a base class for results from hypothesis tests.
virtual Double_t CLb() const
Convert NullPValue into a "confidence level".
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
virtual Double_t CLs() const
is simply (not a method, but a quantity)
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
Basic string class.
Definition TString.h:136
const char * Data() const
Definition TString.h:369
const Int_t n
Definition legend1.C:16
Namespace for the RooStats classes.
Definition Asimov.h:19
Double_t Sqrt(Double_t x)
Definition TMath.h:641