Logo ROOT   6.08/07
Reference Guide
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 /*
15 BEGIN_HTML
16 <p>
17 A factory for building PDFs and data for a number counting combination.
18 The factory produces a PDF for N channels with uncorrelated background
19 uncertainty. Correlations can be added by extending this PDF with additional terms.
20 The factory relates the signal in each channel to a master signal strength times the
21 expected signal in each channel. Thus, the final test is performed on the master signal strength.
22 This yields a more powerful test than letting signal in each channel be independent.
23 </p>
24 <p>
25 The problem has been studied in these references:
26 <ul>
27  <li> http://arxiv.org/abs/physics/0511028</li>
28  <li> http://arxiv.org/abs/physics/0702156</li>
29  <li> http://cdsweb.cern.ch/record/1099969?ln=en</li>
30 </ul>
31 </p>
32 
33 <p>
34 One can incorporate uncertainty on the expected signal by adding additional terms.
35 For the future, perhaps this factory should be extended to include the efficiency terms automatically.
36 </p>
37 END_HTML
38 */
39 
40 #ifndef RooStats_NumberCountingPdfFactory
42 #endif
43 
44 #ifndef RooStats_RooStatsUtils
45 #include "RooStats/RooStatsUtils.h"
46 #endif
47 
48 #include "RooRealVar.h"
49 #include "RooAddition.h"
50 #include "RooProduct.h"
51 #include "RooDataSet.h"
52 #include "RooProdPdf.h"
53 #include "RooFitResult.h"
54 #include "RooPoisson.h"
55 #include "RooGlobalFunc.h"
56 #include "RooCmdArg.h"
57 #include "RooWorkspace.h"
58 #include "RooMsgService.h"
59 #include "TTree.h"
60 #include <sstream>
61 
62 
63 
65 
66 
67 using namespace RooStats;
68 using namespace RooFit;
69 using namespace std;
70 
71 
72 ////////////////////////////////////////////////////////////////////////////////
73 /// constructor
74 
76 }
77 
78 ////////////////////////////////////////////////////////////////////////////////
79 /// destructor
80 
82 }
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 
87  Int_t nbins,
88  RooWorkspace* ws,
89  const char* pdfName, const char* muName) {
90 
91 // This method produces a PDF for N channels with uncorrelated background
92 // uncertainty. It relates the signal in each channel to a master signal strength times the
93 // expected signal in each channel.
94 //
95 // For the future, perhaps this method should be extended to include the efficiency terms automatically.
96 
97  using namespace RooFit;
98  using std::vector;
99 
100  TList likelihoodFactors;
101 
102  // Double_t MaxSigma = 8; // Needed to set ranges for varaibles.
103 
104  RooRealVar* masterSignal =
105  new RooRealVar(muName,"masterSignal",1., 0., 3.);
106 
107 
108  // loop over individual channels
109  for(Int_t i=0; i<nbins; ++i){
110  // need to name the variables dynamically, so please forgive the string manipulation and focus on values & ranges.
111  std::stringstream str;
112  str<<"_"<<i;
113  RooRealVar* expectedSignal =
114  new RooRealVar(("expected_s"+str.str()).c_str(),("expected_s"+str.str()).c_str(),sig[i], 0., 2*sig[i]);
115  expectedSignal->setConstant(kTRUE);
116 
117  RooProduct* s =
118  new RooProduct(("s"+str.str()).c_str(),("s"+str.str()).c_str(), RooArgSet(*masterSignal, *expectedSignal));
119 
120  RooRealVar* b =
121  new RooRealVar(("b"+str.str()).c_str(),("b"+str.str()).c_str(), .5, 0.,1.);
122  RooRealVar* tau =
123  new RooRealVar(("tau"+str.str()).c_str(),("tau"+str.str()).c_str(), .5, 0., 1.);
124  tau->setConstant(kTRUE);
125 
126  RooAddition* splusb =
127  new RooAddition(("splusb"+str.str()).c_str(),("s"+str.str()+"+"+"b"+str.str()).c_str(),
128  RooArgSet(*s,*b));
129  RooProduct* bTau =
130  new RooProduct(("bTau"+str.str()).c_str(),("b*tau"+str.str()).c_str(), RooArgSet(*b, *tau));
131  RooRealVar* x =
132  new RooRealVar(("x"+str.str()).c_str(),("x"+str.str()).c_str(), 0.5 , 0., 1.);
133  RooRealVar* y =
134  new RooRealVar(("y"+str.str()).c_str(),("y"+str.str()).c_str(), 0.5, 0., 1.);
135 
136 
137  RooPoisson* sigRegion =
138  new RooPoisson(("sigRegion"+str.str()).c_str(),("sigRegion"+str.str()).c_str(), *x,*splusb);
139 
140  //LM: need to set noRounding since y can take non integer values
141  RooPoisson* sideband =
142  new RooPoisson(("sideband"+str.str()).c_str(),("sideband"+str.str()).c_str(), *y,*bTau,true);
143 
144  likelihoodFactors.Add(sigRegion);
145  likelihoodFactors.Add(sideband);
146 
147  }
148 
149  RooArgSet likelihoodFactorSet(likelihoodFactors);
150  RooProdPdf joint(pdfName,"joint", likelihoodFactorSet );
151  // joint.printCompactTree();
152 
153  // add this PDF to workspace.
154  // Need to do import into workspace now to get all the structure imported as well.
155  // Just returning the WS will loose the rest of the structure b/c it will go out of scope
157  ws->import(joint);
159 }
160 
161 ////////////////////////////////////////////////////////////////////////////////
162 
164  Double_t* back,
165  Double_t* back_syst,
166  Int_t nbins,
167  RooWorkspace* ws, const char* dsName) {
168  // Arguements are an array of expected signal, expected background, and relative
169  // background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
170 
171  std::vector<Double_t> mainMeas(nbins);
172 
173  // loop over channels
174  for(Int_t i=0; i<nbins; ++i){
175  mainMeas[i] = sig[i] + back[i];
176  }
177  return AddData(&mainMeas[0], back, back_syst, nbins, ws, dsName);
178 }
179 
180 ////////////////////////////////////////////////////////////////////////////////
181 
183  Double_t* backExp,
184  Double_t* tau,
185  Int_t nbins,
186  RooWorkspace* ws, const char* dsName) {
187  // Arguements are an array of expected signal, expected background, and relative
188  // ratio of background expected in the sideband to that expected in signal region, and the number of channels.
189 
190  std::vector<Double_t> mainMeas(nbins);
191  std::vector<Double_t> sideband(nbins);
192  for(Int_t i=0; i<nbins; ++i){
193  mainMeas[i] = sigExp[i] + backExp[i];
194  sideband[i] = backExp[i]*tau[i];
195  }
196  return AddDataWithSideband(&mainMeas[0], &sideband[0], tau, nbins, ws, dsName);
197 
198 }
199 
200 ////////////////////////////////////////////////////////////////////////////////
201 /// need to be careful here that the range of observable in the dataset is consistent with the one in the workspace
202 /// don't rescale unless necessary. If it is necessary, then rescale by x10 or a defined maximum.
203 
205  return SafeObservableCreation(ws, varName, value, 10.*value);
206 
207 }
208 
209 ////////////////////////////////////////////////////////////////////////////////
210 /// need to be careful here that the range of observable in the dataset is consistent with the one in the workspace
211 /// don't rescale unless necessary. If it is necessary, then rescale by x10 or a defined maximum.
212 
214  Double_t value, Double_t maximum) {
215  RooRealVar* x = ws->var( varName );
216  if( !x )
217  x = new RooRealVar(varName, varName, value, 0, maximum );
218  if( x->getMax() < value )
219  x->setMax( max(x->getMax(), 10*value ) );
220  x->setVal( value );
221 
222  return x;
223 }
224 
225 
226 ////////////////////////////////////////////////////////////////////////////////
227 
229  Double_t* back,
230  Double_t* back_syst,
231  Int_t nbins,
232  RooWorkspace* ws, const char* dsName) {
233  // Arguments are an array of results from a main measurement, a measured background,
234  // and relative background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
235 
236  using namespace RooFit;
237  using std::vector;
238 
239  Double_t MaxSigma = 8; // Needed to set ranges for varaibles.
240 
241  TList observablesCollection;
242 
243  TTree* tree = new TTree();
244  std::vector<Double_t> xForTree(nbins);
245  std::vector<Double_t> yForTree(nbins);
246 
247  // loop over channels
248  for(Int_t i=0; i<nbins; ++i){
249  std::stringstream str;
250  str<<"_"<<i;
251 
252  //Double_t _tau = 1./back[i]/back_syst[i]/back_syst[i];
253  // LM: compute tau correctly for the Gamma distribution : mode = tau*b and variance is (tau*b+1)
254  Double_t err = back_syst[i];
255  Double_t _tau = (1.0 + sqrt(1 + 4 * err * err))/ (2. * err * err)/ back[i];
256 
257  RooRealVar* tau = SafeObservableCreation(ws, ("tau"+str.str()).c_str(), _tau );
258 
259  oocoutW(ws,ObjectHandling) << "NumberCountingPdfFactory: changed value of " << tau->GetName() << " to " << tau->getVal() <<
260  " to be consistent with background and its uncertainty. " <<
261  " Also stored these values of tau into workspace with name . " << (string(tau->GetName())+string(dsName)).c_str() <<
262  " if you test with a different dataset, you should adjust tau appropriately.\n"<< endl;
264  ws->import(*((RooRealVar*) tau->clone( (string(tau->GetName())+string(dsName)).c_str() ) ) );
266 
267  // need to be careful
268  RooRealVar* x = SafeObservableCreation(ws, ("x"+str.str()).c_str(), mainMeas[i]);
269 
270  // need to be careful
271  RooRealVar* y = SafeObservableCreation(ws, ("y"+str.str()).c_str(), back[i]*_tau );
272 
273  observablesCollection.Add(x);
274  observablesCollection.Add(y);
275 
276  xForTree[i] = mainMeas[i];
277  yForTree[i] = back[i]*_tau;
278 
279  tree->Branch(("x"+str.str()).c_str(), &xForTree[i] ,("x"+str.str()+"/D").c_str());
280  tree->Branch(("y"+str.str()).c_str(), &yForTree[i] ,("y"+str.str()+"/D").c_str());
281 
282  ws->var(("b"+str.str()).c_str())->setMax( 1.2*back[i]+MaxSigma*(sqrt(back[i])+back[i]*back_syst[i]) );
283  ws->var(("b"+str.str()).c_str())->setVal( back[i] );
284 
285  }
286  tree->Fill();
287  // tree->Print();
288  // tree->Scan();
289 
290  RooArgList* observableList = new RooArgList(observablesCollection);
291 
292  // observableSet->Print();
293  // observableList->Print();
294 
295  RooDataSet* data = new RooDataSet(dsName,"Number Counting Data", tree, *observableList); // one experiment
296  // data->Scan();
297 
298 
299  // import hypothetical data
301  ws->import(*data);
303 
304 }
305 
306 ////////////////////////////////////////////////////////////////////////////////
307 
309  Double_t* sideband,
310  Double_t* tauForTree,
311  Int_t nbins,
312  RooWorkspace* ws, const char* dsName) {
313  // Arguements are an array of expected signal, expected background, and relative
314  // background uncertainty (eg. 0.1 for 10% uncertainty), and the number of channels.
315 
316  using namespace RooFit;
317  using std::vector;
318 
319  Double_t MaxSigma = 8; // Needed to set ranges for varaibles.
320 
321  TList observablesCollection;
322 
323  TTree* tree = new TTree();
324 
325  std::vector<Double_t> xForTree(nbins);
326  std::vector<Double_t> yForTree(nbins);
327 
328 
329  // loop over channels
330  for(Int_t i=0; i<nbins; ++i){
331  std::stringstream str;
332  str<<"_"<<i;
333 
334  Double_t _tau = tauForTree[i];
335  Double_t back_syst = 1./sqrt(sideband[i]);
336  Double_t back = (sideband[i]/_tau);
337 
338 
339  RooRealVar* tau = SafeObservableCreation(ws, ("tau"+str.str()).c_str(), _tau );
340 
341  oocoutW(ws,ObjectHandling) << "NumberCountingPdfFactory: changed value of " << tau->GetName() << " to " << tau->getVal() <<
342  " to be consistent with background and its uncertainty. " <<
343  " Also stored these values of tau into workspace with name . " << (string(tau->GetName())+string(dsName)).c_str() <<
344  " if you test with a different dataset, you should adjust tau appropriately.\n"<< endl;
346  ws->import(*((RooRealVar*) tau->clone( (string(tau->GetName())+string(dsName)).c_str() ) ) );
348 
349  // need to be careful
350  RooRealVar* x = SafeObservableCreation(ws, ("x"+str.str()).c_str(), mainMeas[i]);
351 
352  // need to be careful
353  RooRealVar* y = SafeObservableCreation(ws, ("y"+str.str()).c_str(), sideband[i] );
354 
355 
356  observablesCollection.Add(x);
357  observablesCollection.Add(y);
358 
359  xForTree[i] = mainMeas[i];
360  yForTree[i] = sideband[i];
361 
362  tree->Branch(("x"+str.str()).c_str(), &xForTree[i] ,("x"+str.str()+"/D").c_str());
363  tree->Branch(("y"+str.str()).c_str(), &yForTree[i] ,("y"+str.str()+"/D").c_str());
364 
365  ws->var(("b"+str.str()).c_str())->setMax( 1.2*back+MaxSigma*(sqrt(back)+back*back_syst) );
366  ws->var(("b"+str.str()).c_str())->setVal( back );
367 
368  }
369  tree->Fill();
370  // tree->Print();
371  // tree->Scan();
372 
373  RooArgList* observableList = new RooArgList(observablesCollection);
374 
375  // observableSet->Print();
376  // observableList->Print();
377 
378  RooDataSet* data = new RooDataSet(dsName,"Number Counting Data", tree, *observableList); // one experiment
379  // data->Scan();
380 
381 
382  // import hypothetical data
384  ws->import(*data);
386 
387 }
388 
389 
390 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
Poisson pdf.
Definition: RooPoisson.h:19
virtual Double_t getMax(const char *name=0) const
virtual TObject * clone(const char *newname) const
Definition: RooRealVar.h:48
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
int Int_t
Definition: RtypesCore.h:41
int nbins[3]
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
void AddExpData(Double_t *sigExp, Double_t *bkgExp, Double_t *db, Int_t nbins, RooWorkspace *ws, const char *dsName="ExpectedNumberCountingData")
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:418
void AddExpDataWithSideband(Double_t *sigExp, Double_t *bkgExp, Double_t *tau, Int_t nbins, RooWorkspace *ws, const char *dsName="NumberCountingData")
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t tau
Definition: TRolke.cxx:630
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
A doubly linked list.
Definition: TList.h:47
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
void setConstant(Bool_t value=kTRUE)
void setGlobalKillBelow(RooFit::MsgLevel level)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
RooProduct a RooAbsReal implementation that represent the product of a given set of other RooAbsReal ...
Definition: RooProduct.h:32
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:279
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
void AddDataWithSideband(Double_t *mainMeas, Double_t *sideband, Double_t *tau, Int_t nbins, RooWorkspace *ws, const char *dsName="ExpectedNumberCountingData")
#define oocoutW(o, a)
Definition: RooMsgService.h:47
void AddData(Double_t *mainMeas, Double_t *bkgMeas, Double_t *db, Int_t nbins, RooWorkspace *ws, const char *dsName="NumberCountingData")
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
virtual void Add(TObject *obj)
Definition: TList.h:81
RooRealVar * SafeObservableCreation(RooWorkspace *ws, const char *varName, Double_t value)
need to be careful here that the range of observable in the dataset is consistent with the one in the...
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Definition: tree.py:1
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets...
Definition: RooAddition.h:26
const Bool_t kTRUE
Definition: Rtypes.h:91
void AddModel(Double_t *sigExp, Int_t nchan, RooWorkspace *ws, const char *pdfName="CombinedPdf", const char *masterSignalName="masterSignal")
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:42