Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 "RooFitResult.h"
24#include "RooMsgService.h"
30#include "RooRealVar.h"
31
33
34using namespace RooStats ;
35
36////////////////////////////////////////////////////////////////////////////////
37
39 RooAbsMCStudyModule(Form("UpperLimitMCSModule_%s",poi->first()->GetName()),Form("UpperLimitMCSModule_%s",poi->first()->GetName())),
40 _parName(poi->first()->GetName()),
41 _plc(nullptr),_ul(nullptr),_poi(nullptr), _data(nullptr),_cl(CL), _model(nullptr)
42{
43 std::cout<<"RooUpperLimitConstructor ParName:"<<_parName<<std::endl;
44 std::cout<<"RooUpperLimitConstructor CL:"<<_cl<<std::endl;
45 // Constructor of module with parameter to be interpreted as nSignal and the value of the
46 // null hypothesis for nSignal (usually zero)
47}
48
49
50
51////////////////////////////////////////////////////////////////////////////////
52/// Copy constructor
53
56 _parName(other._poi->first()->GetName()),
57 _plc(nullptr),_ul(nullptr),_poi(other._poi), _data(nullptr), _cl(other._cl), _model(other._model)
58{
59}
60
61////////////////////////////////////////////////////////////////////////////////
62/// Destructor
63
65{
66
67 if (_plc) {
68 delete _plc ;
69 }
70 if (_data) {
71 delete _data ;
72 }
73 if(_ul){
74 delete _ul;
75 }
76 if(_poi){
77 delete _poi;
78 }
79 if (_model){
80 delete _model;
81 }
82}
83
84////////////////////////////////////////////////////////////////////////////////
85/// Initialize module after attachment to RooMCStudy object
86
88{
89 // Check that parameter is also present in fit parameter list of RooMCStudy object
90 if (!fitParams()->find(_parName.c_str())) {
91 coutE(InputArguments) << "UpperLimitMCSModule::initializeInstance:: ERROR: No parameter named " << _parName << " in RooMCStudy!" << std::endl ;
92 return false ;
93 }
94
95 //Construct the ProfileLikelihoodCalculator
96 _poi=new RooArgSet(*(fitParams()->find(_parName.c_str())));
97 std::cout<<"RooUpperLimit Initialize Instance: POI Set:"<<std::endl;
98 _poi->Print("v");
99 std::cout<<"RooUpperLimit Initialize Instance: End:"<<std::endl;
100
101
102
103 std::string ulName = "ul_" + _parName;
104 std::string ulTitle = "UL for parameter " + _parName;
105 _ul = new RooRealVar(ulName.c_str(),ulTitle.c_str(),0) ;
106
107
108 // Create new dataset to be merged with RooMCStudy::fitParDataSet
109 _data = new RooDataSet("ULSigData","Additional data for UL study",RooArgSet(*_ul)) ;
110
111 return true ;
112}
113
114////////////////////////////////////////////////////////////////////////////////
115/// Initialize module at beginning of RooCMStudy run
116
118{
119 _data->reset() ;
120 return true ;
121}
122
123////////////////////////////////////////////////////////////////////////////////
124/// Return auxiliary dataset with results of delta(-log(L))
125/// calculations of this module so that it is merged with
126/// RooMCStudy::fitParDataSet() by RooMCStudy
127
129{
130 return _data ;
131}
132
133////////////////////////////////////////////////////////////////////////////////
134
135// bool UpperLimitMCSModule::processAfterFit(Int_t /*sampleNum*/)
136// {
137// // Save likelihood from nominal fit, fix chosen parameter to its
138// // null hypothesis value and rerun fit Save difference in likelihood
139// // and associated Gaussian significance in auxiliary dataset
140
141// RooRealVar* par = static_cast<RooRealVar*>(fitParams()->find(_parName.c_str())) ;
142// par->setVal(_nullValue) ;
143// par->setConstant(true) ;
144// RooFitResult* frnull = refit() ;
145// par->setConstant(false) ;
146
147// _nll0h->setVal(frnull->minNll()) ;
148
149// double deltaLL = (frnull->minNll() - nllVar()->getVal()) ;
150// double signif = deltaLL>0 ? sqrt(2*deltaLL) : -sqrt(-2*deltaLL) ;
151// _sig0h->setVal(signif) ;
152// _dll0h->setVal(deltaLL) ;
153
154
155// _data->add(RooArgSet(*_nll0h,*_dll0h,*_sig0h)) ;
156
157// delete frnull ;
158// return true ;
159
160// }
161
162////////////////////////////////////////////////////////////////////////////////
163
165 std::cout<<"after generation Test"<<std::endl;
166
167 if (!fitInitParams() || !genSample() || !fitParams() || !fitModel() ) return false;
168
169 static_cast<RooRealVar*>(_poi->first())->setVal(static_cast<RooRealVar*>(fitInitParams()->find(_parName.c_str()))->getVal());
170
171 //_poi->first()->Print();
172 static_cast<RooRealVar*>(_poi->first())->setBins(1000);
173 //fitModel()->Print("v");
174
175 std::cout<<"generated Entries:"<<genSample()->numEntries()<<std::endl;
176
178
179 //PLC calculates intervals. for one sided ul multiply testsize by two
180 plc.SetTestSize(2*(1-_cl));
182
183 if (!pllint) return false;
184
185 std::cout<<"poi value: "<<(static_cast<RooRealVar*>(_poi->first()))->getVal()<<std::endl;
186 std::cout<<(static_cast<RooRealVar*>((fitParams()->find(_parName.c_str()))))->getVal()<<std::endl;
187 std::cout<<(static_cast<RooStats::LikelihoodInterval*>(pllint))->UpperLimit(static_cast<RooRealVar&>(*(_poi->first())))<<std::endl;
188
189
190 //Go to the fit Value for zour POI to make sure upper limit works correct.
191 //fitModel()->fitTo(*genSample());
192
193
194
195 _ul->setVal((static_cast<RooStats::LikelihoodInterval*>(pllint))->UpperLimit(static_cast<RooRealVar&>(*(fitParams()->find(_parName.c_str())))));
196
198 std::cout<<"UL:"<<_ul->getVal()<<std::endl;
199// if (_ul->getVal()<1){
200
201// RooStats::LikelihoodIntervalPlot plotpll((RooStats::LikelihoodInterval*) pllint);
202// TCanvas c1;
203// plotpll.Draw();
204// c1.Print("test.ps");
205// std::cout<<" UL<1 whats going on here?"<<std::endl;
206// abort();
207// }
208
209 delete pllint;
210
211
212 return true;
213}
RooAbsData * _data
Pointer to original input dataset.
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:377
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
RooAbsArg * first() const
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
virtual void reset()
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Base class for add-on modules to RooMCStudy that can perform additional calculations on each generate...
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:103
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Container class to hold unbinned data.
Definition RooDataSet.h:57
void add(const RooArgSet &row, double weight, double weightError)
Add one ore more rows of data.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
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.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
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.
Namespace for the RooStats classes.
Definition Asimov.h:19