Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooConvGenContext.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 RooConvGenContext.cxx
19\class RooConvGenContext
20\ingroup Roofitcore
21
22Efficient implementation of the generator context
23specific for RooAbsAnaConvPdf objects. The physics model is generated
24with a truth resolution model and the requested resolution model is generated
25separately as a PDF. The convolution variable of the physics model is
26subsequently explicitly smeared with the resolution model distribution.
27**/
28
29#include "RooMsgService.h"
30#include "RooErrorHandler.h"
31#include "RooConvGenContext.h"
32#include "RooAbsAnaConvPdf.h"
33#include "RooNumConvPdf.h"
34#include "RooFFTConvPdf.h"
35#include "RooProdPdf.h"
36#include "RooDataSet.h"
37#include "RooArgSet.h"
38#include "RooTruthModel.h"
39#include "Riostream.h"
40
41
42using std::cout, std::endl, std::ostream;
43
45
46////////////////////////////////////////////////////////////////////////////////
47/// Constructor for specialized generator context for analytical convolutions.
48///
49/// Builds a generator for the physics PDF convoluted with the truth model
50/// and a generator for the resolution model as PDF. Events are generated
51/// by sampling events from the p.d.f and smearings from the resolution model
52/// and adding these to obtain a distribution of events consistent with the
53/// convolution of these two. The advantage of this procedure is so that
54/// both p.d.f and resolution model can take advantage of any internal
55/// generators that may be defined.
56
58 const RooDataSet *prototype, const RooArgSet* auxProto, bool verbose) :
59 RooAbsGenContext(model,vars,prototype,auxProto,verbose)
60{
61 cxcoutI(Generation) << "RooConvGenContext::ctor() setting up special generator context for analytical convolution p.d.f. " << model.GetName()
62 << " for generation of observable(s) " << vars << endl ;
63
64 // Clone PDF and change model to internal truth model
65 _pdfCloneSet = std::make_unique<RooArgSet>();
66 RooArgSet(model).snapshot(*_pdfCloneSet, true);
67 if (!_pdfCloneSet) {
68 coutE(Generation) << "RooConvGenContext::RooConvGenContext(" << GetName() << ") Couldn't deep-clone PDF, abort," << endl ;
70 }
71
72 RooAbsAnaConvPdf* pdfClone = static_cast<RooAbsAnaConvPdf*>(_pdfCloneSet->find(model.GetName())) ;
73 RooTruthModel truthModel("truthModel","Truth resolution model",*pdfClone->convVar()) ;
74 pdfClone->changeModel(truthModel) ;
75 auto convV = dynamic_cast<RooRealVar*>(pdfClone->convVar());
76 if (!convV) {
77 throw std::runtime_error("RooConvGenContext only works with convolution variables of type RooRealVar.");
78 }
79 convV->removeRange();
80
81 // Create generator for physics X truth model
82 _pdfVars = std::unique_ptr<RooArgSet>{pdfClone->getObservables(&vars)};
83 _pdfGen.reset(pdfClone->genContext(*_pdfVars,prototype,auxProto,verbose));
84
85 // Clone resolution model and use as normal PDF
86 _modelCloneSet = std::make_unique<RooArgSet>();
87 RooArgSet(*model._convSet.at(0)).snapshot(*_modelCloneSet, true);
88 if (!_modelCloneSet) {
89 coutE(Generation) << "RooConvGenContext::RooConvGenContext(" << GetName() << ") Couldn't deep-clone resolution model, abort," << std::endl;
91 }
92 std::unique_ptr<RooResolutionModel> modelClone{static_cast<RooResolutionModel*>(_modelCloneSet->find(model._convSet.at(0)->GetName())->Clone("smearing"))};
93 modelClone->changeBasis(nullptr) ;
94 convV = dynamic_cast<RooRealVar*>(&modelClone->convVar());
95 if (!convV) {
96 throw std::runtime_error("RooConvGenContext only works with convolution variables of type RooRealVar.");
97 }
98 convV->removeRange();
99
100 // Create generator for resolution model as PDF
101 _modelVars = std::unique_ptr<RooArgSet>{modelClone->getObservables(&vars)};
102
103 _modelVars->add(modelClone->convVar()) ;
104 _convVarName = modelClone->convVar().GetName() ;
105 _modelGen.reset(modelClone->genContext(*_modelVars,prototype,auxProto,verbose));
106
107 _modelCloneSet->addOwned(std::move(modelClone));
108
109 if (prototype) {
110 _pdfVars->add(*prototype->get()) ;
111 _modelVars->add(*prototype->get()) ;
112 }
113
114 // WVE ADD FOR DEBUGGING
115 if (auxProto) {
116 _pdfVars->add(*auxProto) ;
117 _modelVars->add(*auxProto) ;
118 }
119}
120
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Constructor for specialized generator context for numerical convolutions.
125///
126/// Builds a generator for the physics PDF convoluted with the truth model
127/// and a generator for the resolution model as PDF. Events are generated
128/// by sampling events from the p.d.f and smearings from the resolution model
129/// and adding these to obtain a distribution of events consistent with the
130/// convolution of these two. The advantage of this procedure is so that
131/// both p.d.f and resolution model can take advantage of any internal
132/// generators that may be defined.
133
135 const RooDataSet *prototype, const RooArgSet* auxProto, bool verbose) :
136 RooAbsGenContext(model,vars,prototype,auxProto,verbose)
137{
138 cxcoutI(Generation) << "RooConvGenContext::ctor() setting up special generator context for numeric convolution p.d.f. " << model.GetName()
139 << " for generation of observable(s) " << vars << endl ;
140
141 // Create generator for physics X truth model
142 {
143 RooArgSet clonedPdfObservables;
144 model.conv().clonePdf().getObservables(&vars, clonedPdfObservables);
145 _pdfVarsOwned = std::make_unique<RooArgSet>();
146 clonedPdfObservables.snapshot(*_pdfVarsOwned, true);
147 }
148 _pdfVars = std::make_unique<RooArgSet>(*_pdfVarsOwned) ;
149 _pdfGen.reset(static_cast<RooAbsPdf&>(model.conv().clonePdf()).genContext(*_pdfVars,prototype,auxProto,verbose));
150
151 // Create generator for resolution model as PDF
152 {
153 RooArgSet clonedModelObservables;
154 model.conv().cloneModel().getObservables(&vars, clonedModelObservables);
155 _modelVarsOwned = std::make_unique<RooArgSet>();
156 clonedModelObservables.snapshot(*_modelVarsOwned, true);
157 }
158 _modelVars = std::make_unique<RooArgSet>(*_modelVarsOwned) ;
159 _convVarName = model.conv().cloneVar().GetName() ;
160 _modelGen.reset(static_cast<RooAbsPdf&>(model.conv().cloneModel()).genContext(*_modelVars,prototype,auxProto,verbose));
161 _modelCloneSet = std::make_unique<RooArgSet>();
162 _modelCloneSet->add(model.conv().cloneModel()) ;
163
164 if (prototype) {
165 _pdfVars->add(*prototype->get()) ;
166 _modelVars->add(*prototype->get()) ;
167 }
168}
169
170
171
172////////////////////////////////////////////////////////////////////////////////
173/// Constructor for specialized generator context for FFT numerical convolutions.
174///
175/// Builds a generator for the physics PDF convoluted with the truth model
176/// and a generator for the resolution model as PDF. Events are generated
177/// by sampling events from the p.d.f and smearings from the resolution model
178/// and adding these to obtain a distribution of events consistent with the
179/// convolution of these two. The advantage of this procedure is so that
180/// both p.d.f and resolution model can take advantage of any internal
181/// generators that may be defined.
182
184 const RooDataSet *prototype, const RooArgSet* auxProto, bool verbose) :
185 RooAbsGenContext(model,vars,prototype,auxProto,verbose)
186{
187 cxcoutI(Generation) << "RooConvGenContext::ctor() setting up special generator context for fft convolution p.d.f. " << model.GetName()
188 << " for generation of observable(s) " << vars << endl ;
189
190 _convVarName = model._x.arg().GetName() ;
191
192 // Create generator for physics model
193 _pdfCloneSet = std::make_unique<RooArgSet>();
194 RooArgSet(model._pdf1.arg()).snapshot(*_pdfCloneSet, true);
195 RooAbsPdf* pdfClone = static_cast<RooAbsPdf*>(_pdfCloneSet->find(model._pdf1.arg().GetName())) ;
196 RooRealVar* cvPdf = static_cast<RooRealVar*>(_pdfCloneSet->find(model._x.arg().GetName())) ;
197 cvPdf->removeRange() ;
198 RooArgSet tmp1;
199 pdfClone->getObservables(&vars, tmp1) ;
200 _pdfVarsOwned = std::make_unique<RooArgSet>();
201 tmp1.snapshot(*_pdfVarsOwned, true);
202 _pdfVars = std::make_unique<RooArgSet>(*_pdfVarsOwned) ;
203 _pdfGen.reset(pdfClone->genContext(*_pdfVars,prototype,auxProto,verbose));
204
205 // Create generator for resolution model
206 _modelCloneSet = std::make_unique<RooArgSet>();
207 RooArgSet(model._pdf2.arg()).snapshot(*_modelCloneSet, true);
208 RooAbsPdf* modelClone = static_cast<RooAbsPdf*>(_modelCloneSet->find(model._pdf2.arg().GetName())) ;
209 RooRealVar* cvModel = static_cast<RooRealVar*>(_modelCloneSet->find(model._x.arg().GetName())) ;
210 cvModel->removeRange() ;
211 RooArgSet tmp2;
212 modelClone->getObservables(&vars, tmp2) ;
213 _modelVarsOwned = std::make_unique<RooArgSet>();
214 tmp2.snapshot(*_modelVarsOwned, true);
215 _modelVars = std::make_unique<RooArgSet>(*_modelVarsOwned) ;
216 _modelGen.reset(modelClone->genContext(*_pdfVars,prototype,auxProto,verbose));
217
218 if (prototype) {
219 _pdfVars->add(*prototype->get()) ;
220 _modelVars->add(*prototype->get()) ;
221 }
222}
223
224
225////////////////////////////////////////////////////////////////////////////////
226/// Attach given set of arguments to internal clones of
227/// pdf and resolution model
228
230{
231 // Find convolution variable in input and output sets
232 auto* cvModel = static_cast<RooRealVar*>(_modelVars->find(_convVarName));
233 auto* cvPdf = static_cast<RooRealVar*>(_pdfVars->find(_convVarName));
234
235 // Replace all servers in _pdfVars and _modelVars with those in theEvent, except for the convolution variable
236 std::unique_ptr<RooArgSet> pdfCommon{static_cast<RooArgSet*>(args.selectCommon(*_pdfVars))};
237 pdfCommon->remove(*cvPdf,true,true) ;
238
239 std::unique_ptr<RooArgSet> modelCommon{static_cast<RooArgSet*>(args.selectCommon(*_modelVars))};
240 modelCommon->remove(*cvModel,true,true) ;
241
242 _pdfGen->attach(*pdfCommon) ;
243 _modelGen->attach(*modelCommon) ;
244}
245
246
247////////////////////////////////////////////////////////////////////////////////
248/// One-time initialization of generator context, attaches
249/// the context to the supplied event container
250
252{
253 // Find convolution variable in input and output sets
254 _cvModel = static_cast<RooRealVar*>(_modelVars->find(_convVarName)) ;
255 _cvPdf = static_cast<RooRealVar*>(_pdfVars->find(_convVarName)) ;
256 _cvOut = static_cast<RooRealVar*>(theEvent.find(_convVarName)) ;
257
258 // Replace all servers in _pdfVars and _modelVars with those in theEvent, except for the convolution variable
259 std::unique_ptr<RooArgSet> pdfCommon{static_cast<RooArgSet*>(theEvent.selectCommon(*_pdfVars))};
260 pdfCommon->remove(*_cvPdf,true,true) ;
261 _pdfVars->replace(*pdfCommon) ;
262
263 std::unique_ptr<RooArgSet> modelCommon{static_cast<RooArgSet*>(theEvent.selectCommon(*_modelVars))};
264 modelCommon->remove(*_cvModel,true,true) ;
265 _modelVars->replace(*modelCommon) ;
266
267 // Initialize component generators
268 _pdfGen->initGenerator(*_pdfVars) ;
269 _modelGen->initGenerator(*_modelVars) ;
270}
271
272
273
274////////////////////////////////////////////////////////////////////////////////
275/// Generate a single event
276
278{
279 while(true) {
280
281 // Generate pdf and model data
282 _modelGen->generateEvent(*_modelVars,remaining) ;
283 _pdfGen->generateEvent(*_pdfVars,remaining) ;
284
285 // Construct smeared convolution variable
286 double convValSmeared = _cvPdf->getVal() + _cvModel->getVal() ;
287 if (_cvOut->isValidReal(convValSmeared)) {
288 // Smeared value in acceptance range, transfer values to output set
289 theEvent.assign(*_modelVars) ;
290 theEvent.assign(*_pdfVars) ;
291 _cvOut->setVal(convValSmeared) ;
292 return ;
293 }
294 }
295}
296
297
298
299////////////////////////////////////////////////////////////////////////////////
300/// Set the traversal order for events in the prototype dataset
301/// The argument is a array of integers with a size identical
302/// to the number of events in the prototype dataset. Each element
303/// should contain an integer in the range 1-N.
304
306{
308 _modelGen->setProtoDataOrder(lut) ;
309 _pdfGen->setProtoDataOrder(lut) ;
310}
311
312
313////////////////////////////////////////////////////////////////////////////////
314/// Print the details of this generator context
315
316void RooConvGenContext::printMultiline(ostream &os, Int_t content, bool verbose, TString indent) const
317{
318 RooAbsGenContext::printMultiline(os,content,verbose,indent) ;
319 os << indent << "--- RooConvGenContext ---" << endl ;
320 os << indent << "List of component generators" << endl ;
321
322 TString indent2(indent) ;
323 indent2.Append(" ") ;
324
325 _modelGen->printMultiline(os,content,verbose,indent2);
326 _pdfGen->printMultiline(os,content,verbose,indent2);
327}
#define cxcoutI(a)
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:377
static void indent(ostringstream &buf, int indent_level)
Base class for PDFs that represent a physics model that can be analytically convolved with a resoluti...
virtual bool changeModel(const RooResolutionModel &newModel)
Change the current resolution model to newModel.
RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const override
Create a generator context for this p.d.f.
RooAbsRealLValue * convVar()
Retrieve the convolution variable.
RooListProxy _convSet
Set of (resModel (x) basisFunc) convolution objects.
RooFit::OwningPtr< RooArgSet > getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
bool selectCommon(const RooAbsCollection &refColl, RooAbsCollection &outColl) const
Create a subset of the current collection, consisting only of those elements that are contained as we...
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract base class for generator contexts of RooAbsPdf objects.
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Interface for multi-line printing.
virtual void setProtoDataOrder(Int_t *lut)
Set the traversal order of prototype data to that in the lookup tables passed as argument.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const
Interface function to create a generator context from a p.d.f.
bool isValidReal(double value, bool printError=false) const override
Check if given value is valid.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:178
Efficient implementation of the generator context specific for RooAbsAnaConvPdf objects.
std::unique_ptr< RooAbsGenContext > _modelGen
Resolution model generator context.
RooConvGenContext(const RooFFTConvPdf &model, const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool _verbose=false)
Constructor for specialized generator context for FFT numerical convolutions.
std::unique_ptr< RooArgSet > _pdfCloneSet
Owner of PDF clone.
void setProtoDataOrder(Int_t *lut) override
Set the traversal order for events in the prototype dataset The argument is a array of integers with ...
std::unique_ptr< RooArgSet > _pdfVars
Holder of PDF x truth event.
std::unique_ptr< RooArgSet > _modelVars
Holder of resModel event.
RooRealVar * _cvPdf
Convolution variable in PDFxTruth event.
std::unique_ptr< RooArgSet > _modelVarsOwned
Owning version of modelVars ;.
void attach(const RooArgSet &params) override
Attach given set of arguments to internal clones of pdf and resolution model.
std::unique_ptr< RooArgSet > _modelCloneSet
Owner of resModel clone.
RooRealVar * _cvOut
Convolution variable in output event.
RooRealVar * _cvModel
Convolution variable in resModel event.
void generateEvent(RooArgSet &theEvent, Int_t remaining) override
Generate a single event.
void initGenerator(const RooArgSet &theEvent) override
One-time initialization of generator context, attaches the context to the supplied event container.
std::unique_ptr< RooArgSet > _pdfVarsOwned
Owning version of pdfVars ;.
TString _convVarName
Name of convolution variable.
void printMultiline(std::ostream &os, Int_t content, bool verbose=false, TString indent="") const override
Print the details of this generator context.
std::unique_ptr< RooAbsGenContext > _pdfGen
Physics model generator context.
Container class to hold unbinned data.
Definition RooDataSet.h:57
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
static void softAbort()
Soft abort function that interrupts macro execution but doesn't kill ROOT.
PDF for the numerical (FFT) convolution of two PDFs.
RooRealProxy _pdf1
First input p.d.f.
RooRealProxy _x
Convolution observable.
RooRealProxy _pdf2
Second input p.d.f.
Numeric 1-dimensional convolution operator PDF.
RooNumConvolution & conv() const
RooAbsReal & cloneModel() const
RooRealVar & cloneVar() const
RooAbsReal & clonePdf() const
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
void removeRange(const char *name=nullptr)
Remove range limits for binning with given name. Empty name means default range.
RooResolutionModel is the base class for PDFs that represent a resolution model that can be convolute...
virtual void changeBasis(RooFormulaVar *basis)
Change the basis function we convolute with.
const T & arg() const
Return reference to object held in proxy.
Implements a RooResolution model that corresponds to a delta function.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Basic string class.
Definition TString.h:139
TString & Append(const char *cs)
Definition TString.h:572