Logo ROOT   6.08/07
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 ////////////////////////////////////////////////////////////////////////////////
17 
18 
19 
20 #include "RooDataHist.h"
21 #include "RooDataSet.h"
22 #include "RooGlobalFunc.h" // for RooFit::Extended()
23 #include "RooNLLVar.h"
24 #include "RooRealVar.h"
25 #include "RooAbsData.h"
26 
27 #include "RooStats/HybridResult.h"
28 #include "RooStats/HybridPlot.h"
29 
30 
31 /// ClassImp for building the THtml documentation of the class
32 using namespace std;
33 
35 
36 using namespace RooStats;
37 
38 ///////////////////////////////////////////////////////////////////////////
39 
40 HybridResult::HybridResult( const char *name) :
41  HypoTestResult(name),
42  fTestStat_data(-999.),
43  fComputationsNulDoneFlag(false),
44  fComputationsAltDoneFlag(false),
45  fSumLargerValues(false)
46 {
47  // HybridResult default constructor (with name )
48 }
49 
50 ///////////////////////////////////////////////////////////////////////////
51 
53  const std::vector<double>& testStat_sb_vals,
54  const std::vector<double>& testStat_b_vals,
55  bool sumLargerValues ) :
56  HypoTestResult(name,0,0),
57  fTestStat_data(-999.),
60  fSumLargerValues(sumLargerValues)
61 {
62  // HybridResult constructor (with name, title and vectors of S+B and B values)
63 
64  int vector_size_sb = testStat_sb_vals.size();
65  assert(vector_size_sb>0);
66 
67  int vector_size_b = testStat_b_vals.size();
68  assert(vector_size_b>0);
69 
70  fTestStat_sb.reserve(vector_size_sb);
71  for (int i=0;i<vector_size_sb;++i)
72  fTestStat_sb.push_back(testStat_sb_vals[i]);
73 
74  fTestStat_b.reserve(vector_size_b);
75  for (int i=0;i<vector_size_b;++i)
76  fTestStat_b.push_back(testStat_b_vals[i]);
77 }
78 
79 
80 ///////////////////////////////////////////////////////////////////////////
81 /// HybridResult destructor
82 
84 {
85 
86  fTestStat_sb.clear();
87  fTestStat_b.clear();
88 }
89 
90 ///////////////////////////////////////////////////////////////////////////
91 /// set the value of the test statistics on data
92 
93 void HybridResult::SetDataTestStatistics(double testStat_data_val)
94 {
97  fTestStat_data = testStat_data_val;
98  return;
99 }
100 
101 ///////////////////////////////////////////////////////////////////////////
102 /// Returns \f$1 - CL_{b}\f$ : the B p-value
103 
105 {
106  if (fComputationsNulDoneFlag==false) {
107  int nToys = fTestStat_b.size();
108  if (nToys==0) {
109  std::cout << "Error: no toy data present. Returning -1.\n";
110  return -1;
111  }
112 
113  double larger_than_measured=0;
114  if (fSumLargerValues) {
115  for (int iToy=0;iToy<nToys;++iToy)
116  if ( fTestStat_b[iToy] >= fTestStat_data ) ++larger_than_measured;
117  } else {
118  for (int iToy=0;iToy<nToys;++iToy)
119  if ( fTestStat_b[iToy] <= fTestStat_data ) ++larger_than_measured;
120  }
121 
122  if (larger_than_measured==0) std::cout << "Warning: CLb = 0 ... maybe more toys are needed!\n";
123 
125  fNullPValue = 1-larger_than_measured/nToys;
126  }
127 
128  return fNullPValue;
129 }
130 
131 ///////////////////////////////////////////////////////////////////////////
132 /// Returns \f$CL_{s+b}\f$ : the S+B p-value
133 
135 {
136  if (fComputationsAltDoneFlag==false) {
137  int nToys = fTestStat_b.size();
138  if (nToys==0) {
139  std::cout << "Error: no toy data present. Returning -1.\n";
140  return -1;
141  }
142 
143  double larger_than_measured=0;
144  if (fSumLargerValues) {
145  for (int iToy=0;iToy<nToys;++iToy)
146  if ( fTestStat_sb[iToy] >= fTestStat_data ) ++larger_than_measured;
147  } else {
148  for (int iToy=0;iToy<nToys;++iToy)
149  if ( fTestStat_sb[iToy] <= fTestStat_data ) ++larger_than_measured;
150  }
151 
152  if (larger_than_measured==0) std::cout << "Warning: CLsb = 0 ... maybe more toys are needed!\n";
153 
155  fAlternatePValue = larger_than_measured/nToys;
156  }
157 
158  return fAlternatePValue;
159 }
160 
161 ///////////////////////////////////////////////////////////////////////////
162 /// Returns an estimate of the error on \f$CL_{b}\f$ assuming a binomial
163 /// error on \f$CL_{b}\f$:
164 /// \f[
165 /// \sigma_{CL_{b}} = \sqrt{CL_{b} \left( 1 - CL_{b} \right) / n_{toys}}
166 /// \f]
167 
169 {
170  unsigned const int n = fTestStat_b.size();
171  return TMath::Sqrt(CLb() * (1. - CLb()) / n);
172 }
173 
174 ///////////////////////////////////////////////////////////////////////////
175 /// Returns an estimate of the error on \f$CL_{s+b}\f$ assuming a binomial
176 /// error on \f$CL_{s+b}\f$:
177 /// \f[
178 /// \sigma_{CL_{s+b}} = \sqrt{CL_{s+b} \left( 1 - CL_{s+b} \right) / n_{toys}}
179 /// \f]
180 
182 {
183  unsigned const int n = fTestStat_sb.size();
184  return TMath::Sqrt(CLsplusb() * (1. - CLsplusb()) / n);
185 }
186 
187 ///////////////////////////////////////////////////////////////////////////
188 /// Returns an estimate of the error on \f$CL_{s}\f$ through combination
189 /// of the errors on \f$CL_{b}\f$ and \f$CL_{s+b}\f$:
190 /// \f[
191 /// \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}
192 /// \f]
193 
195 {
196  unsigned const int n_b = fTestStat_b.size();
197  unsigned const int n_sb = fTestStat_sb.size();
198 
199  if (CLb() == 0 || CLsplusb() == 0)
200  return 0;
201 
202  double cl_b_err = (1. - CLb()) / (n_b * CLb());
203  double cl_sb_err = (1. - CLsplusb()) / (n_sb * CLsplusb());
204 
205  return CLs() * TMath::Sqrt(cl_b_err + cl_sb_err);
206 }
207 
208 ///////////////////////////////////////////////////////////////////////////
209 /// add additional toy-MC experiments to the current results
210 /// use the data test statistics of the added object if none is already present
211 /// (otherwise, ignore the new one)
212 
214 {
215 
216  int other_size_sb = other->GetTestStat_sb().size();
217  for (int i=0;i<other_size_sb;++i)
218  fTestStat_sb.push_back(other->GetTestStat_sb()[i]);
219 
220  int other_size_b = other->GetTestStat_b().size();
221  for (int i=0;i<other_size_b;++i)
222  fTestStat_b.push_back(other->GetTestStat_b()[i]);
223 
224  // if no data is present use the other's HybridResult's data
225  if (fTestStat_data==-999.)
226  fTestStat_data = other->GetTestStat_data();
227 
228  fComputationsAltDoneFlag = false;
229  fComputationsNulDoneFlag = false;
230 
231  return;
232 }
233 
234 ///////////////////////////////////////////////////////////////////////////
235 /// prepare a plot showing a result and return a pointer to a HybridPlot object
236 /// the needed arguments are: an object name, a title and the number of bins in the plot
237 
238 HybridPlot* HybridResult::GetPlot(const char* name,const char* title, int n_bins)
239 {
240  // default plot name
241  TString plot_name;
242  if ( TString(name)=="" ) {
243  plot_name += GetName();
244  plot_name += "_plot";
245  } else plot_name = name;
246 
247  // default plot title
248  TString plot_title;
249  if ( TString(title)=="" ) {
250  plot_title += GetTitle();
251  plot_title += "_plot (";
252  plot_title += fTestStat_b.size();
253  plot_title += " toys)";
254  } else plot_title = title;
255 
256  HybridPlot* plot = new HybridPlot( plot_name.Data(),
257  plot_title.Data(),
258  fTestStat_sb,
259  fTestStat_b,
261  n_bins,
262  true );
263  return plot;
264 }
265 
266 ///////////////////////////////////////////////////////////////////////////
267 /// Print out some information about the results
268 
269 void HybridResult::PrintMore(const char* /*options */)
270 {
271  std::cout << "\nResults " << GetName() << ":\n"
272  << " - Number of S+B toys: " << fTestStat_b.size() << std::endl
273  << " - Number of B toys: " << fTestStat_sb.size() << std::endl
274  << " - test statistics evaluated on data: " << fTestStat_data << std::endl
275  << " - CL_b " << CLb() << std::endl
276  << " - CL_s+b " << CLsplusb() << std::endl
277  << " - CL_s " << CLs() << std::endl;
278 
279  return;
280 }
281 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
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:53
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:279
double Double_t
Definition: RtypesCore.h:55
double GetTestStat_data()
Get test statistics value for data.
Definition: HybridResult.h:74
std::vector< double > fTestStat_b
Definition: HybridResult.h:93
std::vector< double > GetTestStat_sb()
Get test statistics values for the sb model.
Definition: HybridResult.h:68
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:71
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
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:52
std::vector< double > fTestStat_sb
Definition: HybridResult.h:94