Logo ROOT   6.16/01
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
40
41using namespace RooStats;
42using namespace std;
43
44////////////////////////////////////////////////////////////////////////////////
45/// Default constructor
46
47ConfidenceBelt::ConfidenceBelt() :
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 return &(fSamplingSummaries.at(index).GetAcceptanceRegion());
233 }
234 else if( tree ){
235 // need a way to get index for given point
236 // RooStats::SetParameters(&parameterPoint, tree->get()); // set tree's parameters to desired values
237 Int_t index = 0; //need something like tree->calcTreeIndex();
238 const RooArgSet* thisPoint = 0;
239 for(index=0; index<tree->numEntries(); ++index){
240 thisPoint = tree->get(index);
241 bool samePoint = true;
242 TIter it = parameterPoint.createIterator();
243 RooRealVar *myarg;
244 while ( samePoint && (myarg = (RooRealVar *)it.Next())) {
245 if(myarg->getVal() != thisPoint->getRealValue(myarg->GetName()))
246 samePoint = false;
247 }
248 if(samePoint)
249 break;
250 }
251
252 return &(fSamplingSummaries.at(index).GetAcceptanceRegion());
253 }
254 else {
255 std::cout << "dataset is not initialized properly" << std::endl;
256 }
257
258 return 0;
259
260}
261
262////////////////////////////////////////////////////////////////////////////////
263/// returns list of parameters
264
266{
267 return new RooArgSet(*(fParameterPoints->get()));
268}
269
270////////////////////////////////////////////////////////////////////////////////
271
273{
274 if (parameterPoint.getSize() != fParameterPoints->get()->getSize() ) {
275 std::cout << "size is wrong, parameters don't match" << std::endl;
276 return false;
277 }
278 if ( ! parameterPoint.equals( *(fParameterPoints->get() ) ) ) {
279 std::cout << "size is ok, but parameters don't match" << std::endl;
280 return false;
281 }
282 return true;
283}
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:363
Int_t getSize() const
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically named contents.
TIterator * createIterator(Bool_t dir=kIterForward) const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
virtual const RooArgSet * get() const
Definition: RooAbsData.h:79
Double_t getVal(const RooArgSet *set=0) const
Evaluate object. Returns either cached value or triggers a recalculation.
Definition: RooAbsReal.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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.
Definition: RooArgSet.cxx:472
RooDataSet is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
virtual Int_t numEntries() const
Return the number of bins.
Int_t getIndex(const RooArgSet &coord, Bool_t fast=kFALSE)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
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()
Definition: TCollection.h:249
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
@(#)root/roostats:$Id$ Author: George Lewis, Kyle Cranmer
Definition: Asimov.h:20
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:889
STL namespace.
Definition: tree.py:1