Logo ROOT   6.10/09
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  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 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TIterator * createIterator(Bool_t dir=kIterForward) const
ConfidenceBelt()
Default constructor.
virtual const RooArgSet * get() const
Definition: RooAbsData.h:77
AcceptanceRegion * GetAcceptanceRegion(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
Method to determine if a parameter point is in the interval.
Double_t QuietNaN()
Definition: TMath.h:784
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
void Add(Double_t cl, Double_t leftside)
Int_t getIndex(const RooArgSet &coord, Bool_t fast=kFALSE)
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
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.)
RooDataSet 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 fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
Int_t getSize() const
Int_t GetLookupIndex(Double_t cl, Double_t leftside)
TObject * Next()
Definition: TCollection.h:153
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0)
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:336
virtual Int_t numEntries() const
Return the number of bins.
double Double_t
Definition: RtypesCore.h:55
std::vector< Double_t > ConfidenceLevels() const
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TNamed.cxx:65
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:527
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
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:269
Double_t GetAcceptanceRegionMin(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)