Logo ROOT  
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/** \class RooStats::UpperLimitMCSModule
12 \ingroup Roostats
13
14This class allow to compute in the ToyMcStudy framework the ProfileLikelihood
15upper limit for each toy-MC sample generated
16
17*/
18
19#include "Riostream.h"
20
21#include "RooDataSet.h"
22//#include "RooRealVar.h"
23#include "TString.h"
24//#include "RooFit.h"
25#include "RooFitResult.h"
27#include "RooMsgService.h"
33#include "TCanvas.h"
34#include "RooMinuit.h"
35#include "RooNLLVar.h"
36#include "RooCmdArg.h"
37#include "RooRealVar.h"
38
39using namespace std;
40
42
43
44using namespace RooStats ;
45
46////////////////////////////////////////////////////////////////////////////////
47
48UpperLimitMCSModule::UpperLimitMCSModule(const RooArgSet* poi, Double_t CL) :
49 RooAbsMCStudyModule(Form("UpperLimitMCSModule_%s",poi->first()->GetName()),Form("UpperLimitMCSModule_%s",poi->first()->GetName())),
50 _parName(poi->first()->GetName()),
51 _plc(0),_ul(0),_poi(0), _data(0),_cl(CL), _model(0)
52{
53 std::cout<<"RooUpperLimitConstructor ParName:"<<_parName<<std::endl;
54 std::cout<<"RooUpperLimitConstructor CL:"<<_cl<<std::endl;
55 // Constructor of module with parameter to be interpreted as nSignal and the value of the
56 // null hypothesis for nSignal (usually zero)
57}
58
59
60
61////////////////////////////////////////////////////////////////////////////////
62/// Copy constructor
63
66 _parName(other._poi->first()->GetName()),
67 _plc(0),_ul(0),_poi(other._poi), _data(0), _cl(other._cl), _model(other._model)
68{
69}
70
71////////////////////////////////////////////////////////////////////////////////
72/// Destructor
73
75{
76
77 if (_plc) {
78 delete _plc ;
79 }
80 if (_data) {
81 delete _data ;
82 }
83 if(_ul){
84 delete _ul;
85 }
86 if(_poi){
87 delete _poi;
88 }
89 if (_model){
90 delete _model;
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/// Initialize module at beginning of RooCMStudy run
126
128{
129 _data->reset() ;
130 return kTRUE ;
131}
132
133////////////////////////////////////////////////////////////////////////////////
134/// Return auxiliary dataset with results of delta(-log(L))
135/// calculations of this module so that it is merged with
136/// RooMCStudy::fitParDataSet() by RooMCStudy
137
139{
140 return _data ;
141}
142
143////////////////////////////////////////////////////////////////////////////////
144
145// Bool_t UpperLimitMCSModule::processAfterFit(Int_t /*sampleNum*/)
146// {
147// // Save likelihood from nominal fit, fix chosen parameter to its
148// // null hypothesis value and rerun fit Save difference in likelihood
149// // and associated Gaussian significance in auxiliary dataset
150
151// RooRealVar* par = static_cast<RooRealVar*>(fitParams()->find(_parName.c_str())) ;
152// par->setVal(_nullValue) ;
153// par->setConstant(kTRUE) ;
154// RooFitResult* frnull = refit() ;
155// par->setConstant(kFALSE) ;
156
157// _nll0h->setVal(frnull->minNll()) ;
158
159// Double_t deltaLL = (frnull->minNll() - nllVar()->getVal()) ;
160// Double_t signif = deltaLL>0 ? sqrt(2*deltaLL) : -sqrt(-2*deltaLL) ;
161// _sig0h->setVal(signif) ;
162// _dll0h->setVal(deltaLL) ;
163
164
165// _data->add(RooArgSet(*_nll0h,*_dll0h,*_sig0h)) ;
166
167// delete frnull ;
168// return kTRUE ;
169
170// }
171
172////////////////////////////////////////////////////////////////////////////////
173
175 std::cout<<"after generation Test"<<std::endl;
176
177 if (!fitInitParams() || !genSample() || !fitParams() || !fitModel() ) return kFALSE;
178
179 static_cast<RooRealVar*>(_poi->first())->setVal(static_cast<RooRealVar*>(fitInitParams()->find(_parName.c_str()))->getVal());
180
181 //_poi->first()->Print();
182 static_cast<RooRealVar*>(_poi->first())->setBins(1000);
183 //fitModel()->Print("v");
184
185 std::cout<<"generated Entries:"<<genSample()->numEntries()<<std::endl;
186
188
189 //PLC calculates intervals. for one sided ul multiply testsize by two
190 plc.SetTestSize(2*(1-_cl));
192
193 if (!pllint) return kFALSE;
194
195 std::cout<<"poi value: "<<((RooRealVar*)( _poi->first()))->getVal()<<std::endl;
196 std::cout<<(static_cast<RooRealVar*>((fitParams()->find(_parName.c_str()))))->getVal()<<std::endl;
197 std::cout<<((RooStats::LikelihoodInterval*)pllint)->UpperLimit((RooRealVar&)*(_poi->first()))<<std::endl;
198
199
200 //Go to the fit Value for zour POI to make sure upper limit works correct.
201 //fitModel()->fitTo(*genSample());
202
203
204
205 _ul->setVal(((RooStats::LikelihoodInterval*)pllint)->UpperLimit(static_cast<RooRealVar&>(*(fitParams()->find(_parName.c_str())))));
206
208 std::cout<<"UL:"<<_ul->getVal()<<std::endl;
209// if (_ul->getVal()<1){
210
211// RooStats::LikelihoodIntervalPlot plotpll((RooStats::LikelihoodInterval*) pllint);
212// TCanvas c1;
213// plotpll.Draw();
214// c1.Print("test.ps");
215// std::cout<<" UL<1 whats going on here?"<<std::endl;
216// abort();
217// }
218
219 delete pllint;
220
221
222 return kTRUE;
223}
#define coutE(a)
Definition: RooMsgService.h:34
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
char * Form(const char *fmt,...)
RooAbsArg * first() const
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual void reset()
Definition: RooAbsData.cxx:314
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:307
RooAbsMCStudyModule is a base class for add-on modules to RooMCStudy that can perform additional calc...
RooAbsData * genSample()
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:278
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)
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
virtual LikelihoodInterval * GetInterval() const
Return a likelihood interval.
This class allow to compute in the ToyMcStudy framework the ProfileLikelihood upper limit for each to...
Bool_t initializeRun(Int_t)
Initialize module at beginning of RooCMStudy run.
Bool_t initializeInstance()
Initialize module after attachment to RooMCStudy object.
UpperLimitMCSModule(const RooArgSet *poi, Double_t CL=0.95)
virtual ~UpperLimitMCSModule()
Destructor.
RooStats::ProfileLikelihoodCalculator * _plc
RooDataSet * finalizeRun()
Return auxiliary dataset with results of delta(-log(L)) calculations of this module so that it is mer...
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:150
@ InputArguments
Definition: RooGlobalFunc.h:68
Namespace for the RooStats classes.
Definition: Asimov.h:20
Definition: first.py:1