Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
40
41using namespace RooStats;
42using namespace std;
43
44////////////////////////////////////////////////////////////////////////////////
45/// Default constructor
46
48 TNamed(), fParameterPoints(0)
49{
50}
51
52////////////////////////////////////////////////////////////////////////////////
53/// Alternate constructor
54
56 TNamed(name,name), fParameterPoints(0)
57{
58}
59
60////////////////////////////////////////////////////////////////////////////////
61/// Alternate constructor
62
63ConfidenceBelt::ConfidenceBelt(const char* name, const char* title) :
64 TNamed(name,title), fParameterPoints(0)
65{
66}
67
68////////////////////////////////////////////////////////////////////////////////
69/// Alternate constructor
70
72 TNamed(name,name), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
73{
74}
75
76////////////////////////////////////////////////////////////////////////////////
77/// Alternate constructor
78
79ConfidenceBelt::ConfidenceBelt(const char* name, const char* title, RooAbsData& data) :
80 TNamed(name,title), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
81{
82}
83
84////////////////////////////////////////////////////////////////////////////////
85/// Destructor
86
88{
89}
90
91////////////////////////////////////////////////////////////////////////////////
92
94 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
95 AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
96 return (region) ? region->GetLowerLimit() : TMath::QuietNaN();
97}
98
99////////////////////////////////////////////////////////////////////////////////
100
102 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
103 AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
104 return (region) ? region->GetUpperLimit() : TMath::QuietNaN();
105}
106
107////////////////////////////////////////////////////////////////////////////////
108
109vector<Double_t> ConfidenceBelt::ConfidenceLevels() const {
110 vector<Double_t> levels;
111 return levels;
112}
113
114////////////////////////////////////////////////////////////////////////////////
115
117 Double_t lower, Double_t upper,
118 Double_t cl, Double_t leftside){
119 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
120
121 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
122 RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
123
124 // cout << "add: " << tree << " " << hist << endl;
125
126 if( !this->CheckParameters(parameterPoint) )
127 std::cout << "problem with parameters" << std::endl;
128
129 Int_t luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
130 // cout << "lookup index = " << luIndex << endl;
131 if(luIndex <0 ) {
132 fSamplingSummaryLookup.Add(cl,leftside);
133 luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
134 cout << "lookup index = " << luIndex << endl;
135 }
136 AcceptanceRegion* thisRegion = new AcceptanceRegion(luIndex, lower, upper);
137
138 if( hist ) {
139 // need a way to get index for given point
140 // Can do this by setting hist's internal parameters to desired values
141 // need a better way
142 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
143 // int index = hist->calcTreeIndex(); // get index
144 int index = hist->getIndex(parameterPoint); // get index
145 // cout << "hist index = " << index << endl;
146
147 // allocate memory if necessary. numEntries is overkill?
148 if((Int_t)fSamplingSummaries.size() <= index) {
149 fSamplingSummaries.reserve( hist->numEntries() );
150 fSamplingSummaries.resize( hist->numEntries() );
151 }
152
153 // set the region for this point (check for duplicate?)
154 fSamplingSummaries.at(index) = *thisRegion;
155 }
156 else if( tree ){
157 // int index = tree->getIndex(parameterPoint);
158 int index = dsIndex;
159 // cout << "tree index = " << index << endl;
160
161 // allocate memory if necessary. numEntries is overkill?
162 if((Int_t)fSamplingSummaries.size() <= index){
163 fSamplingSummaries.reserve( tree->numEntries() );
164 fSamplingSummaries.resize( tree->numEntries() );
165 }
166
167 // set the region for this point (check for duplicate?)
168 fSamplingSummaries.at( index ) = *thisRegion;
169 }
170}
171
172////////////////////////////////////////////////////////////////////////////////
173
175 Double_t cl, Double_t leftside){
176 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
177
178 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
179 RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
180
181 if( !this->CheckParameters(parameterPoint) )
182 std::cout << "problem with parameters" << std::endl;
183
184
185 if( hist ) {
186 // need a way to get index for given point
187 // Can do this by setting hist's internal parameters to desired values
188 // need a better way
189 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
190 // int index = hist->calcTreeIndex(); // get index
191 int index = hist->getIndex(parameterPoint); // get index
192
193 // allocate memory if necessary. numEntries is overkill?
194 if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( hist->numEntries() );
195
196 // set the region for this point (check for duplicate?)
197 fSamplingSummaries.at(index) = region;
198 }
199 else if( tree ){
200 tree->add( parameterPoint ); // assume it's unique for now
201 int index = tree->numEntries() - 1; //check that last point added has index nEntries -1
202 // allocate memory if necessary. numEntries is overkill?
203 if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( tree->numEntries() );
204
205 // set the region for this point (check for duplicate?)
206 fSamplingSummaries.at( index ) = region;
207 }
208}
209
210////////////////////////////////////////////////////////////////////////////////
211/// Method to determine if a parameter point is in the interval
212
214{
215 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
216
217 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
218 RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
219
220 if( !this->CheckParameters(parameterPoint) ){
221 std::cout << "problem with parameters" << std::endl;
222 return 0;
223 }
224
225 if( hist ) {
226 // need a way to get index for given point
227 // Can do this by setting hist's internal parameters to desired values
228 // need a better way
229 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
230 // Int_t index = hist->calcTreeIndex(); // get index
231 int index = hist->getIndex(parameterPoint); // get index
232 if (index >= (int)fSamplingSummaries.size())
233 throw std::runtime_error("ConfidenceBelt::GetAcceptanceRegion: Sampling summaries are not filled yet. Switch on NeymanConstruction::CreateConfBelt() or FeldmanCousins::CreateConfBelt().");
234
235 return &(fSamplingSummaries[index].GetAcceptanceRegion());
236 }
237 else if( tree ){
238 // need a way to get index for given point
239 // RooStats::SetParameters(&parameterPoint, tree->get()); // set tree's parameters to desired values
240 Int_t index = 0; //need something like tree->calcTreeIndex();
241 const RooArgSet* thisPoint = 0;
242 for(index=0; index<tree->numEntries(); ++index){
243 thisPoint = tree->get(index);
244 bool samePoint = true;
245 TIter it = parameterPoint.createIterator();
246 RooRealVar *myarg;
247 while ( samePoint && (myarg = (RooRealVar *)it.Next())) {
248 if(myarg->getVal() != thisPoint->getRealValue(myarg->GetName()))
249 samePoint = false;
250 }
251 if(samePoint)
252 break;
253 }
254
255 if (index >= (int)fSamplingSummaries.size())
256 throw std::runtime_error("ConfidenceBelt::GetAcceptanceRegion: Sampling summaries are not filled yet. Switch on NeymanConstruction::CreateConfBelt() or FeldmanCousins::CreateConfBelt().");
257
258 return &(fSamplingSummaries[index].GetAcceptanceRegion());
259 }
260 else {
261 std::cout << "dataset is not initialized properly" << std::endl;
262 }
263
264 return 0;
265
266}
267
268////////////////////////////////////////////////////////////////////////////////
269/// returns list of parameters
270
272{
273 return new RooArgSet(*(fParameterPoints->get()));
274}
275
276////////////////////////////////////////////////////////////////////////////////
277
279{
280 if (parameterPoint.getSize() != fParameterPoints->get()->getSize() ) {
281 std::cout << "size is wrong, parameters don't match" << std::endl;
282 return false;
283 }
284 if ( ! parameterPoint.equals( *(fParameterPoints->get() ) ) ) {
285 std::cout << "size is ok, but parameters don't match" << std::endl;
286 return false;
287 }
288 return true;
289}
#define ClassImp(name)
Definition Rtypes.h:364
char name[80]
Definition TGX11.cxx:110
Int_t getSize() const
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:49
virtual const RooArgSet * get() const
Definition RooAbsData.h:92
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:91
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
The RooDataHist is a container class to hold N-dimensional binned data.
Definition RooDataHist.h:37
Int_t getIndex(const RooArgSet &coord, Bool_t fast=false) const
Calculate bin number of the given coordinates.
Int_t numEntries() const override
Return the number of bins.
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:33
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
std::vector< SamplingSummary > fSamplingSummaries
std::vector< Double_t > ConfidenceLevels() const
ConfidenceBelt()
Default constructor.
virtual RooArgSet * GetParameters() const
returns list of parameters
virtual ~ConfidenceBelt()
Destructor.
AcceptanceRegion * GetAcceptanceRegion(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
Method to determine if a parameter point is in the interval.
Double_t GetAcceptanceRegionMin(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
SamplingSummaryLookup fSamplingSummaryLookup
Bool_t CheckParameters(RooArgSet &) const
Double_t GetAcceptanceRegionMax(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
void AddAcceptanceRegion(RooArgSet &, AcceptanceRegion region, Double_t cl=-1., Double_t leftside=-1.)
void Add(Double_t cl, Double_t leftside)
Int_t GetLookupIndex(Double_t cl, Double_t leftside)
TObject * Next()
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
Namespace for the RooStats classes.
Definition Asimov.h:19
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition TMath.h:901
Definition tree.py:1