Logo ROOT   6.14/05
Reference Guide
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 
19 Class encapsulating the result of the HybridCalculatorOriginal.
20 This class is a fresh rewrite in RooStats of
21 RooStatsCms/LimitResults developed by D. Piparo and G. Schott
22 New contributions to this class have been written by Matthias Wolf (error estimation)
23 
24 The objects of this class store and access with lightweight methods the
25 information calculated by LimitResults through a Lent calculation using
26 MC toy experiments.
27 In some ways can be considered an extended and extensible implementation of the
28 TConfidenceLevel 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 
40 #include "RooStats/HybridResult.h"
41 #include "RooStats/HybridPlot.h"
42 
43 
44 /// ClassImp for building the THtml documentation of the class
45 using namespace std;
46 
48 
49 using namespace RooStats;
50 
51 ////////////////////////////////////////////////////////////////////////////////
52 /// Construtor
53 
54 HybridResult::HybridResult( const char *name) :
55  HypoTestResult(name),
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 ) :
71  HypoTestResult(name,0,0),
72  fTestStat_data(-999.),
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 
108 void HybridResult::SetDataTestStatistics(double testStat_data_val)
109 {
110  fComputationsAltDoneFlag = false;
111  fComputationsNulDoneFlag = false;
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.)
241  fTestStat_data = other->GetTestStat_data();
242 
243  fComputationsAltDoneFlag = false;
244  fComputationsNulDoneFlag = false;
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 
253 HybridPlot* 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(),
273  fTestStat_sb,
274  fTestStat_b,
276  n_bins,
277  true );
278  return plot;
279 }
280 
281 ////////////////////////////////////////////////////////////////////////////////
282 /// Print out some information about the results
283 
284 void 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 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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:...
virtual Double_t CLb() const
Convert NullPValue into a "confidence level".
Double_t CLbError() const
The error on the "confidence level" of the null hypothesis.
void Add(HybridResult *other)
add additional toy-MC experiments to the current results use the data test statistics of the added ob...
HypoTestResult is a base class for results from hypothesis tests.
Double_t CLsError() const
The error on the ratio .
STL namespace.
virtual ~HybridResult()
Destructor of HybridResult.
void SetDataTestStatistics(double testStat_data_val)
set the value of the test statistics on data
This class provides the plots for the result of a study performed with the HybridCalculatorOriginal c...
Definition: HybridPlot.h:36
HybridResult(const char *name=0)
Default constructor.
Double_t NullPValue() const
Returns : the B p-value.
Double_t CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
virtual Double_t CLs() const
is simply (not a method, but a quantity)
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
double GetTestStat_data()
Get test statistics value for data.
Definition: HybridResult.h:57
std::vector< double > fTestStat_b
Definition: HybridResult.h:76
std::vector< double > GetTestStat_sb()
Get test statistics values for the sb model.
Definition: HybridResult.h:51
void PrintMore(const char *options)
Print out some information about the results.
std::vector< double > GetTestStat_b()
Get test statistics values for the b model.
Definition: HybridResult.h:54
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Class encapsulating the result of the HybridCalculatorOriginal.
Definition: HybridResult.h:25
const Int_t n
Definition: legend1.C:16
char name[80]
Definition: TGX11.cxx:109
Double_t AlternatePValue() const
Returns : the S+B p-value.
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
std::vector< double > fTestStat_sb
Definition: HybridResult.h:77