Logo ROOT  
Reference Guide
ConfidenceBelt.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11/*****************************************************************************
12 * Project: RooStats
13 * Package: RooFit/RooStats
14 * @(#)root/roofit/roostats:$Id$
15 * Original Author: Kyle Cranmer
16 * Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
17 *
18 *****************************************************************************/
19
20/** \class RooStats::ConfidenceBelt
21 \ingroup Roostats
22
23ConfidenceBelt is a concrete implementation of the ConfInterval interface.
24It implements simple general purpose interval of arbitrary dimensions and shape.
25It does not assume the interval is connected.
26It uses either a RooDataSet (eg. a list of parameter points in the interval) or
27a RooDataHist (eg. a Histogram-like object for small regions of the parameter space) to
28store the interval.
29
30*/
31
33
34#include "RooDataSet.h"
35#include "RooDataHist.h"
36
38
39#include <TMath.h>
40
42
43using namespace RooStats;
44using namespace std;
45
46////////////////////////////////////////////////////////////////////////////////
47/// Default constructor
48
49ConfidenceBelt::ConfidenceBelt() :
50 TNamed(), fParameterPoints(0)
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55/// Alternate constructor
56
58 TNamed(name,name), fParameterPoints(0)
59{
60}
61
62////////////////////////////////////////////////////////////////////////////////
63/// Alternate constructor
64
65ConfidenceBelt::ConfidenceBelt(const char* name, const char* title) :
66 TNamed(name,title), fParameterPoints(0)
67{
68}
69
70////////////////////////////////////////////////////////////////////////////////
71/// Alternate constructor
72
74 TNamed(name,name), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
75{
76}
77
78////////////////////////////////////////////////////////////////////////////////
79/// Alternate constructor
80
81ConfidenceBelt::ConfidenceBelt(const char* name, const char* title, RooAbsData& data) :
82 TNamed(name,title), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
83{
84}
85
86////////////////////////////////////////////////////////////////////////////////
87/// Destructor
88
90{
91}
92
93////////////////////////////////////////////////////////////////////////////////
94
95double ConfidenceBelt::GetAcceptanceRegionMin(RooArgSet& parameterPoint, double cl, double leftside) {
96 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
97 AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
98 return (region) ? region->GetLowerLimit() : TMath::QuietNaN();
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
103double ConfidenceBelt::GetAcceptanceRegionMax(RooArgSet& parameterPoint, double cl, double leftside) {
104 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
105 AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
106 return (region) ? region->GetUpperLimit() : TMath::QuietNaN();
107}
108
109////////////////////////////////////////////////////////////////////////////////
110
111vector<double> ConfidenceBelt::ConfidenceLevels() const {
112 vector<double> levels;
113 return levels;
114}
115
116////////////////////////////////////////////////////////////////////////////////
117
119 double lower, double upper,
120 double cl, double leftside){
121 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
122
123 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
124 RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
125
126 // cout << "add: " << tree << " " << hist << endl;
127
128 if( !this->CheckParameters(parameterPoint) )
129 std::cout << "problem with parameters" << std::endl;
130
131 Int_t luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
132 // cout << "lookup index = " << luIndex << endl;
133 if(luIndex <0 ) {
134 fSamplingSummaryLookup.Add(cl,leftside);
135 luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
136 cout << "lookup index = " << luIndex << endl;
137 }
138 AcceptanceRegion* thisRegion = new AcceptanceRegion(luIndex, lower, upper);
139
140 if( hist ) {
141 // need a way to get index for given point
142 // Can do this by setting hist's internal parameters to desired values
143 // need a better way
144 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
145 // int index = hist->calcTreeIndex(); // get index
146 int index = hist->getIndex(parameterPoint); // get index
147 // cout << "hist index = " << index << endl;
148
149 // allocate memory if necessary. numEntries is overkill?
150 if((Int_t)fSamplingSummaries.size() <= index) {
151 fSamplingSummaries.reserve( hist->numEntries() );
152 fSamplingSummaries.resize( hist->numEntries() );
153 }
154
155 // set the region for this point (check for duplicate?)
156 fSamplingSummaries.at(index) = *thisRegion;
157 }
158 else if( tree ){
159 // int index = tree->getIndex(parameterPoint);
160 int index = dsIndex;
161 // cout << "tree index = " << index << endl;
162
163 // allocate memory if necessary. numEntries is overkill?
164 if((Int_t)fSamplingSummaries.size() <= index){
165 fSamplingSummaries.reserve( tree->numEntries() );
166 fSamplingSummaries.resize( tree->numEntries() );
167 }
168
169 // set the region for this point (check for duplicate?)
170 fSamplingSummaries.at( index ) = *thisRegion;
171 }
172}
173
174////////////////////////////////////////////////////////////////////////////////
175
177 double cl, double leftside){
178 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
179
180 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
181 RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
182
183 if( !this->CheckParameters(parameterPoint) )
184 std::cout << "problem with parameters" << std::endl;
185
186
187 if( hist ) {
188 // need a way to get index for given point
189 // Can do this by setting hist's internal parameters to desired values
190 // need a better way
191 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
192 // int index = hist->calcTreeIndex(); // get index
193 int index = hist->getIndex(parameterPoint); // get index
194
195 // allocate memory if necessary. numEntries is overkill?
196 if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( hist->numEntries() );
197
198 // set the region for this point (check for duplicate?)
199 fSamplingSummaries.at(index) = region;
200 }
201 else if( tree ){
202 tree->add( parameterPoint ); // assume it's unique for now
203 int index = tree->numEntries() - 1; //check that last point added has index nEntries -1
204 // allocate memory if necessary. numEntries is overkill?
205 if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( tree->numEntries() );
206
207 // set the region for this point (check for duplicate?)
208 fSamplingSummaries.at( index ) = region;
209 }
210}
211
212////////////////////////////////////////////////////////////////////////////////
213/// Method to determine if a parameter point is in the interval
214
215AcceptanceRegion* ConfidenceBelt::GetAcceptanceRegion(RooArgSet &parameterPoint, double cl, double leftside)
216{
217 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
218
219 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
220 RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
221
222 if( !this->CheckParameters(parameterPoint) ){
223 std::cout << "problem with parameters" << std::endl;
224 return 0;
225 }
226
227 if( hist ) {
228 // need a way to get index for given point
229 // Can do this by setting hist's internal parameters to desired values
230 // need a better way
231 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
232 // Int_t index = hist->calcTreeIndex(); // get index
233 int index = hist->getIndex(parameterPoint); // get index
234 if (index >= (int)fSamplingSummaries.size())
235 throw std::runtime_error("ConfidenceBelt::GetAcceptanceRegion: Sampling summaries are not filled yet. Switch on NeymanConstruction::CreateConfBelt() or FeldmanCousins::CreateConfBelt().");
236
237 return &(fSamplingSummaries[index].GetAcceptanceRegion());
238 }
239 else if( tree ){
240 // need a way to get index for given point
241 // RooStats::SetParameters(&parameterPoint, tree->get()); // set tree's parameters to desired values
242 Int_t index = 0; //need something like tree->calcTreeIndex();
243 const RooArgSet* thisPoint = 0;
244 for(index=0; index<tree->numEntries(); ++index){
245 thisPoint = tree->get(index);
246 bool samePoint = true;
247 for (auto const *myarg : static_range_cast<RooRealVar *>(parameterPoint)) {
248 if (samePoint == false)
249 break;
250 if(myarg->getVal() != thisPoint->getRealValue(myarg->GetName()))
251 samePoint = false;
252 }
253 if(samePoint)
254 break;
255 }
256
257 if (index >= static_cast<int>(fSamplingSummaries.size()))
258 throw std::runtime_error("ConfidenceBelt::GetAcceptanceRegion: Sampling summaries are not filled yet. Switch on NeymanConstruction::CreateConfBelt() or FeldmanCousins::CreateConfBelt().");
259
260 return &(fSamplingSummaries[index].GetAcceptanceRegion());
261 }
262 else {
263 std::cout << "dataset is not initialized properly" << std::endl;
264 }
265
266 return 0;
267
268}
269
270////////////////////////////////////////////////////////////////////////////////
271/// returns list of parameters
272
274{
275 return new RooArgSet(*(fParameterPoints->get()));
276}
277
278////////////////////////////////////////////////////////////////////////////////
279
281{
282 if (parameterPoint.getSize() != fParameterPoints->get()->getSize() ) {
283 std::cout << "size is wrong, parameters don't match" << std::endl;
284 return false;
285 }
286 if ( ! parameterPoint.equals( *(fParameterPoints->get() ) ) ) {
287 std::cout << "size is ok, but parameters don't match" << std::endl;
288 return false;
289 }
290 return true;
291}
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
char name[80]
Definition: TGX11.cxx:110
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
Int_t getSize() const
Return the number of elements in the collection.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:61
virtual const RooArgSet * get() const
Definition: RooAbsData.h:105
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
Int_t getIndex(const RooAbsCollection &coord, bool fast=false) const
Calculate bin number of the given coordinates.
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
std::vector< SamplingSummary > fSamplingSummaries
ConfidenceBelt()
Default constructor.
virtual RooArgSet * GetParameters() const
do we want it to return list of parameters
AcceptanceRegion * GetAcceptanceRegion(RooArgSet &, double cl=-1., double leftside=-1.)
Method to determine if a parameter point is in the interval.
double GetAcceptanceRegionMax(RooArgSet &, double cl=-1., double leftside=-1.)
void AddAcceptanceRegion(RooArgSet &, AcceptanceRegion region, double cl=-1., double leftside=-1.)
add after creating a region
SamplingSummaryLookup fSamplingSummaryLookup
std::vector< double > ConfidenceLevels() const
double GetAcceptanceRegionMin(RooArgSet &, double cl=-1., double leftside=-1.)
~ConfidenceBelt() override
Destructor.
bool CheckParameters(RooArgSet &) const
check if parameters are correct. (dummy implementation to start)
void Add(double cl, double leftside)
Int_t GetLookupIndex(double cl, double leftside)
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Namespace for the RooStats classes.
Definition: Asimov.h:19
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754.
Definition: TMath.h:900
Definition: tree.py:1