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 "TString.h"
23#include "RooFitResult.h"
25#include "RooMsgService.h"
31#include "TCanvas.h"
32#include "RooCmdArg.h"
33#include "RooRealVar.h"
34
35using namespace std;
36
38
39
40using namespace RooStats ;
41
42////////////////////////////////////////////////////////////////////////////////
43
44UpperLimitMCSModule::UpperLimitMCSModule(const RooArgSet* poi, double 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
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/// Destructor
69
71{
72
73 if (_plc) {
74 delete _plc ;
75 }
76 if (_data) {
77 delete _data ;
78 }
79 if(_ul){
80 delete _ul;
81 }
82 if(_poi){
83 delete _poi;
84 }
85 if (_model){
86 delete _model;
87 }
88}
89
90////////////////////////////////////////////////////////////////////////////////
91/// Initialize module after attachment to RooMCStudy object
92
94{
95 // Check that parameter is also present in fit parameter list of RooMCStudy object
96 if (!fitParams()->find(_parName.c_str())) {
97 coutE(InputArguments) << "UpperLimitMCSModule::initializeInstance:: ERROR: No parameter named " << _parName << " in RooMCStudy!" << endl ;
98 return false ;
99 }
100
101 //Construct the ProfileLikelihoodCalculator
102 _poi=new RooArgSet(*(fitParams()->find(_parName.c_str())));
103 std::cout<<"RooUpperLimit Initialize Instance: POI Set:"<<std::endl;
104 _poi->Print("v");
105 std::cout<<"RooUpperLimit Initialize Instance: End:"<<std::endl;
106
107
108
109 TString ulName = Form("ul_%s",_parName.c_str()) ;
110 TString ulTitle = Form("UL for parameter %s",_parName.c_str()) ;
111 _ul = new RooRealVar(ulName.Data(),ulTitle.Data(),0) ;
112
113
114 // Create new dataset to be merged with RooMCStudy::fitParDataSet
115 _data = new RooDataSet("ULSigData","Additional data for UL study",RooArgSet(*_ul)) ;
116
117 return true ;
118}
119
120////////////////////////////////////////////////////////////////////////////////
121/// Initialize module at beginning of RooCMStudy run
122
124{
125 _data->reset() ;
126 return true ;
127}
128
129////////////////////////////////////////////////////////////////////////////////
130/// Return auxiliary dataset with results of delta(-log(L))
131/// calculations of this module so that it is merged with
132/// RooMCStudy::fitParDataSet() by RooMCStudy
133
135{
136 return _data ;
137}
138
139////////////////////////////////////////////////////////////////////////////////
140
141// bool UpperLimitMCSModule::processAfterFit(Int_t /*sampleNum*/)
142// {
143// // Save likelihood from nominal fit, fix chosen parameter to its
144// // null hypothesis value and rerun fit Save difference in likelihood
145// // and associated Gaussian significance in auxiliary dataset
146
147// RooRealVar* par = static_cast<RooRealVar*>(fitParams()->find(_parName.c_str())) ;
148// par->setVal(_nullValue) ;
149// par->setConstant(true) ;
150// RooFitResult* frnull = refit() ;
151// par->setConstant(false) ;
152
153// _nll0h->setVal(frnull->minNll()) ;
154
155// double deltaLL = (frnull->minNll() - nllVar()->getVal()) ;
156// double signif = deltaLL>0 ? sqrt(2*deltaLL) : -sqrt(-2*deltaLL) ;
157// _sig0h->setVal(signif) ;
158// _dll0h->setVal(deltaLL) ;
159
160
161// _data->add(RooArgSet(*_nll0h,*_dll0h,*_sig0h)) ;
162
163// delete frnull ;
164// return true ;
165
166// }
167
168////////////////////////////////////////////////////////////////////////////////
169
171 std::cout<<"after generation Test"<<std::endl;
172
173 if (!fitInitParams() || !genSample() || !fitParams() || !fitModel() ) return false;
174
175 static_cast<RooRealVar*>(_poi->first())->setVal(static_cast<RooRealVar*>(fitInitParams()->find(_parName.c_str()))->getVal());
176
177 //_poi->first()->Print();
178 static_cast<RooRealVar*>(_poi->first())->setBins(1000);
179 //fitModel()->Print("v");
180
181 std::cout<<"generated Entries:"<<genSample()->numEntries()<<std::endl;
182
184
185 //PLC calculates intervals. for one sided ul multiply testsize by two
186 plc.SetTestSize(2*(1-_cl));
188
189 if (!pllint) return false;
190
191 std::cout<<"poi value: "<<((RooRealVar*)( _poi->first()))->getVal()<<std::endl;
192 std::cout<<(static_cast<RooRealVar*>((fitParams()->find(_parName.c_str()))))->getVal()<<std::endl;
193 std::cout<<((RooStats::LikelihoodInterval*)pllint)->UpperLimit((RooRealVar&)*(_poi->first()))<<std::endl;
194
195
196 //Go to the fit Value for zour POI to make sure upper limit works correct.
197 //fitModel()->fitTo(*genSample());
198
199
200
201 _ul->setVal(((RooStats::LikelihoodInterval*)pllint)->UpperLimit(static_cast<RooRealVar&>(*(fitParams()->find(_parName.c_str())))));
202
204 std::cout<<"UL:"<<_ul->getVal()<<std::endl;
205// if (_ul->getVal()<1){
206
207// RooStats::LikelihoodIntervalPlot plotpll((RooStats::LikelihoodInterval*) pllint);
208// TCanvas c1;
209// plotpll.Draw();
210// c1.Print("test.ps");
211// std::cout<<" UL<1 whats going on here?"<<std::endl;
212// abort();
213// }
214
215 delete pllint;
216
217
218 return true;
219}
#define coutE(a)
Definition: RooMsgService.h:37
#define ClassImp(name)
Definition: Rtypes.h:375
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2456
RooAbsArg * first() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
virtual void reset()
Definition: RooAbsData.cxx:381
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:374
RooAbsMCStudyModule is a base class for add-on modules to RooMCStudy that can perform additional calc...
RooAbsData * genSample()
Return generate sample.
RooAbsPdf * fitModel()
Return fit model.
RooArgSet * fitParams()
Return current value of parameters of fit model.
RooArgSet * fitInitParams()
Return initial value of parameters of fit model.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
void add(const RooArgSet &row, double weight=1.0, double weightError=0.0) override
Add one ore more rows of data.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:254
void SetTestSize(double size) override
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.
double UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
LikelihoodInterval * GetInterval() const override
Return a likelihood interval.
This class allow to compute in the ToyMcStudy framework the ProfileLikelihood upper limit for each to...
const RooArgSet * _poi
parameters of interest
UpperLimitMCSModule(const RooArgSet *poi, double CL=0.95)
RooDataSet * finalizeRun() override
Return auxiliary dataset with results of delta(-log(L)) calculations of this module so that it is mer...
~UpperLimitMCSModule() override
Destructor.
bool processBetweenGenAndFit(Int_t) override
Method called after generation of toy data sample and resetting of fit parameters to initial values a...
RooStats::ProfileLikelihoodCalculator * _plc
bool initializeRun(Int_t) override
Initialize module at beginning of RooCMStudy run.
std::string _parName
Name of Nsignal parameter.
RooDataSet * _data
Summary dataset to store results.
bool initializeInstance() override
Initialize module after attachment to RooMCStudy object.
Basic string class.
Definition: TString.h:136
const char * Data() const
Definition: TString.h:369
@ InputArguments
Definition: RooGlobalFunc.h:62
Namespace for the RooStats classes.
Definition: Asimov.h:19
Definition: first.py:1