Logo ROOT   6.08/07
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 
21 
22 
23 
24 #ifndef RooStats_ConfidenceBelt
26 #endif
27 
28 #include "RooDataSet.h"
29 #include "RooDataHist.h"
30 
31 #include "RooStats/RooStatsUtils.h"
32 
34 
35 using namespace RooStats;
36 using namespace std;
37 
38 ////////////////////////////////////////////////////////////////////////////////
39 /// Default constructor
40 
42  TNamed(), fParameterPoints(0)
43 {
44 }
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 /// Alternate constructor
48 
50  TNamed(name,name), fParameterPoints(0)
51 {
52 }
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 /// Alternate constructor
56 
57 ConfidenceBelt::ConfidenceBelt(const char* name, const char* title) :
58  TNamed(name,title), fParameterPoints(0)
59 {
60 }
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 /// Alternate constructor
64 
66  TNamed(name,name), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
67 {
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// Alternate constructor
72 
73 ConfidenceBelt::ConfidenceBelt(const char* name, const char* title, RooAbsData& data) :
74  TNamed(name,title), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
75 {
76 }
77 
78 
79 
80 ////////////////////////////////////////////////////////////////////////////////
81 /// Destructor
82 
84 {
85 }
86 
87 
88 ////////////////////////////////////////////////////////////////////////////////
89 
91  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
92  AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
93  return (region) ? region->GetLowerLimit() : TMath::QuietNaN();
94 }
95 
96 ////////////////////////////////////////////////////////////////////////////////
97 
99  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
100  AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside);
101  return (region) ? region->GetUpperLimit() : TMath::QuietNaN();
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 
106 vector<Double_t> ConfidenceBelt::ConfidenceLevels() const {
107  vector<Double_t> levels;
108  return levels;
109 }
110 
111 ////////////////////////////////////////////////////////////////////////////////
112 
114  Double_t lower, Double_t upper,
115  Double_t cl, Double_t leftside){
116  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
117 
118  RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
119  RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
120 
121  // cout << "add: " << tree << " " << hist << endl;
122 
123  if( !this->CheckParameters(parameterPoint) )
124  std::cout << "problem with parameters" << std::endl;
125 
126  Int_t luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
127  // cout << "lookup index = " << luIndex << endl;
128  if(luIndex <0 ) {
129  fSamplingSummaryLookup.Add(cl,leftside);
130  luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
131  cout << "lookup index = " << luIndex << endl;
132  }
133  AcceptanceRegion* thisRegion = new AcceptanceRegion(luIndex, lower, upper);
134 
135  if( hist ) {
136  // need a way to get index for given point
137  // Can do this by setting hist's internal parameters to desired values
138  // need a better way
139  // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
140  // int index = hist->calcTreeIndex(); // get index
141  int index = hist->getIndex(parameterPoint); // get index
142  // cout << "hist index = " << index << endl;
143 
144  // allocate memory if necessary. numEntries is overkill?
145  if((Int_t)fSamplingSummaries.size() <= index) {
146  fSamplingSummaries.reserve( hist->numEntries() );
147  fSamplingSummaries.resize( hist->numEntries() );
148  }
149 
150  // set the region for this point (check for duplicate?)
151  fSamplingSummaries.at(index) = *thisRegion;
152  }
153  else if( tree ){
154  // int index = tree->getIndex(parameterPoint);
155  int index = dsIndex;
156  // cout << "tree index = " << index << endl;
157 
158  // allocate memory if necessary. numEntries is overkill?
159  if((Int_t)fSamplingSummaries.size() <= index){
160  fSamplingSummaries.reserve( tree->numEntries() );
161  fSamplingSummaries.resize( tree->numEntries() );
162  }
163 
164  // set the region for this point (check for duplicate?)
165  fSamplingSummaries.at( index ) = *thisRegion;
166  }
167 }
168 
169 ////////////////////////////////////////////////////////////////////////////////
170 
172  Double_t cl, Double_t leftside){
173  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
174 
175  RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
176  RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
177 
178  if( !this->CheckParameters(parameterPoint) )
179  std::cout << "problem with parameters" << std::endl;
180 
181 
182  if( hist ) {
183  // need a way to get index for given point
184  // Can do this by setting hist's internal parameters to desired values
185  // need a better way
186  // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
187  // int index = hist->calcTreeIndex(); // get index
188  int index = hist->getIndex(parameterPoint); // get index
189 
190  // allocate memory if necessary. numEntries is overkill?
191  if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( hist->numEntries() );
192 
193  // set the region for this point (check for duplicate?)
194  fSamplingSummaries.at(index) = region;
195  }
196  else if( tree ){
197  tree->add( parameterPoint ); // assume it's unique for now
198  int index = tree->numEntries() - 1; //check that last point added has index nEntries -1
199  // allocate memory if necessary. numEntries is overkill?
200  if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( tree->numEntries() );
201 
202  // set the region for this point (check for duplicate?)
203  fSamplingSummaries.at( index ) = region;
204  }
205 }
206 
207 ////////////////////////////////////////////////////////////////////////////////
208 /// Method to determine if a parameter point is in the interval
209 
211 {
212  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
213 
214  RooDataSet* tree = dynamic_cast<RooDataSet*>( fParameterPoints );
215  RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
216 
217  if( !this->CheckParameters(parameterPoint) ){
218  std::cout << "problem with parameters" << std::endl;
219  return 0;
220  }
221 
222  if( hist ) {
223  // need a way to get index for given point
224  // Can do this by setting hist's internal parameters to desired values
225  // need a better way
226  // RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get()));
227  // Int_t index = hist->calcTreeIndex(); // get index
228  int index = hist->getIndex(parameterPoint); // get index
229  return &(fSamplingSummaries.at(index).GetAcceptanceRegion());
230  }
231  else if( tree ){
232  // need a way to get index for given point
233  // RooStats::SetParameters(&parameterPoint, tree->get()); // set tree's parameters to desired values
234  Int_t index = 0; //need something like tree->calcTreeIndex();
235  const RooArgSet* thisPoint = 0;
236  for(index=0; index<tree->numEntries(); ++index){
237  thisPoint = tree->get(index);
238  bool samePoint = true;
239  TIter it = parameterPoint.createIterator();
240  RooRealVar *myarg;
241  while ( samePoint && (myarg = (RooRealVar *)it.Next())) {
242  if(myarg->getVal() != thisPoint->getRealValue(myarg->GetName()))
243  samePoint = false;
244  }
245  if(samePoint)
246  break;
247  }
248 
249  return &(fSamplingSummaries.at(index).GetAcceptanceRegion());
250  }
251  else {
252  std::cout << "dataset is not initialized properly" << std::endl;
253  }
254 
255  return 0;
256 
257 }
258 
259 ////////////////////////////////////////////////////////////////////////////////
260 /// returns list of parameters
261 
263 {
264  return new RooArgSet(*(fParameterPoints->get()));
265 }
266 
267 ////////////////////////////////////////////////////////////////////////////////
268 
270 {
271  if (parameterPoint.getSize() != fParameterPoints->get()->getSize() ) {
272  std::cout << "size is wrong, parameters don't match" << std::endl;
273  return false;
274  }
275  if ( ! parameterPoint.equals( *(fParameterPoints->get() ) ) ) {
276  std::cout << "size is ok, but parameters don't match" << std::endl;
277  return false;
278  }
279  return true;
280 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
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:622
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:33
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
SamplingSummaryLookup fSamplingSummaryLookup
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
Int_t getSize() const
Int_t GetLookupIndex(Double_t cl, Double_t leftside)
TObject * Next()
Definition: TCollection.h:158
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:279
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
char name[80]
Definition: TGX11.cxx:109
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:269
Double_t GetAcceptanceRegionMin(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)