Logo ROOT  
Reference Guide
RooBinnedGenContext.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooBinnedGenContext.cxx
19\class RooBinnedGenContext
20\ingroup Roofitcore
21
22RooBinnedGenContext is an efficient implementation of the
23generator context specific for binned pdfs.
24**/
25
26#include "Riostream.h"
27
28
29#include "RooMsgService.h"
30#include "RooBinnedGenContext.h"
31#include "RooAbsPdf.h"
32#include "RooRealVar.h"
33#include "RooDataHist.h"
34#include "RooDataSet.h"
35#include "RooRandom.h"
36
37using namespace std;
38
40;
41
42
43////////////////////////////////////////////////////////////////////////////////
44/// Constructor
45
47 const RooDataSet *prototype, const RooArgSet* auxProto,
49 RooAbsGenContext(model,vars,prototype,auxProto,verbose)
50{
51 cxcoutI(Generation) << "RooBinnedGenContext::ctor() setting up event special generator context for sum p.d.f. " << model.GetName()
52 << " for generation of observable(s) " << vars ;
53 if (prototype) ccxcoutI(Generation) << " with prototype data for " << *prototype->get() ;
54 if (auxProto && auxProto->getSize()>0) ccxcoutI(Generation) << " with auxiliary prototypes " << *auxProto ;
55 ccxcoutI(Generation) << endl ;
56
57 // Constructor. Build an array of generator contexts for each product component PDF
59 _pdf = (RooAbsPdf*) _pdfSet->find(model.GetName()) ;
61
62 // Fix normalization set of this RooAddPdf
63 if (prototype)
64 {
65 RooArgSet coefNSet(vars) ;
66 coefNSet.add(*prototype->get()) ;
67 _pdf->fixAddCoefNormalization(coefNSet) ;
68 }
69
71 _vars = _pdf->getObservables(vars) ;
72
73 // If pdf has boundary definitions, follow those for the binning
74 RooFIter viter = _vars->fwdIterator() ;
75 RooAbsArg* var ;
76 while((var=viter.next())) {
77 RooRealVar* rvar = dynamic_cast<RooRealVar*>(var) ;
78 if (rvar) {
79 list<Double_t>* binb = model.binBoundaries(*rvar,rvar->getMin(),rvar->getMax()) ;
80 delete binb ;
81 }
82 }
83
84
85 // Create empty RooDataHist
86 _hist = new RooDataHist("genData","genData",*_vars) ;
87
89}
90
91
92////////////////////////////////////////////////////////////////////////////////
93/// Destructor. Delete all owned subgenerator contexts
94
96{
97 delete _vars ;
98 delete _pdfSet ;
99 delete _hist ;
100}
101
102
103
104////////////////////////////////////////////////////////////////////////////////
105/// Attach given set of variables to internal p.d.f. clone
106
108{
110}
111
112
113
114////////////////////////////////////////////////////////////////////////////////
115/// One-time initialization of generator contex. Attach theEvent
116/// to internal p.d.f clone and forward initialization call to
117/// the component generators
118
120{
121 _pdf->recursiveRedirectServers(theEvent) ;
122
123}
124
125
126////////////////////////////////////////////////////////////////////////////////
127
129{
130 _expectedData = flag ;
131}
132
133
134////////////////////////////////////////////////////////////////////////////////
135
137{
138 // Scale to number of events and introduce Poisson fluctuations
139 _hist->reset() ;
140
141 Double_t nEvents = nEvt ;
142
143 if (nEvents<=0) {
144 if (!_pdf->canBeExtended()) {
145 coutE(InputArguments) << "RooAbsPdf::generateBinned(" << GetName()
146 << ") ERROR: No event count provided and p.d.f does not provide expected number of events" << endl ;
147 return 0 ;
148 } else {
149 // Don't round in expectedData mode
150 if (_expectedData || extended) {
151 nEvents = _pdf->expectedEvents(_vars) ;
152 } else {
153 nEvents = Int_t(_pdf->expectedEvents(_vars)+0.5) ;
154 }
155 }
156 }
157
158 // Sample p.d.f. distribution
160
161 // Output container
162 RooRealVar weight("weight","weight",0,1e9) ;
163 RooArgSet tmp(*_vars) ;
164 tmp.add(weight) ;
165 RooDataSet* wudata = new RooDataSet("wu","wu",tmp,RooFit::WeightVar("weight")) ;
166
167 vector<int> histOut(_hist->numEntries()) ;
168 Double_t histMax(-1) ;
169 Int_t histOutSum(0) ;
170 for (int i=0 ; i<_hist->numEntries() ; i++) {
171 _hist->get(i) ;
172 if (_expectedData) {
173
174 // Expected data, multiply p.d.f by nEvents
175 Double_t w=_hist->weight()*nEvents ;
176 wudata->add(*_hist->get(),w) ;
177
178 } else if (extended) {
179
180 // Extended mode, set contents to Poisson(pdf*nEvents)
182 wudata->add(*_hist->get(),w) ;
183
184 } else {
185
186 // Regular mode, fill array of weights with Poisson(pdf*nEvents), but to not fill
187 // histogram yet.
188 if (_hist->weight()>histMax) {
189 histMax = _hist->weight() ;
190 }
191 histOut[i] = RooRandom::randomGenerator()->Poisson(_hist->weight()*nEvents) ;
192 histOutSum += histOut[i] ;
193 }
194 }
195
196
197 if (!_expectedData && !extended) {
198
199 // Second pass for regular mode - Trim/Extend dataset to exact number of entries
200
201 // Calculate difference between what is generated so far and what is requested
202 Int_t nEvtExtra = abs(Int_t(nEvents)-histOutSum) ;
203 Int_t wgt = (histOutSum>nEvents) ? -1 : 1 ;
204
205 // Perform simple binned accept/reject procedure to get to exact event count
206 while(nEvtExtra>0) {
207
209 _hist->get(ibinRand) ;
210 Double_t ranY = RooRandom::randomGenerator()->Uniform(histMax) ;
211
212 if (ranY<_hist->weight()) {
213 if (wgt==1) {
214 histOut[ibinRand]++ ;
215 } else {
216 // If weight is negative, prior bin content must be at least 1
217 if (histOut[ibinRand]>0) {
218 histOut[ibinRand]-- ;
219 } else {
220 continue ;
221 }
222 }
223 nEvtExtra-- ;
224 }
225 }
226
227 // Transfer working array to histogram
228 for (int i=0 ; i<_hist->numEntries() ; i++) {
229 _hist->get(i) ;
230 wudata->add(*_hist->get(),histOut[i]) ;
231 }
232
233 }
234
235 return wudata ;
236
237}
238
239
240
241////////////////////////////////////////////////////////////////////////////////
242/// this method is not implemented for this context
243
245{
246 assert(0) ;
247}
248
249
250
251////////////////////////////////////////////////////////////////////////////////
252/// Print the details of the context
253
255{
257 os << indent << "--- RooBinnedGenContext ---" << endl ;
258 os << indent << "Using PDF ";
260}
#define cxcoutI(a)
Definition: RooMsgService.h:85
#define coutE(a)
Definition: RooMsgService.h:33
#define ccxcoutI(a)
Definition: RooMsgService.h:86
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:375
static void indent(ostringstream &buf, int indent_level)
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:78
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition: RooAbsArg.h:318
void setOperMode(OperMode mode, Bool_t recurseADirty=kTRUE)
Set the operation mode of this node.
Definition: RooAbsArg.cxx:1876
Bool_t recursiveRedirectServers(const RooAbsCollection &newServerList, Bool_t mustReplaceAll=kFALSE, Bool_t nameChange=kFALSE, Bool_t recurseInNewSet=kTRUE)
Recursively replace all servers with the new servers in newSet.
Definition: RooAbsArg.cxx:1197
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Int_t getSize() const
Return the number of elements in the collection.
RooFIter fwdIterator() const
One-time forward iterator.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const override
Interface for multi-line printing.
RooArgSet _theEvent
Pointer to observable event being generated.
Bool_t canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:263
virtual Double_t expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Definition: RooAbsPdf.cxx:3231
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Double_t getMin(const char *name=0) const
Get minimum of currently defined range.
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &obs, Double_t xlo, Double_t xhi) const
Retrieve bin boundaries if this distribution is binned in obs.
virtual void fixAddCoefNormalization(const RooArgSet &addNormSet=RooArgSet(), Bool_t force=kTRUE)
Fix the interpretation of the coefficient of any RooAddPdf component in the expression tree headed by...
RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, Double_t scaleFactor, Bool_t correctForBinVolume=kFALSE, Bool_t showProgress=kFALSE) const
Fill a RooDataHist with values sampled from this function at the bin centers.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:180
RooBinnedGenContext is an efficient implementation of the generator context specific for binned pdfs.
RooAbsPdf * _pdf
Pointer to cloned p.d.f.
void setExpectedData(Bool_t) override
RooBinnedGenContext(const RooAbsPdf &model, const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t _verbose=kFALSE)
Constructor.
~RooBinnedGenContext() override
Destructor. Delete all owned subgenerator contexts.
void generateEvent(RooArgSet &theEvent, Int_t remaining) override
this method is not implemented for this context
const RooArgSet * _vars
RooDataHist * _hist
Histogram.
RooDataSet * generate(Double_t nEvents=0, Bool_t skipInit=kFALSE, Bool_t extendedMode=kFALSE) override
Generate the specified number of events with nEvents>0 and and return a dataset containing the genera...
void initGenerator(const RooArgSet &theEvent) override
One-time initialization of generator contex.
Bool_t _expectedData
Asimov?
void printMultiline(std::ostream &os, Int_t content, Bool_t verbose=kFALSE, TString indent="") const override
Print the details of the context.
void attach(const RooArgSet &params) override
Attach given set of variables to internal p.d.f. clone.
RooArgSet * _pdfSet
Set owned all nodes of internal clone of p.d.f.
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:45
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:110
Int_t numEntries() const override
Return the number of bins.
void reset() override
Reset all bin weights to zero.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:82
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add one ore more rows of data.
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
RooAbsArg * next()
Return next element or nullptr if at end.
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:51
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:402
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:672
virtual UInt_t Integer(UInt_t imax)
Returns a random integer uniformly distributed on the interval [ 0, imax-1 ].
Definition: TRandom.cxx:360
Basic string class.
Definition: TString.h:136
RooCmdArg WeightVar(const char *name, Bool_t reinterpretAsWeight=kFALSE)
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1739
@ Generation
Definition: RooGlobalFunc.h:63
@ InputArguments
Definition: RooGlobalFunc.h:64