ROOT  6.06/09
Reference Guide
UpperLimitMCSModule.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Nils Ruthmann
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 
13 #include "Riostream.h"
14 
15 #include "RooDataSet.h"
16 //#include "RooRealVar.h"
17 #include "TString.h"
18 //#include "RooFit.h"
19 #include "RooFitResult.h"
21 #include "RooMsgService.h"
22 #include "RooStats/ConfInterval.h"
27 #include "TCanvas.h"
28 #include "RooMinuit.h"
29 #include "RooNLLVar.h"
30 #include "RooCmdArg.h"
31 #include "RooRealVar.h"
32 
33 
34 
35 using namespace std;
36 
38 
39 
40 using namespace RooStats ;
41 
42 ////////////////////////////////////////////////////////////////////////////////
43 
44 UpperLimitMCSModule::UpperLimitMCSModule(const RooArgSet* poi, Double_t CL) :
45  RooAbsMCStudyModule(Form("UpperLimitMCSModule_%s",poi->first()->GetName()),Form("UpperLimitMCSModule_%s",poi->first()->GetName())),
46  _parName(poi->first()->GetName()),
47  _plc(0),_ul(0),_poi(0), _data(0),_cl(CL), _model(0)
48 {
49  std::cout<<"RooUpperLimitConstructor ParName:"<<_parName<<std::endl;
50  std::cout<<"RooUpperLimitConstructor CL:"<<_cl<<std::endl;
51  // Constructor of module with parameter to be interpreted as nSignal and the value of the
52  // null hypothesis for nSignal (usually zero)
53 }
54 
55 
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// Copy constructor
59 
61  RooAbsMCStudyModule(other),
62  _parName(other._poi->first()->GetName()),
63  _plc(0),_ul(0),_poi(other._poi), _data(0), _cl(other._cl), _model(other._model)
64 {
65 }
66 
67 
68 
69 ////////////////////////////////////////////////////////////////////////////////
70 /// Destructor
71 
73 {
74 
75  if (_plc) {
76  delete _plc ;
77  }
78  if (_data) {
79  delete _data ;
80  }
81  if(_ul){
82  delete _ul;
83  }
84  if(_poi){
85  delete _poi;
86  }
87  if (_model){
88  delete _model;
89  }
90 }
91 
92 
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// Initialize module after attachment to RooMCStudy object
96 
98 {
99  // Check that parameter is also present in fit parameter list of RooMCStudy object
100  if (!fitParams()->find(_parName.c_str())) {
101  coutE(InputArguments) << "UpperLimitMCSModule::initializeInstance:: ERROR: No parameter named " << _parName << " in RooMCStudy!" << endl ;
102  return kFALSE ;
103  }
104 
105  //Construct the ProfileLikelihoodCalculator
106  _poi=new RooArgSet(*(fitParams()->find(_parName.c_str())));
107  std::cout<<"RooUpperLimit Initialize Instance: POI Set:"<<std::endl;
108  _poi->Print("v");
109  std::cout<<"RooUpperLimit Initialize Instance: End:"<<std::endl;
110 
111 
112 
113  TString ulName = Form("ul_%s",_parName.c_str()) ;
114  TString ulTitle = Form("UL for parameter %s",_parName.c_str()) ;
115  _ul = new RooRealVar(ulName.Data(),ulTitle.Data(),0) ;
116 
117 
118  // Create new dataset to be merged with RooMCStudy::fitParDataSet
119  _data = new RooDataSet("ULSigData","Additional data for UL study",RooArgSet(*_ul)) ;
120 
121  return kTRUE ;
122 }
123 
124 
125 
126 ////////////////////////////////////////////////////////////////////////////////
127 /// Initialize module at beginning of RooCMStudy run
128 
130 {
131  _data->reset() ;
132  return kTRUE ;
133 }
134 
135 
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 /// Return auxiliary dataset with results of delta(-log(L))
139 /// calculations of this module so that it is merged with
140 /// RooMCStudy::fitParDataSet() by RooMCStudy
141 
143 {
144  return _data ;
145 }
146 
147 
148 
149 //_____________________________________________________________________________
150 
151 // Bool_t UpperLimitMCSModule::processAfterFit(Int_t /*sampleNum*/)
152 // {
153 // // Save likelihood from nominal fit, fix chosen parameter to its
154 // // null hypothesis value and rerun fit Save difference in likelihood
155 // // and associated Gaussian significance in auxilary dataset
156 
157 // RooRealVar* par = static_cast<RooRealVar*>(fitParams()->find(_parName.c_str())) ;
158 // par->setVal(_nullValue) ;
159 // par->setConstant(kTRUE) ;
160 // RooFitResult* frnull = refit() ;
161 // par->setConstant(kFALSE) ;
162 
163 // _nll0h->setVal(frnull->minNll()) ;
164 
165 // Double_t deltaLL = (frnull->minNll() - nllVar()->getVal()) ;
166 // Double_t signif = deltaLL>0 ? sqrt(2*deltaLL) : -sqrt(-2*deltaLL) ;
167 // _sig0h->setVal(signif) ;
168 // _dll0h->setVal(deltaLL) ;
169 
170 
171 // _data->add(RooArgSet(*_nll0h,*_dll0h,*_sig0h)) ;
172 
173 // delete frnull ;
174 // return kTRUE ;
175 
176 // }
177 
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 
182  std::cout<<"after generation Test"<<std::endl;
183 
184  if (!fitInitParams() || !genSample() || !fitParams() || !fitModel() ) return kFALSE;
185 
186  static_cast<RooRealVar*>(_poi->first())->setVal(static_cast<RooRealVar*>(fitInitParams()->find(_parName.c_str()))->getVal());
187 
188  //_poi->first()->Print();
189  static_cast<RooRealVar*>(_poi->first())->setBins(1000);
190  //fitModel()->Print("v");
191 
192  std::cout<<"generated Entries:"<<genSample()->numEntries()<<std::endl;
193 
195 
196  //PLC calculates intervals. for one sided ul multiply testsize by two
197  plc.SetTestSize(2*(1-_cl));
198  RooStats::ConfInterval* pllint=plc.GetInterval();
199 
200  if (!pllint) return kFALSE;
201 
202  std::cout<<"poi value: "<<((RooRealVar*)( _poi->first()))->getVal()<<std::endl;
203  std::cout<<(static_cast<RooRealVar*>((fitParams()->find(_parName.c_str()))))->getVal()<<std::endl;
204  std::cout<<((RooStats::LikelihoodInterval*)pllint)->UpperLimit((RooRealVar&)*(_poi->first()))<<std::endl;
205 
206 
207  //Go to the fit Value for zour POI to make sure upperlimit works correct.
208  //fitModel()->fitTo(*genSample());
209 
210 
211 
212  _ul->setVal(((RooStats::LikelihoodInterval*)pllint)->UpperLimit(static_cast<RooRealVar&>(*(fitParams()->find(_parName.c_str())))));
213 
214  _data->add(RooArgSet(*_ul));
215  std::cout<<"UL:"<<_ul->getVal()<<std::endl;
216 // if (_ul->getVal()<1){
217 
218 // RooStats::LikelihoodIntervalPlot plotpll((RooStats::LikelihoodInterval*) pllint);
219 // TCanvas c1;
220 // plotpll.Draw();
221 // c1.Print("test.ps");
222 // std::cout<<" UL<1 whats going on here?"<<std::endl;
223 // abort();
224 // }
225 
226  delete pllint;
227 
228 
229  return kTRUE;
230 }
#define coutE(a)
Definition: RooMsgService.h:35
RooAbsData * genSample()
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual void SetTestSize(Double_t size)
set the size of the test (rate of Type I error) ( Eg. 0.05 for a 95% Confidence Interval) ...
ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface class f...
Bool_t initializeRun(Int_t)
Initialize module at beginning of RooCMStudy run.
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
STL namespace.
This class allow to compute in the ToyMcStudy framework the ProfileLikelihood upper limit for each to...
virtual LikelihoodInterval * GetInterval() const
Return a likelihood interval.
const char * Data() const
Definition: TString.h:349
RooAbsArg * first() const
virtual void reset()
Definition: RooAbsData.cxx:302
virtual ~UpperLimitMCSModule()
Destructor.
RooDataSet * finalizeRun()
Return auxiliary dataset with results of delta(-log(L)) calculations of this module so that it is mer...
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:202
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:291
char * Form(const char *fmt,...)
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 'data' argset, to the data set...
RooStats::ProfileLikelihoodCalculator * _plc
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:44
Namespace for the RooStats classes.
Definition: Asimov.h:20
ClassImp(RooStats::UpperLimitMCSModule)
double Double_t
Definition: RtypesCore.h:55
Bool_t initializeInstance()
Initialize module after attachment to RooMCStudy object.
UpperLimitMCSModule(const RooArgSet *poi, Double_t CL=0.95)
const Bool_t kTRUE
Definition: Rtypes.h:91