Logo ROOT   master
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 
23 ConfidenceBelt is a concrete implementation of the ConfInterval interface.
24 It implements simple general purpose interval of arbitrary dimensions and shape.
25 It does not assume the interval is connected.
26 It uses either a RooDataSet (eg. a list of parameter points in the interval) or
27 a RooDataHist (eg. a Histogram-like object for small regions of the parameter space) to
28 store the interval.
29 
30 */
31 
33 
34 #include "RooDataSet.h"
35 #include "RooDataHist.h"
36 
37 #include "RooStats/RooStatsUtils.h"
38 
40 
41 using namespace RooStats;
42 using 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 
63 ConfidenceBelt::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 
79 ConfidenceBelt::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 
109 vector<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 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
ConfidenceBelt()
Default constructor.
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
AcceptanceRegion * GetAcceptanceRegion(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
Method to determine if a parameter point is in the interval.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:891
void Add(Double_t cl, Double_t leftside)
Int_t getIndex(const RooArgSet &coord, Bool_t fast=kFALSE)
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
virtual ~ConfidenceBelt()
Destructor.
STL namespace.
void AddAcceptanceRegion(RooArgSet &, AcceptanceRegion region, Double_t cl=-1., Double_t leftside=-1.)
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
virtual RooArgSet * GetParameters() const
returns list of parameters
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
SamplingSummaryLookup fSamplingSummaryLookup
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
Int_t getSize() const
Int_t GetLookupIndex(Double_t cl, Double_t leftside)
TObject * Next()
Definition: TCollection.h:249
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
Namespace for the RooStats classes.
Definition: Asimov.h:19
#define ClassImp(name)
Definition: Rtypes.h:361
virtual Int_t numEntries() const
Return the number of bins.
std::vector< Double_t > ConfidenceLevels() 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.
Definition: RooArgSet.cxx:472
Bool_t CheckParameters(RooArgSet &) const
Definition: tree.py:1
Double_t GetAcceptanceRegionMax(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
std::vector< SamplingSummary > fSamplingSummaries
char name[80]
Definition: TGX11.cxx:109
Double_t GetAcceptanceRegionMin(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)