Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
NumberCountingPdfFactory.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer 28/07/2008
3
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12////////////////////////////////////////////////////////////////////////////////
13
14/** \class RooStats::NumberCountingPdfFactory
15 \ingroup Roostats
16
17A factory for building PDFs and data for a number counting combination.
18The factory produces a PDF for N channels with uncorrelated background
19uncertainty. Correlations can be added by extending this PDF with additional terms.
20The factory relates the signal in each channel to a master signal strength times the
21expected signal in each channel. Thus, the final test is performed on the master signal strength.
22This yields a more powerful test than letting signal in each channel be independent.
23
24The problem has been studied in these references:
25
26 - http://arxiv.org/abs/physics/0511028
27 - http://arxiv.org/abs/physics/0702156
28 - http://cdsweb.cern.ch/record/1099969?ln=en
29
30One can incorporate uncertainty on the expected signal by adding additional terms.
31For the future, perhaps this factory should be extended to include the efficiency terms automatically.
32
33*/
34
36
38
39#include "RooRealVar.h"
40#include "RooAddition.h"
41#include "RooProduct.h"
42#include "RooDataSet.h"
43#include "RooProdPdf.h"
44#include "RooFitResult.h"
45#include "RooPoisson.h"
46#include "RooGlobalFunc.h"
47#include "RooCmdArg.h"
48#include "RooWorkspace.h"
49#include "RooMsgService.h"
50#include "TTree.h"
51#include <sstream>
52
53
54
56
57
58////////////////////////////////////////////////////////////////////////////////
59/// This method produces a PDF for N channels with uncorrelated background
60/// uncertainty. It relates the signal in each channel to a master signal strength times the
61/// expected signal in each channel.
62///
63/// For the future, perhaps this method should be extended to include the efficiency terms automatically.
64
65void RooStats::NumberCountingPdfFactory::AddModel(double *sig, Int_t nbins, RooWorkspace *ws, const char *pdfName,
66 const char *muName)
67{
68 TList likelihoodFactors;
69
70 // double MaxSigma = 8; // Needed to set ranges for variables.
71
72 RooRealVar* masterSignal =
73 new RooRealVar(muName,"masterSignal",1., 0., 3.);
74
75
76 // loop over individual channels
77 for(Int_t i=0; i<nbins; ++i){
78 // need to name the variables dynamically, so please forgive the string manipulation and focus on values & ranges.
79 std::stringstream str;
80 str<<"_"<<i;
81 RooRealVar* expectedSignal =
82 new RooRealVar(("expected_s"+str.str()).c_str(),("expected_s"+str.str()).c_str(),sig[i], 0., 2*sig[i]);
83 expectedSignal->setConstant(true);
84
85 RooProduct* s =
86 new RooProduct(("s"+str.str()).c_str(),("s"+str.str()).c_str(), RooArgSet(*masterSignal, *expectedSignal));
87
88 RooRealVar* b =
89 new RooRealVar(("b"+str.str()).c_str(),("b"+str.str()).c_str(), .5, 0.,1.);
90 RooRealVar* tau =
91 new RooRealVar(("tau"+str.str()).c_str(),("tau"+str.str()).c_str(), .5, 0., 1.);
92 tau->setConstant(true);
93
94 RooAddition* splusb =
95 new RooAddition(("splusb"+str.str()).c_str(),("s"+str.str()+"+"+"b"+str.str()).c_str(),
96 RooArgSet(*s,*b));
97 RooProduct* bTau =
98 new RooProduct(("bTau"+str.str()).c_str(),("b*tau"+str.str()).c_str(), RooArgSet(*b, *tau));
99 RooRealVar* x =
100 new RooRealVar(("x"+str.str()).c_str(),("x"+str.str()).c_str(), 0.5 , 0., 1.);
101 RooRealVar* y =
102 new RooRealVar(("y"+str.str()).c_str(),("y"+str.str()).c_str(), 0.5, 0., 1.);
103
104
105 RooPoisson* sigRegion =
106 new RooPoisson(("sigRegion"+str.str()).c_str(),("sigRegion"+str.str()).c_str(), *x,*splusb);
107
108 //LM: need to set noRounding since y can take non integer values
109 RooPoisson* sideband =
110 new RooPoisson(("sideband"+str.str()).c_str(),("sideband"+str.str()).c_str(), *y,*bTau,true);
111
112 likelihoodFactors.Add(sigRegion);
113 likelihoodFactors.Add(sideband);
114
115 }
116
117 RooArgSet likelihoodFactorSet(likelihoodFactors);
118 RooProdPdf joint(pdfName,"joint", likelihoodFactorSet );
119 // joint.printCompactTree();
120
121 // add this PDF to workspace.
122 // Need to do import into workspace now to get all the structure imported as well.
123 // Just returning the WS will loose the rest of the structure b/c it will go out of scope
125 ws->import(joint);
127}
128
129////////////////////////////////////////////////////////////////////////////////
130/// Arguments are an array of expected signal, expected background, and relative
131/// background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
132
133void RooStats::NumberCountingPdfFactory::AddExpData(double *sig, double *back, double *back_syst, Int_t nbins,
134 RooWorkspace *ws, const char *dsName)
135{
136 std::vector<double> mainMeas(nbins);
137
138 // loop over channels
139 for(Int_t i=0; i<nbins; ++i){
140 mainMeas[i] = sig[i] + back[i];
141 }
142 return AddData(&mainMeas[0], back, back_syst, nbins, ws, dsName);
143}
144
145////////////////////////////////////////////////////////////////////////////////
146/// Arguments are an array of expected signal, expected background, and relative
147/// ratio of background expected in the sideband to that expected in signal region,
148/// and the number of channels.
149
150void RooStats::NumberCountingPdfFactory::AddExpDataWithSideband(double *sigExp, double *backExp, double *tau,
151 Int_t nbins, RooWorkspace *ws, const char *dsName)
152{
153 std::vector<double> mainMeas(nbins);
154 std::vector<double> sideband(nbins);
155 for(Int_t i=0; i<nbins; ++i){
156 mainMeas[i] = sigExp[i] + backExp[i];
157 sideband[i] = backExp[i]*tau[i];
158 }
159 return AddDataWithSideband(&mainMeas[0], &sideband[0], tau, nbins, ws, dsName);
160
161}
162
163////////////////////////////////////////////////////////////////////////////////
164/// need to be careful here that the range of observable in the dataset is consistent with the one in the workspace
165/// don't rescale unless necessary. If it is necessary, then rescale by x10 or a defined maximum.
166
169{
170 return SafeObservableCreation(ws, varName, value, 10. * value);
171}
172
173////////////////////////////////////////////////////////////////////////////////
174/// need to be careful here that the range of observable in the dataset is consistent with the one in the workspace
175/// don't rescale unless necessary. If it is necessary, then rescale by x10 or a defined maximum.
176
178 double value, double maximum)
179{
180 RooRealVar* x = ws->var( varName );
181 if( !x )
182 x = new RooRealVar(varName, varName, value, 0, maximum );
183 if( x->getMax() < value )
184 x->setMax( std::max(x->getMax(), 10*value ) );
185 x->setVal( value );
186
187 return x;
188}
189
190
191////////////////////////////////////////////////////////////////////////////////
192/// Arguments are an array of results from a main measurement, a measured background,
193/// and relative background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
194
195void RooStats::NumberCountingPdfFactory::AddData(double *mainMeas, double *back, double *back_syst, Int_t nbins,
196 RooWorkspace *ws, const char *dsName)
197{
198 double MaxSigma = 8; // Needed to set ranges for variables.
199
200 TList observablesCollection;
201
202 TTree* tree = new TTree();
203 std::vector<double> xForTree(nbins);
204 std::vector<double> yForTree(nbins);
205
206 // loop over channels
207 for(Int_t i=0; i<nbins; ++i){
208 std::stringstream str;
209 str<<"_"<<i;
210
211 //double _tau = 1./back[i]/back_syst[i]/back_syst[i];
212 // LM: compute tau correctly for the Gamma distribution : mode = tau*b and variance is (tau*b+1)
213 double err = back_syst[i];
214 double _tau = (1.0 + sqrt(1 + 4 * err * err))/ (2. * err * err)/ back[i];
215
216 RooRealVar* tau = SafeObservableCreation(ws, ("tau"+str.str()).c_str(), _tau );
217
218 oocoutW(ws,ObjectHandling) << "NumberCountingPdfFactory: changed value of " << tau->GetName() << " to " << tau->getVal() <<
219 " to be consistent with background and its uncertainty. " <<
220 " Also stored these values of tau into workspace with name . " << (std::string{tau->GetName()}+dsName).c_str() <<
221 " if you test with a different dataset, you should adjust tau appropriately.\n"<< std::endl;
223 ws->import(*(static_cast<RooRealVar*>(tau->clone( (std::string{tau->GetName()}+dsName).c_str() )) ) );
225
226 // need to be careful
227 RooRealVar* x = SafeObservableCreation(ws, ("x"+str.str()).c_str(), mainMeas[i]);
228
229 // need to be careful
230 RooRealVar* y = SafeObservableCreation(ws, ("y"+str.str()).c_str(), back[i]*_tau );
231
232 observablesCollection.Add(x);
233 observablesCollection.Add(y);
234
235 xForTree[i] = mainMeas[i];
236 yForTree[i] = back[i]*_tau;
237
238 tree->Branch(("x"+str.str()).c_str(), &xForTree[i] ,("x"+str.str()+"/D").c_str());
239 tree->Branch(("y"+str.str()).c_str(), &yForTree[i] ,("y"+str.str()+"/D").c_str());
240
241 ws->var("b"+str.str())->setMax( 1.2*back[i]+MaxSigma*(sqrt(back[i])+back[i]*back_syst[i]) );
242 ws->var("b"+str.str())->setVal( back[i] );
243
244 }
245 tree->Fill();
246 // tree->Print();
247 // tree->Scan();
248
249 RooArgList* observableList = new RooArgList(observablesCollection);
250
251 // observableSet->Print();
252 // observableList->Print();
253
254 RooDataSet data{dsName,"Number Counting Data", *observableList, RooFit::Import(*tree)}; // one experiment
255
256 // import hypothetical data
258 ws->import(data);
260
261}
262
263////////////////////////////////////////////////////////////////////////////////
264/// Arguments are an array of expected signal, expected background, and relative
265/// background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
266
267void RooStats::NumberCountingPdfFactory::AddDataWithSideband(double *mainMeas, double *sideband, double *tauForTree,
268 Int_t nbins, RooWorkspace *ws, const char *dsName)
269{
270 double MaxSigma = 8; // Needed to set ranges for variables.
271
272 TList observablesCollection;
273
274 TTree* tree = new TTree();
275
276 std::vector<double> xForTree(nbins);
277 std::vector<double> yForTree(nbins);
278
279
280 // loop over channels
281 for(Int_t i=0; i<nbins; ++i){
282 std::stringstream str;
283 str<<"_"<<i;
284
285 double _tau = tauForTree[i];
286 double back_syst = 1./sqrt(sideband[i]);
287 double back = (sideband[i]/_tau);
288
289
290 RooRealVar* tau = SafeObservableCreation(ws, ("tau"+str.str()).c_str(), _tau );
291
292 oocoutW(ws,ObjectHandling) << "NumberCountingPdfFactory: changed value of " << tau->GetName() << " to " << tau->getVal() <<
293 " to be consistent with background and its uncertainty. " <<
294 " Also stored these values of tau into workspace with name . " << (std::string{tau->GetName()}+dsName).c_str() <<
295 " if you test with a different dataset, you should adjust tau appropriately.\n"<< std::endl;
297 ws->import(*(static_cast<RooRealVar*>(tau->clone( (std::string{tau->GetName()}+dsName).c_str() )) ) );
299
300 // need to be careful
301 RooRealVar* x = SafeObservableCreation(ws, ("x"+str.str()).c_str(), mainMeas[i]);
302
303 // need to be careful
304 RooRealVar* y = SafeObservableCreation(ws, ("y"+str.str()).c_str(), sideband[i] );
305
306
307 observablesCollection.Add(x);
308 observablesCollection.Add(y);
309
310 xForTree[i] = mainMeas[i];
311 yForTree[i] = sideband[i];
312
313 tree->Branch(("x"+str.str()).c_str(), &xForTree[i] ,("x"+str.str()+"/D").c_str());
314 tree->Branch(("y"+str.str()).c_str(), &yForTree[i] ,("y"+str.str()+"/D").c_str());
315
316 ws->var("b"+str.str())->setMax( 1.2*back+MaxSigma*(sqrt(back)+back*back_syst) );
317 ws->var("b"+str.str())->setVal( back );
318
319 }
320 tree->Fill();
321
322 RooArgList observableList{observablesCollection};
323
324 RooDataSet data{dsName,"Number Counting Data", observableList, RooFit::Import(*tree)}; // one experiment
325
326 // import hypothetical data
328 ws->import(data);
330}
#define b(i)
Definition RSha256.hxx:100
#define oocoutW(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
void setConstant(bool value=true)
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
Calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
Definition RooAddition.h:27
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
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
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
Poisson pdf.
Definition RooPoisson.h:19
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:33
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
TObject * clone(const char *newname) const override
Definition RooRealVar.h:48
void setMax(const char *name, double value)
Set maximum of name range to given value.
A factory for building PDFs and data for a number counting combination.
void AddDataWithSideband(double *mainMeas, double *sideband, double *tau, Int_t nbins, RooWorkspace *ws, const char *dsName="ExpectedNumberCountingData")
Arguments are an array of expected signal, expected background, and relative background uncertainty (...
void AddModel(double *sigExp, Int_t nchan, RooWorkspace *ws, const char *pdfName="CombinedPdf", const char *masterSignalName="masterSignal")
This method produces a PDF for N channels with uncorrelated background uncertainty.
void AddData(double *mainMeas, double *bkgMeas, double *db, Int_t nbins, RooWorkspace *ws, const char *dsName="NumberCountingData")
Arguments are an array of results from a main measurement, a measured background, and relative backgr...
void AddExpData(double *sigExp, double *bkgExp, double *db, Int_t nbins, RooWorkspace *ws, const char *dsName="ExpectedNumberCountingData")
Arguments are an array of expected signal, expected background, and relative background uncertainty (...
RooRealVar * SafeObservableCreation(RooWorkspace *ws, const char *varName, double value)
need to be careful here that the range of observable in the dataset is consistent with the one in the...
void AddExpDataWithSideband(double *sigExp, double *bkgExp, double *tau, Int_t nbins, RooWorkspace *ws, const char *dsName="NumberCountingData")
Arguments are an array of expected signal, expected background, and relative ratio of background expe...
Persistable container for RooFit projects.
bool import(const RooAbsArg &arg, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={})
Import a RooAbsArg object, e.g.
RooRealVar * var(RooStringView name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
A doubly linked list.
Definition TList.h:38
void Add(TObject *obj) override
Definition TList.h:83
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
A TTree represents a columnar dataset.
Definition TTree.h:79
RooCmdArg Import(const char *state, TH1 &histo)
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17