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
39#include <TMath.h>
40
42
43using namespace RooStats;
44using std::cout, std::endl, std::vector;
45
46////////////////////////////////////////////////////////////////////////////////
47/// Alternate constructor
48
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55/// Alternate constructor
56
57ConfidenceBelt::ConfidenceBelt(const char* name, const char* title) :
58 TNamed(name,title)
59{
60}
61
62////////////////////////////////////////////////////////////////////////////////
63/// Alternate constructor
64
66 TNamed(name,name), fParameterPoints(static_cast<RooAbsData*>(data.Clone("PointsToTestForBelt")))
67{
68}
69
70////////////////////////////////////////////////////////////////////////////////
71/// Alternate constructor
72
73ConfidenceBelt::ConfidenceBelt(const char* name, const char* title, RooAbsData& data) :
74 TNamed(name,title), fParameterPoints(static_cast<RooAbsData*>(data.Clone("PointsToTestForBelt")))
75{
76}
77
78////////////////////////////////////////////////////////////////////////////////
79
80double ConfidenceBelt::GetAcceptanceRegionMin(RooArgSet& parameterPoint, double cl, double leftside) {
81 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
82 AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
83 return (region) ? region->GetLowerLimit() : TMath::QuietNaN();
84}
85
86////////////////////////////////////////////////////////////////////////////////
87
88double ConfidenceBelt::GetAcceptanceRegionMax(RooArgSet& parameterPoint, double cl, double leftside) {
89 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
90 AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
91 return (region) ? region->GetUpperLimit() : TMath::QuietNaN();
92}
93
94////////////////////////////////////////////////////////////////////////////////
95
96vector<double> ConfidenceBelt::ConfidenceLevels() const {
97 vector<double> levels;
98 return levels;
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
104 double lower, double upper,
105 double cl, double leftside){
106 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
107
108 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
109 RooDataHist* hist = dynamic_cast<RooDataHist*>(fParameterPoints );
110
111 // cout << "add: " << tree << " " << hist << endl;
112
113 if( !this->CheckParameters(parameterPoint) )
114 std::cout << "problem with parameters" << std::endl;
115
116 Int_t luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
117 // cout << "lookup index = " << luIndex << endl;
118 if(luIndex <0 ) {
119 fSamplingSummaryLookup.Add(cl,leftside);
120 luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
121 cout << "lookup index = " << luIndex << endl;
122 }
123 AcceptanceRegion* thisRegion = new AcceptanceRegion(luIndex, lower, upper);
124
125 if( hist ) {
126 // need a way to get index for given point
127 // Can do this by setting hist's internal parameters to desired values
128 // need a better way
129 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
130 // int index = hist->calcTreeIndex(); // get index
131 int index = hist->getIndex(parameterPoint); // get index
132 // cout << "hist index = " << index << endl;
133
134 // allocate memory if necessary. numEntries is overkill?
135 if((Int_t)fSamplingSummaries.size() <= index) {
136 fSamplingSummaries.reserve( hist->numEntries() );
137 fSamplingSummaries.resize( hist->numEntries() );
138 }
139
140 // set the region for this point (check for duplicate?)
141 fSamplingSummaries.at(index) = *thisRegion;
142 }
143 else if( tree ){
144 // int index = tree->getIndex(parameterPoint);
145 int index = dsIndex;
146 // cout << "tree index = " << index << endl;
147
148 // allocate memory if necessary. numEntries is overkill?
149 if((Int_t)fSamplingSummaries.size() <= index){
150 fSamplingSummaries.reserve( tree->numEntries() );
151 fSamplingSummaries.resize( tree->numEntries() );
152 }
153
154 // set the region for this point (check for duplicate?)
155 fSamplingSummaries.at( index ) = *thisRegion;
156 }
157}
158
159////////////////////////////////////////////////////////////////////////////////
160
162 double cl, double leftside){
163 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
164
165 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
166 RooDataHist* hist = dynamic_cast<RooDataHist*>(fParameterPoints );
167
168 if( !this->CheckParameters(parameterPoint) )
169 std::cout << "problem with parameters" << std::endl;
170
171
172 if( hist ) {
173 // need a way to get index for given point
174 // Can do this by setting hist's internal parameters to desired values
175 // need a better way
176 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
177 // int index = hist->calcTreeIndex(); // get index
178 int index = hist->getIndex(parameterPoint); // get index
179
180 // allocate memory if necessary. numEntries is overkill?
181 if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( hist->numEntries() );
182
183 // set the region for this point (check for duplicate?)
184 fSamplingSummaries.at(index) = region;
185 }
186 else if( tree ){
187 tree->add( parameterPoint ); // assume it's unique for now
188 int index = tree->numEntries() - 1; //check that last point added has index nEntries -1
189 // allocate memory if necessary. numEntries is overkill?
190 if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( tree->numEntries() );
191
192 // set the region for this point (check for duplicate?)
193 fSamplingSummaries.at( index ) = region;
194 }
195}
196
197////////////////////////////////////////////////////////////////////////////////
198/// Method to determine if a parameter point is in the interval
199
200AcceptanceRegion* ConfidenceBelt::GetAcceptanceRegion(RooArgSet &parameterPoint, double cl, double leftside)
201{
202 if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
203
204 RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
205 RooDataHist* hist = dynamic_cast<RooDataHist*>(fParameterPoints );
206
207 if( !this->CheckParameters(parameterPoint) ){
208 std::cout << "problem with parameters" << std::endl;
209 return nullptr;
210 }
211
212 if( hist ) {
213 // need a way to get index for given point
214 // Can do this by setting hist's internal parameters to desired values
215 // need a better way
216 // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
217 // Int_t index = hist->calcTreeIndex(); // get index
218 int index = hist->getIndex(parameterPoint); // get index
219 if (index >= (int)fSamplingSummaries.size())
220 throw std::runtime_error("ConfidenceBelt::GetAcceptanceRegion: Sampling summaries are not filled yet. Switch on NeymanConstruction::CreateConfBelt() or FeldmanCousins::CreateConfBelt().");
221
222 return &(fSamplingSummaries[index].GetAcceptanceRegion());
223 }
224 else if( tree ){
225 // need a way to get index for given point
226 // RooStats::SetParameters(&parameterPoint, tree->get()); // set tree's parameters to desired values
227 Int_t index = 0; //need something like tree->calcTreeIndex();
228 const RooArgSet* thisPoint = nullptr;
229 for(index=0; index<tree->numEntries(); ++index){
230 thisPoint = tree->get(index);
231 bool samePoint = true;
232 for (auto const *myarg : static_range_cast<RooRealVar *>(parameterPoint)) {
233 if (samePoint == false)
234 break;
235 if(myarg->getVal() != thisPoint->getRealValue(myarg->GetName()))
236 samePoint = false;
237 }
238 if(samePoint)
239 break;
240 }
241
242 if (index >= static_cast<int>(fSamplingSummaries.size()))
243 throw std::runtime_error("ConfidenceBelt::GetAcceptanceRegion: Sampling summaries are not filled yet. Switch on NeymanConstruction::CreateConfBelt() or FeldmanCousins::CreateConfBelt().");
244
245 return &(fSamplingSummaries[index].GetAcceptanceRegion());
246 }
247 else {
248 std::cout << "dataset is not initialized properly" << std::endl;
249 }
250
251 return nullptr;
252
253}
254
255////////////////////////////////////////////////////////////////////////////////
256/// returns list of parameters
257
259{
260 return new RooArgSet(*(fParameterPoints->get()));
261}
262
263////////////////////////////////////////////////////////////////////////////////
264
266{
267 if (parameterPoint.size() != fParameterPoints->get()->size() ) {
268 std::cout << "size is wrong, parameters don't match" << std::endl;
269 return false;
270 }
271 if ( ! parameterPoint.equals( *(fParameterPoints->get() ) ) ) {
272 std::cout << "size is ok, but parameters don't match" << std::endl;
273 return false;
274 }
275 return true;
276}
#define ClassImp(name)
Definition Rtypes.h:377
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.
Storage_t::size_type size() const
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:40
Int_t getIndex(const RooAbsCollection &coord, bool fast=false) const
Calculate bin number of the given coordinates.
Container class to hold unbinned data.
Definition RooDataSet.h:57
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
std::vector< SamplingSummary > fSamplingSummaries
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.)
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:902