Logo ROOT   master
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 
14 This class allow to compute in the ToyMcStudy framework the ProfileLikelihood
15 upper 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"
28 #include "RooStats/ConfInterval.h"
33 #include "TCanvas.h"
34 #include "RooMinuit.h"
35 #include "RooNLLVar.h"
36 #include "RooCmdArg.h"
37 #include "RooRealVar.h"
38 
39 using namespace std;
40 
42 
43 
44 using namespace RooStats ;
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 
48 UpperLimitMCSModule::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 
65  RooAbsMCStudyModule(other),
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));
191  RooStats::ConfInterval* pllint=plc.GetInterval();
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 
207  _data->add(RooArgSet(*_ul));
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 }
std::string GetName(const std::string &scope_name)
Definition: Cppyy.cxx:150
#define coutE(a)
Definition: RooMsgService.h:33
RooAbsData * genSample()
RooAbsMCStudyModule is a base class for add-on modules to RooMCStudy that can perform additional calc...
virtual LikelihoodInterval * GetInterval() const
Return a likelihood interval.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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) ...
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
Bool_t initializeRun(Int_t)
Initialize module at beginning of RooCMStudy run.
Basic string class.
Definition: TString.h:131
STL namespace.
This class allow to compute in the ToyMcStudy framework the ProfileLikelihood upper limit for each to...
virtual void reset()
Definition: RooAbsData.cxx:313
virtual ~UpperLimitMCSModule()
Destructor.
RooDataSet * finalizeRun()
Return auxiliary dataset with results of delta(-log(L)) calculations of this module so that it is mer...
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
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 &#39;value&#39;.
Definition: RooRealVar.cxx:225
RooAbsArg * first() const
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
char * Form(const char *fmt,...)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
const Bool_t kFALSE
Definition: RtypesCore.h:90
RooStats::ProfileLikelihoodCalculator * _plc
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Namespace for the RooStats classes.
Definition: Asimov.h:19
#define ClassImp(name)
Definition: Rtypes.h:361
RooAbsArg * find(const char *name) const
Find object with given name in list.
Bool_t initializeInstance()
Initialize module after attachment to RooMCStudy object.
UpperLimitMCSModule(const RooArgSet *poi, Double_t CL=0.95)
Definition: first.py:1
const Bool_t kTRUE
Definition: RtypesCore.h:89
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
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 &#39;data&#39; argset, to the data set...
const char * Data() const
Definition: TString.h:364