Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooFFTConvPdf.cxx
Go to the documentation of this file.
1
2 /*****************************************************************************
3 * Project: RooFit *
4 * *
5 * Copyright (c) 2000-2005, Regents of the University of California *
6 * and Stanford University. All rights reserved. *
7 * *
8 * Redistribution and use in source and binary forms, *
9 * with or without modification, are permitted according to the terms *
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
11 *****************************************************************************/
12
13//////////////////////////////////////////////////////////////////////////////
14/// \class RooFFTConvPdf
15/// \ingroup Roofitcore
16///
17/// This class implements a generic one-dimensional numeric convolution of two PDFs,
18/// and can convolve any two RooAbsPdfs. The class exploits the convolution theorem
19/// \f[
20/// f(x) * g(x) \rightarrow F(k_i) \cdot G(k_i)
21/// \f]
22/// to calculate the convolution by calculating a Real->Complex FFT of both input PDFs,
23/// multiplying the complex coefficients and performing the reverse Complex->Real FFT
24/// to get the result in the input space. This class uses the ROOT FFT interface to
25/// the (free) FFTW3 package (www.fftw.org), and requires that your ROOT installation is
26/// compiled with the `fftw3=ON` (default). Instructions for manually installing fftw below.
27///
28/// Note that the performance in terms of speed and stability of RooFFTConvPdf is
29/// vastly superior to that of RooNumConvPdf.
30///
31/// An important feature of FFT convolutions is that the observable is assumed to be
32/// cyclical. This is correct for cyclical observables such as angles,
33/// but does not hold in general. For non-cyclical variables, wrap-around artifacts may be
34/// encountered, *e.g.* if the PDF is zero at xMin and non-zero at xMax. A rising tail may appear at xMin.
35/// This is inevitable when using FFTs. A distribution with 3 bins therefore looks like:
36/// ```
37/// ... 0 1 2 0 1 2 0 1 2 ...
38/// ```
39///
40/// Therefore, if bins 0 and 2 are not equal, the FFT sees a cyclical function with a step at the 2|0 boundary, which causes
41/// artifacts in Fourier space.
42///
43/// The spillover or discontinuity can be reduced or eliminated by
44/// introducing a buffer zone in the FFT calculation. If this feature is activated (on by default),
45/// the sampling array for the FFT calculation is extended in both directions,
46/// and padded with the lowest/highest bin.
47/// Example:
48/// ```
49/// original: -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5
50/// add buffer zones: U U -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5 O O
51/// rotate: 0 +1 +2 +3 +4 +5 O O U U -5 -4 -3 -2 -1
52/// ```
53/// The buffer bins are stripped away when the FFT output values
54/// are transferred back to the p.d.f cache. The default buffer size is 10% of the
55/// observable domain size, and can be changed with the `setBufferFraction()` member function.
56///
57/// The RooFFTConvPdf uses caching inherited from a RooAbsCachedPdf. If it is
58/// evaluated for a particular value of x, the FFT and convolution is calculated
59/// for all bins in the observable space for the given choice of parameters,
60/// which are also stored in the cache. Subsequent evaluations for different values of the convolution observable and
61/// identical parameters will be retrieved from the cache. If one or more
62/// of the parameters change, the cache will be updated, *i.e.*, a new FFT runs.
63///
64/// The sampling density of the FFT is controlled by the binning of the
65/// the convolution observable, which can be changed using RooRealVar::setBins(N).
66/// For good results, N should be large (>=1000). Additional interpolation
67/// between the bins may improve the result if coarse binnings are chosen. These can be
68/// activated in the constructor or by calling `setInterpolationOrder()`.
69/// For N >> 1000, interpolation will not substantially improve the accuracy.
70///
71/// Additionial information on caching can be displayed by monitoring
72/// the message stream with topic "Caching" at the INFO level, *i.e.*
73/// by calling `RooMsgService::instance().addStream(RooMsgService::INFO,Topic("Caching"))`
74/// to see these message on stdout.
75///
76/// Multi-dimensional convolutions are not supported at the moment.
77///
78/// ---
79///
80/// Installing an external version of FFTW on Linux and compiling ROOT to use it
81/// -------
82///
83/// You have two options:
84/// * **Recommended**: ROOT can automatically install FFTW for itself, see `builtin_fftw3` at https://root.cern.ch/building-root
85/// * Install FFTW and let ROOT discover it. `fftw3` is on by default (see https://root.cern.ch/building-root)
86///
87/// 1) Go to www.fftw.org and download the latest stable version (a .tar.gz file)
88///
89/// If you have root access to your machine and want to make a system installation of FFTW
90///
91/// 2) Untar fftw-XXX.tar.gz in /tmp, cd into the untarred directory
92/// and type './configure' followed by 'make install'.
93/// This will install fftw in /usr/local/bin,lib etc...
94///
95/// 3) Start from a source installation of ROOT. ROOT should discover it. See https://root.cern.ch/building-root
96///
97///
98/// If you do not have root access and want to make a private installation of FFTW
99///
100/// 2) Make a private install area for FFTW, e.g. /home/myself/fftw
101///
102/// 3) Untar fftw-XXX.tar.gz in /tmp, cd into the untarred directory
103/// and type './configure --prefix=/home/myself/fftw' followed by 'make install'.
104/// Substitute /home/myself/fftw with a directory of your choice. This
105/// procedure will install FFTW in the location designated by you
106///
107/// 4) Start from a source installation of ROOT.
108/// Look up and set the proper paths for ROOT to discover FFTW. See https://root.cern.ch/building-root
109///
110
111#include "RooFFTConvPdf.h"
112
113#include "RooAbsReal.h"
114#include "RooMsgService.h"
115#include "RooDataHist.h"
116#include "RooHistPdf.h"
117#include "RooRealVar.h"
118#include "RooGenContext.h"
119#include "RooConvGenContext.h"
120#include "RooBinning.h"
121#include "RooLinearVar.h"
122#include "RooCustomizer.h"
123#include "RooGlobalFunc.h"
124#include "RooConstVar.h"
125#include "RooUniformBinning.h"
126
127#include "TClass.h"
128#include "TComplex.h"
129#include "TVirtualFFT.h"
130
131#include <iostream>
132#include <stdexcept>
133
134using namespace std;
135
137
138
139////////////////////////////////////////////////////////////////////////////////
140/// Constructor for numerical (FFT) convolution of PDFs.
141/// \param[in] name Name of this PDF
142/// \param[in] title Title for plotting this PDF
143/// \param[in] convVar Observable to convolve the PDFs in \attention Use a high number of bins (>= 1000) for good accuracy.
144/// \param[in] pdf1 First PDF to be convolved
145/// \param[in] pdf2 Second PDF to be convolved
146/// \param[in] ipOrder Order for interpolation between bins (since FFT is discrete)
147/// The binning used for the FFT sampling is controlled by the binning named "cache" in the convolution observable `convVar`.
148/// If such a binning is not set, the same number of bins as for `convVar` will be used.
149
150RooFFTConvPdf::RooFFTConvPdf(const char *name, const char *title, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder) :
151 RooAbsCachedPdf(name,title,ipOrder),
152 _x("!x","Convolution Variable",this,convVar),
153 _xprime("!xprime","External Convolution Variable",this,false),
154 _pdf1("!pdf1","pdf1",this,pdf1,false),
155 _pdf2("!pdf2","pdf2",this,pdf2,false),
156 _params("!params","effective parameters",this),
157 _bufFrac(0.1),
158 _bufStrat(Extend),
159 _shift1(0),
160 _shift2(0),
161 _cacheObs("!cacheObs","Cached observables",this,false,false)
162{
163 prepareFFTBinning(convVar);
164
165 _shift2 = (convVar.getMax("cache")+convVar.getMin("cache"))/2 ;
166
167 calcParams() ;
168
169}
170
171////////////////////////////////////////////////////////////////////////////////
172/// \copydoc RooFFTConvPdf(const char*, const char*, RooRealVar&, RooAbsPdf&, RooAbsPdf&, Int_t)
173/// \param[in] pdfConvVar If the variable used for convolution is a PDF, itself, pass the PDF here, and pass the convolution variable to
174/// `convVar`. See also rf210_angularconv.C in the <a href="https://root.cern.ch/root/html/tutorials/roofit/index.html.">roofit tutorials</a>
175
176RooFFTConvPdf::RooFFTConvPdf(const char *name, const char *title, RooAbsReal& pdfConvVar, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder) :
177 RooAbsCachedPdf(name,title,ipOrder),
178 _x("!x","Convolution Variable",this,convVar,false,false),
179 _xprime("!xprime","External Convolution Variable",this,pdfConvVar),
180 _pdf1("!pdf1","pdf1",this,pdf1,false),
181 _pdf2("!pdf2","pdf2",this,pdf2,false),
182 _params("!params","effective parameters",this),
183 _bufFrac(0.1),
184 _bufStrat(Extend),
185 _shift1(0),
186 _shift2(0),
187 _cacheObs("!cacheObs","Cached observables",this,false,false)
188{
189 prepareFFTBinning(convVar);
190
191 _shift2 = (convVar.getMax("cache")+convVar.getMin("cache"))/2 ;
192
193 calcParams() ;
194}
195
196
197
198////////////////////////////////////////////////////////////////////////////////
199/// Copy constructor
200
202 RooAbsCachedPdf(other,name),
203 _x("!x",this,other._x),
204 _xprime("!xprime",this,other._xprime),
205 _pdf1("!pdf1",this,other._pdf1),
206 _pdf2("!pdf2",this,other._pdf2),
207 _params("!params",this,other._params),
208 _bufFrac(other._bufFrac),
209 _bufStrat(other._bufStrat),
210 _shift1(other._shift1),
211 _shift2(other._shift2),
212 _cacheObs("!cacheObs",this,other._cacheObs)
213 {
214 }
215
216
217
218////////////////////////////////////////////////////////////////////////////////
219/// Destructor
220
222{
223}
224
225
226////////////////////////////////////////////////////////////////////////////////
227/// Try to improve the binning and inform user if possible.
228/// With a 10% buffer fraction, 930 raw bins yield 1024 FFT bins,
229/// a sweet spot for the speed of FFTW.
230
232 if (!convVar.hasBinning("cache")) {
233 const RooAbsBinning& varBinning = convVar.getBinning();
234 const int optimal = static_cast<Int_t>(1024/(1.+_bufFrac));
235
236 //Can improve precision if binning is uniform
237 if (varBinning.numBins() < optimal && varBinning.isUniform()) {
238 coutI(Caching) << "Changing internal binning of variable '" << convVar.GetName()
239 << "' in FFT '" << fName << "'"
240 << " from " << varBinning.numBins()
241 << " to " << optimal << " to improve the precision of the numerical FFT."
242 << " This can be done manually by setting an additional binning named 'cache'." << std::endl;
243 convVar.setBinning(RooUniformBinning(varBinning.lowBound(), varBinning.highBound(), optimal, "cache"), "cache");
244 } else {
245 coutE(Caching) << "The internal binning of variable " << convVar.GetName()
246 << " is not uniform. The numerical FFT will likely yield wrong results." << std::endl;
247 convVar.setBinning(varBinning, "cache");
248 }
249 }
250}
251
252
253////////////////////////////////////////////////////////////////////////////////
254/// Return base name component for cache components in this case 'PDF1_CONV_PDF2'
255
257{
258 static TString name ;
259 name = _pdf1.arg().GetName() ;
260 name.Append("_CONV_") ;
261 name.Append(_pdf2.arg().GetName()) ;
262 return name.Data() ;
263}
264
265
266
267
268////////////////////////////////////////////////////////////////////////////////
269/// Return specialized cache subclass for FFT calculations
270
272{
273 return new FFTCacheElem(*this,nset) ;
274}
275
276
277
278
279////////////////////////////////////////////////////////////////////////////////
280/// Clone input pdf and attach to dataset
281
283 PdfCacheElem(self,nsetIn)
284{
285 RooAbsPdf* clonePdf1 = (RooAbsPdf*) self._pdf1.arg().cloneTree() ;
286 RooAbsPdf* clonePdf2 = (RooAbsPdf*) self._pdf2.arg().cloneTree() ;
287 clonePdf1->attachDataSet(*hist()) ;
288 clonePdf2->attachDataSet(*hist()) ;
289
290 // Shift observable
291 RooRealVar* convObs = (RooRealVar*) hist()->get()->find(self._x.arg().GetName()) ;
292
293 // Install FFT reference range
294 string refName = Form("refrange_fft_%s",self.GetName()) ;
295 convObs->setRange(refName.c_str(),convObs->getMin(),convObs->getMax()) ;
296
297 if (self._shift1!=0) {
298 RooLinearVar* shiftObs1 = new RooLinearVar(Form("%s_shifted_FFTBuffer1",convObs->GetName()),"shiftObs1",
299 *convObs,RooFit::RooConst(1),RooFit::RooConst(-1*self._shift1)) ;
300
301 RooArgSet clonedBranches1 ;
302 RooCustomizer cust(*clonePdf1,"fft") ;
303 cust.replaceArg(*convObs,*shiftObs1) ;
304
305 pdf1Clone.reset(static_cast<RooAbsPdf*>(cust.build()));
306
307 pdf1Clone->addOwnedComponents(*shiftObs1) ;
308 pdf1Clone->addOwnedComponents(*clonePdf1) ;
309
310 } else {
311 pdf1Clone.reset(clonePdf1) ;
312 }
313
314 if (self._shift2!=0) {
315 RooLinearVar* shiftObs2 = new RooLinearVar(Form("%s_shifted_FFTBuffer2",convObs->GetName()),"shiftObs2",
316 *convObs,RooFit::RooConst(1),RooFit::RooConst(-1*self._shift2)) ;
317
318 RooArgSet clonedBranches2 ;
319 RooCustomizer cust(*clonePdf2,"fft") ;
320 cust.replaceArg(*convObs,*shiftObs2) ;
321
322 pdf1Clone->addOwnedComponents(*shiftObs2) ;
323 pdf1Clone->addOwnedComponents(*clonePdf2) ;
324
325 pdf2Clone.reset(static_cast<RooAbsPdf*>(cust.build()));
326
327 } else {
328 pdf2Clone.reset(clonePdf2) ;
329 }
330
331
332 // Attach cloned pdf to all original parameters of self
333 RooArgSet convObsSet{*convObs};
334 RooArgSet fftParams;
335 self.getParameters(&convObsSet, fftParams) ;
336
337 // Remove all cache histogram from fftParams as these
338 // observable need to remain attached to the histogram
339 fftParams.remove(*hist()->get(),true,true) ;
340
341 pdf1Clone->recursiveRedirectServers(fftParams) ;
342 pdf2Clone->recursiveRedirectServers(fftParams) ;
343 pdf1Clone->fixAddCoefRange(refName.c_str(), true) ;
344 pdf2Clone->fixAddCoefRange(refName.c_str(), true) ;
345
346 // Ensure that coefficients for Add PDFs are only interpreted with respect to the convolution observable
347 RooArgSet convSet(self._x.arg());
348 pdf1Clone->fixAddCoefNormalization(convSet, true);
349 pdf2Clone->fixAddCoefNormalization(convSet, true);
350
351 // Save copy of original histX binning and make alternate binning
352 // for extended range scanning
353
354 const Int_t N = convObs->numBins();
355 if (N < 900) {
356 oocoutW(&self, Eval) << "The FFT convolution '" << self.GetName() << "' will run with " << N
357 << " bins. A decent accuracy for difficult convolutions is typically only reached with n >= 1000. Suggest to increase the number"
358 " of bins of the observable '" << convObs->GetName() << "'." << std::endl;
359 }
360 Int_t Nbuf = static_cast<Int_t>((N*self.bufferFraction())/2 + 0.5) ;
361 double obw = (convObs->getMax() - convObs->getMin())/N ;
362 Int_t N2 = N+2*Nbuf ;
363
364 scanBinning = std::make_unique<RooUniformBinning>(convObs->getMin()-Nbuf*obw,convObs->getMax()+Nbuf*obw,N2);
365 histBinning.reset(convObs->getBinning().clone());
366
367 // Deactivate dirty state propagation on datahist observables
368 // and set all nodes on both pdfs to operMode AlwaysDirty
369 hist()->setDirtyProp(false) ;
370 convObs->setOperMode(ADirty,true) ;
371}
372
373
374////////////////////////////////////////////////////////////////////////////////
375/// Suffix for cache histogram (added in addition to suffix for cache name)
376
378{
379 return TString(Form("_BufFrac%3.1f_BufStrat%d",_bufFrac,_bufStrat)) ;
380}
381
382
383
384////////////////////////////////////////////////////////////////////////////////
385/// Returns all RooAbsArg objects contained in the cache element
386
388{
389 RooArgList ret(PdfCacheElem::containedArgs(a)) ;
390
391 ret.add(*pdf1Clone) ;
392 ret.add(*pdf2Clone) ;
393 if (pdf1Clone->ownedComponents()) {
394 ret.add(*pdf1Clone->ownedComponents()) ;
395 }
396 if (pdf2Clone->ownedComponents()) {
397 ret.add(*pdf2Clone->ownedComponents()) ;
398 }
399
400 return ret ;
401}
402
403
404
405////////////////////////////////////////////////////////////////////////////////
406/// Fill the contents of the cache the FFT convolution output
407
409{
410 RooDataHist& cacheHist = *cache.hist() ;
411
412 ((FFTCacheElem&)cache).pdf1Clone->setOperMode(ADirty,true) ;
413 ((FFTCacheElem&)cache).pdf2Clone->setOperMode(ADirty,true) ;
414
415 // Determine if there other observables than the convolution observable in the cache
416 RooArgSet otherObs ;
417 RooArgSet(*cacheHist.get()).snapshot(otherObs) ;
418
419 RooAbsArg* histArg = otherObs.find(_x.arg().GetName()) ;
420 if (histArg) {
421 otherObs.remove(*histArg,true,true) ;
422 }
423
424 //cout << "RooFFTConvPdf::fillCacheObject() otherObs = " << otherObs << endl ;
425
426 // Handle trivial scenario -- no other observables
427 if (otherObs.empty()) {
429 return ;
430 }
431
432 // Handle cases where there are other cache slices
433 // Iterator over available slice positions and fill each
434
435 // Determine number of bins for each slice position observable
436 Int_t n = otherObs.getSize() ;
437 std::vector<int> binCur(n+1);
438 std::vector<int> binMax(n+1);
439 Int_t curObs = 0 ;
440
441 std::vector<RooAbsLValue*> obsLV(n);
442 Int_t i(0) ;
443 for (auto const& arg : otherObs) {
444 RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(arg) ;
445 obsLV[i] = lvarg ;
446 binCur[i] = 0 ;
447 // coverity[FORWARD_NULL]
448 binMax[i] = lvarg->numBins(binningName())-1 ;
449 i++ ;
450 }
451
452 bool loop(true) ;
453 while(loop) {
454 // Set current slice position
455 for (Int_t j=0 ; j<n ; j++) { obsLV[j]->setBin(binCur[j],binningName()) ; }
456
457// cout << "filling slice: bin of obsLV[0] = " << obsLV[0]->getBin() << endl ;
458
459 // Fill current slice
460 fillCacheSlice((FFTCacheElem&)cache,otherObs) ;
461
462 // Determine which iterator to increment
463 while(binCur[curObs]==binMax[curObs]) {
464
465 // Reset current iterator and consider next iterator ;
466 binCur[curObs]=0 ;
467 curObs++ ;
468
469 // master termination condition
470 if (curObs==n) {
471 loop=false ;
472 break ;
473 }
474 }
475
476 // Increment current iterator
477 binCur[curObs]++ ;
478 curObs=0 ;
479
480 }
481
482}
483
484
485////////////////////////////////////////////////////////////////////////////////
486/// Fill a slice of cachePdf with the output of the FFT convolution calculation
487
488void RooFFTConvPdf::fillCacheSlice(FFTCacheElem& aux, const RooArgSet& slicePos) const
489{
490 // Extract histogram that is the basis of the RooHistPdf
491 RooDataHist& cacheHist = *aux.hist() ;
492
493 // Sample array of input points from both pdfs
494 // Note that returned arrays have optional buffers zones below and above range ends
495 // to reduce cyclical effects and have been cyclically rotated so that bin containing
496 // zero value is at position zero. Example:
497 //
498 // original: -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5
499 // add buffer zones: U U -5 -4 -3 -2 -1 0 +1 +2 +3 +4 +5 O O
500 // rotate: 0 +1 +2 +3 +4 +5 O O U U -5 -4 -3 -2 -1
501 //
502 //
503
504 Int_t N,N2,binShift1,binShift2 ;
505
506 RooRealVar* histX = (RooRealVar*) cacheHist.get()->find(_x.arg().GetName()) ;
507 if (_bufStrat==Extend) histX->setBinning(*aux.scanBinning) ;
508 std::vector<double> input1 = scanPdf((RooRealVar&)_x.arg(),*aux.pdf1Clone,cacheHist,slicePos,N,N2,binShift1,_shift1) ;
509 std::vector<double> input2 = scanPdf((RooRealVar&)_x.arg(),*aux.pdf2Clone,cacheHist,slicePos,N,N2,binShift2,_shift2) ;
510 if (_bufStrat==Extend) histX->setBinning(*aux.histBinning) ;
511
512
513
514
515 // Retrieve previously defined FFT transformation plans
516 if (!aux.fftr2c1) {
517 aux.fftr2c1.reset(TVirtualFFT::FFT(1, &N2, "R2CK"));
518 aux.fftr2c2.reset(TVirtualFFT::FFT(1, &N2, "R2CK"));
519 aux.fftc2r.reset(TVirtualFFT::FFT(1, &N2, "C2RK"));
520
521 if (aux.fftr2c1 == nullptr || aux.fftr2c2 == nullptr || aux.fftc2r == nullptr) {
522 coutF(Eval) << "RooFFTConvPdf::fillCacheSlice(" << GetName() << "Cannot get a handle to fftw. Maybe ROOT was built without it?" << std::endl;
523 throw std::runtime_error("Cannot get a handle to fftw.");
524 }
525 }
526
527 // Real->Complex FFT Transform on p.d.f. 1 sampling
528 aux.fftr2c1->SetPoints(input1.data());
529 aux.fftr2c1->Transform();
530
531 // Real->Complex FFT Transform on p.d.f 2 sampling
532 aux.fftr2c2->SetPoints(input2.data());
533 aux.fftr2c2->Transform();
534
535 // Loop over first half +1 of complex output results, multiply
536 // and set as input of reverse transform
537 for (Int_t i=0 ; i<N2/2+1 ; i++) {
538 double re1,re2,im1,im2 ;
539 aux.fftr2c1->GetPointComplex(i,re1,im1) ;
540 aux.fftr2c2->GetPointComplex(i,re2,im2) ;
541 double re = re1*re2 - im1*im2 ;
542 double im = re1*im2 + re2*im1 ;
543 TComplex t(re,im) ;
544 aux.fftc2r->SetPointComplex(i,t) ;
545 }
546
547 // Reverse Complex->Real FFT transform product
548 aux.fftc2r->Transform() ;
549
550 Int_t totalShift = binShift1 + (N2-N)/2 ;
551
552 // Store FFT result in cache
553
554 std::unique_ptr<TIterator> iter{const_cast<RooDataHist&>(cacheHist).sliceIterator(const_cast<RooAbsReal&>(_x.arg()),slicePos)};
555 for (Int_t i =0 ; i<N ; i++) {
556
557 // Cyclically shift array back so that bin containing zero is back in zeroBin
558 Int_t j = i + totalShift ;
559 while (j<0) j+= N2 ;
560 while (j>=N2) j-= N2 ;
561
562 iter->Next() ;
563 cacheHist.set(aux.fftc2r->GetPointReal(j)) ;
564 }
565
566 // cacheHist.dump2() ;
567}
568
569
570////////////////////////////////////////////////////////////////////////////////
571/// Scan the values of 'pdf' in observable 'obs' using the bin values stored in 'hist' at slice position 'slicePos'
572/// N is filled with the number of bins defined in hist, N2 is filled with N plus the number of buffer bins
573/// The return value is an array of doubles of length N2 with the sampled values. The caller takes ownership
574/// of the array
575
576std::vector<double> RooFFTConvPdf::scanPdf(RooRealVar& obs, RooAbsPdf& pdf, const RooDataHist& hist, const RooArgSet& slicePos,
577 Int_t& N, Int_t& N2, Int_t& zeroBin, double shift) const
578{
579
580 RooRealVar* histX = (RooRealVar*) hist.get()->find(obs.GetName()) ;
581
582 // Calculate number of buffer bins on each size to avoid cyclical flow
583 N = histX->numBins(binningName()) ;
584 Int_t Nbuf = static_cast<Int_t>((N*bufferFraction())/2 + 0.5) ;
585 N2 = N+2*Nbuf ;
586
587
588 // Allocate array of sampling size plus optional buffer zones
589 std::vector<double> array(N2);
590
591 // Set position of non-convolution observable to that of the cache slice that were are processing now
592 hist.get(slicePos) ;
593
594 // Find bin ID that contains zero value
595 zeroBin = 0 ;
596 if (histX->getMax()>=0 && histX->getMin()<=0) {
597 zeroBin = histX->getBinning().binNumber(0) ;
598 } else if (histX->getMin()>0) {
599 double bw = (histX->getMax() - histX->getMin())/N2 ;
600 zeroBin = Int_t(-histX->getMin()/bw) ;
601 } else {
602 double bw = (histX->getMax() - histX->getMin())/N2 ;
603 zeroBin = Int_t(-1*histX->getMax()/bw) ;
604 }
605
606 Int_t binShift = Int_t((N2* shift) / (histX->getMax()-histX->getMin())) ;
607
608 zeroBin += binShift ;
609 while(zeroBin>=N2) zeroBin-= N2 ;
610 while(zeroBin<0) zeroBin+= N2 ;
611
612 // First scan hist into temp array
613 std::vector<double> tmp(N2);
614 Int_t k(0) ;
615 switch(_bufStrat) {
616
617 case Extend:
618 // Sample entire extended range (N2 samples)
619 for (k=0 ; k<N2 ; k++) {
620 histX->setBin(k) ;
621 tmp[k] = pdf.getVal(hist.get()) ;
622 }
623 break ;
624
625 case Flat:
626 // Sample original range (N samples) and fill lower and upper buffer
627 // bins with p.d.f. value at respective boundary
628 {
629 histX->setBin(0) ;
630 double val = pdf.getVal(hist.get()) ;
631 for (k=0 ; k<Nbuf ; k++) {
632 tmp[k] = val ;
633 }
634 for (k=0 ; k<N ; k++) {
635 histX->setBin(k) ;
636 tmp[k+Nbuf] = pdf.getVal(hist.get()) ;
637 }
638 histX->setBin(N-1) ;
639 val = pdf.getVal(hist.get()) ;
640 for (k=0 ; k<Nbuf ; k++) {
641 tmp[N+Nbuf+k] = val ;
642 }
643 }
644 break ;
645
646 case Mirror:
647 // Sample original range (N samples) and fill lower and upper buffer
648 // bins with mirror image of sampled range
649 for (k=0 ; k<N ; k++) {
650 histX->setBin(k) ;
651 tmp[k+Nbuf] = pdf.getVal(hist.get()) ;
652 }
653 for (k=1 ; k<=Nbuf ; k++) {
654 histX->setBin(k) ;
655 tmp[Nbuf-k] = pdf.getVal(hist.get()) ;
656 histX->setBin(N-k) ;
657 tmp[Nbuf+N+k-1] = pdf.getVal(hist.get()) ;
658 }
659 break ;
660 }
661
662 // Scan function and store values in array
663 for (Int_t i=0 ; i<N2 ; i++) {
664 // Cyclically shift writing location by zero bin position
665 Int_t j = i - (zeroBin) ;
666 if (j<0) j+= N2 ;
667 if (j>=N2) j-= N2 ;
668 array[i] = tmp[j] ;
669 }
670
671 return array ;
672}
673
674
675
676////////////////////////////////////////////////////////////////////////////////
677/// Return the observables to be cached given the normalization set nset.
678///
679/// If the cache observable is in nset then this is
680/// - the convolution observable plus
681/// - any member of nset that is either a RooCategory,
682/// - or was previously specified through setCacheObservables().
683///
684/// In case the cache observable is *not* in nset, then it is
685/// - the convolution observable plus
686/// - all member of nset that are observables of this p.d.f.
687///
688
690{
691 // Get complete list of observables
692 auto obs1 = new RooArgSet{};
693 RooArgSet obs2;
694 _pdf1.arg().getObservables(&nset, *obs1) ;
695 _pdf2.arg().getObservables(&nset, obs2) ;
696 obs1->add(obs2, true) ;
697
698 // Check if convolution observable is in nset
699 if (nset.contains(_x.arg())) {
700
701 // Now strip out all non-category observables
702 RooArgSet killList ;
703 for(auto const& arg : *obs1) {
704 if (arg->IsA()->InheritsFrom(RooAbsReal::Class()) && !_cacheObs.find(arg->GetName())) {
705 killList.add(*arg) ;
706 }
707 }
708 obs1->remove(killList) ;
709
710 // And add back the convolution observables
711 obs1->add(_x.arg(),true) ;
712
713 obs1->add(_cacheObs) ;
714
715 } else {
716
717 // If cacheObs was filled, cache only observables in there
718 if (!_cacheObs.empty()) {
719 RooArgSet killList ;
720 for(auto const& arg : *obs1) {
721 if (arg->IsA()->InheritsFrom(RooAbsReal::Class()) && !_cacheObs.find(arg->GetName())) {
722 killList.add(*arg) ;
723 }
724 }
725 obs1->remove(killList) ;
726 }
727
728
729 // Make sure convolution observable is always in there
730 obs1->add(_x.arg(),true) ;
731
732 }
733
734 return RooFit::OwningPtr<RooArgSet>{obs1};
735}
736
737
738
739////////////////////////////////////////////////////////////////////////////////
740/// Return the parameters on which the cache depends given normalization
741/// set nset. For this p.d.f these are the parameters of the input p.d.f.
742/// but never the convolution variable, in case it is not part of nset.
743
745{
746 auto vars = getVariables() ;
747 vars->remove(*std::unique_ptr<RooArgSet>{actualObservables(nset)});
748
749 return RooFit::OwningPtr<RooArgSet>{std::move(vars)};
750}
751
752
753
754////////////////////////////////////////////////////////////////////////////////
755/// Return p.d.f. observable (which can be a function) to substitute given
756/// p.d.f. observable. Substitutes x by xprime if xprime is set.
757
759{
760 if (_xprime.absArg() && string(histObservable.GetName())==_x.absArg()->GetName()) {
761 return (*_xprime.absArg()) ;
762 }
763 return histObservable ;
764}
765
766
767
768////////////////////////////////////////////////////////////////////////////////
769/// Create appropriate generator context for this convolution. If both input p.d.f.s support
770/// internal generation, if it is safe to use them and if no observables other than the convolution
771/// observable are requested for generation, use the specialized convolution generator context
772/// which implements a smearing strategy in the convolution observable. If not return the
773/// regular accept/reject generator context
774
776 const RooArgSet* auxProto, bool verbose) const
777{
778 RooArgSet vars2(vars) ;
779 vars2.remove(_x.arg(),true,true) ;
780 Int_t numAddDep = vars2.getSize() ;
781
782 RooArgSet dummy ;
783 bool pdfCanDir = (((RooAbsPdf&)_pdf1.arg()).getGenerator(_x.arg(),dummy) != 0 && \
785 bool resCanDir = (((RooAbsPdf&)_pdf2.arg()).getGenerator(_x.arg(),dummy) !=0 &&
787
788 if (pdfCanDir) {
789 cxcoutI(Generation) << "RooFFTConvPdf::genContext() input p.d.f " << _pdf1.arg().GetName()
790 << " has internal generator that is safe to use in current context" << endl ;
791 }
792 if (resCanDir) {
793 cxcoutI(Generation) << "RooFFTConvPdf::genContext() input p.d.f. " << _pdf2.arg().GetName()
794 << " has internal generator that is safe to use in current context" << endl ;
795 }
796 if (numAddDep>0) {
797 cxcoutI(Generation) << "RooFFTConvPdf::genContext() generation requested for observables other than the convolution observable " << _x.arg().GetName() << endl ;
798 }
799
800
801 if (numAddDep>0 || !pdfCanDir || !resCanDir) {
802 // Any resolution model with more dependents than the convolution variable
803 // or pdf or resmodel do not support direct generation
804 cxcoutI(Generation) << "RooFFTConvPdf::genContext() selecting accept/reject generator context because one or both of the input "
805 << "p.d.f.s cannot use internal generator and/or "
806 << "observables other than the convolution variable are requested for generation" << endl ;
807 return new RooGenContext(*this,vars,prototype,auxProto,verbose) ;
808 }
809
810 // Any other resolution model: use specialized generator context
811 cxcoutI(Generation) << "RooFFTConvPdf::genContext() selecting specialized convolution generator context as both input "
812 << "p.d.fs are safe for internal generator and only "
813 << "the convolution observables is requested for generation" << endl ;
814 return new RooConvGenContext(*this,vars,prototype,auxProto,verbose) ;
815}
816
817
818
819////////////////////////////////////////////////////////////////////////////////
820/// Change the size of the buffer on either side of the observable range to `frac` times the
821/// size of the range of the convolution observable.
822
824{
825 if (frac<0) {
826 coutE(InputArguments) << "RooFFTConvPdf::setBufferFraction(" << GetName() << ") fraction should be greater than or equal to zero" << endl ;
827 return ;
828 }
829 _bufFrac = frac ;
830
831 // Sterilize the cache as certain partial results depend on buffer fraction
833}
834
835
836////////////////////////////////////////////////////////////////////////////////
837/// Change strategy to fill the overflow buffer on either side of the convolution observable range.
838///
839/// - `Extend` means is that the input p.d.f convolution observable range is widened to include the buffer range
840/// - `Flat` means that the buffer is filled with the p.d.f. value at the boundary of the observable range
841/// - `Mirror` means that the buffer is filled with a mirror image of the p.d.f. around the convolution observable boundary
842///
843/// The default strategy is extend. If one of the input p.d.f.s is a RooAddPdf, it is configured so that the interpretation
844/// range of the fraction coefficients is kept at the nominal convolutions observable range (instead of interpreting coefficients
845/// in the widened range including the buffer).
846
848{
849 _bufStrat = bs ;
850}
851
852
853
854////////////////////////////////////////////////////////////////////////////////
855/// Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the
856/// product operator construction
857
858void RooFFTConvPdf::printMetaArgs(ostream& os) const
859{
860 os << _pdf1.arg().GetName() << "(" << _x.arg().GetName() << ") (*) " << _pdf2.arg().GetName() << "(" << _x.arg().GetName() << ") " ;
861}
862
863
864
865////////////////////////////////////////////////////////////////////////////////
866/// (Re)calculate effective parameters of this p.d.f.
867
869{
870 RooArgSet argSet{_x.arg()};
871 RooArgSet params1;
872 RooArgSet params2;
873 _pdf1.arg().getParameters(&argSet, params1);
874 _pdf2.arg().getParameters(&argSet, params2);
876 _params.add(params1) ;
877 _params.add(params2,true) ;
878}
#define a(i)
Definition RSha256.hxx:99
#define coutI(a)
#define cxcoutI(a)
#define oocoutW(o, a)
#define coutF(a)
#define coutE(a)
int Int_t
Definition RtypesCore.h:45
#define ClassImp(name)
Definition Rtypes.h:377
#define N
char name[80]
Definition TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2467
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:74
void setOperMode(OperMode mode, bool recurseADirty=true)
Set the operation mode of this node.
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
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.
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
virtual RooAbsArg * cloneTree(const char *newname=nullptr) const
Clone tree expression of objects.
void attachDataSet(const RooAbsData &set)
Replace server nodes with names matching the dataset variable names with those data set variables,...
RooAbsBinning is the abstract base class for RooRealVar binning definitions.
int binNumber(double x) const
Returns the bin number corresponding to the value x.
Int_t numBins() const
Return number of bins.
virtual bool isUniform() const
virtual double highBound() const =0
virtual double lowBound() const =0
virtual RooAbsBinning * clone(const char *name=nullptr) const =0
RooAbsCachedPdf is the abstract base class for p.d.f.s that need or want to cache their evaluate() ou...
virtual const char * binningName() const
RooObjCacheManager _cacheMgr
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
void setDirtyProp(bool flag)
Control propagation of dirty flags from observables in dataset.
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
Abstract base class for objects that are lvalues, i.e.
virtual Int_t numBins(const char *rangeName=nullptr) const =0
virtual bool isDirectGenSafe(const RooAbsArg &arg) const
Check if given observable can be safely generated using the pdfs internal generator mechanism (if tha...
virtual Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Int_t numBins(const char *rangeName=nullptr) const override
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
void setBin(Int_t ibin, const char *rangeName=nullptr) override
Set value to center of bin 'ibin' of binning 'rangeName' (or of default binning if no range is specif...
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:91
static TClass * Class()
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * absArg() const
Return pointer to contained argument.
Definition RooArgProxy.h:47
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
void removeAll() override
Remove all argument inset using remove(const RooAbsArg&).
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, bool verbose=false)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
The RooDataHist is a container class to hold N-dimensional binned data.
Definition RooDataHist.h:39
void set(std::size_t binNumber, double weight, double wgtErr)
Set bin content of bin that was last loaded with get(std::size_t).
const RooArgSet * get() const override
Get bin centre of current bin.
Definition RooDataHist.h:76
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:57
std::unique_ptr< TVirtualFFT > fftr2c2
std::unique_ptr< TVirtualFFT > fftc2r
RooArgList containedArgs(Action) override
Returns all RooAbsArg objects contained in the cache element.
std::unique_ptr< RooAbsPdf > pdf2Clone
FFTCacheElem(const RooFFTConvPdf &self, const RooArgSet *nset)
Clone input pdf and attach to dataset.
std::unique_ptr< RooAbsBinning > histBinning
std::unique_ptr< RooAbsBinning > scanBinning
std::unique_ptr< RooAbsPdf > pdf1Clone
std::unique_ptr< TVirtualFFT > fftr2c1
PDF for the numerical (FFT) convolution of two PDFs.
friend class RooConvGenContext
RooSetProxy _params
Effective parameters of this p.d.f.
BufStrat _bufStrat
void calcParams()
(Re)calculate effective parameters of this p.d.f.
void setBufferFraction(double frac)
Change the size of the buffer on either side of the observable range to frac times the size of the ra...
double bufferFraction() const
Return value of buffer fraction applied in FFT calculation array beyond either end of the observable ...
TString histNameSuffix() const override
Suffix for cache histogram (added in addition to suffix for cache name)
void prepareFFTBinning(RooRealVar &convVar) const
Try to improve the binning and inform user if possible.
void fillCacheSlice(FFTCacheElem &cache, const RooArgSet &slicePosition) const
Fill a slice of cachePdf with the output of the FFT convolution calculation.
RooRealProxy _xprime
Input function representing value of convolution observable.
std::vector< double > scanPdf(RooRealVar &obs, RooAbsPdf &pdf, const RooDataHist &hist, const RooArgSet &slicePos, Int_t &N, Int_t &N2, Int_t &zeroBin, double shift) const
Scan the values of 'pdf' in observable 'obs' using the bin values stored in 'hist' at slice position ...
RooRealProxy _pdf1
First input p.d.f.
RooRealProxy _x
Convolution observable.
const char * inputBaseName() const override
Return base name component for cache components in this case 'PDF1_CONV_PDF2'.
RooAbsArg & pdfObservable(RooAbsArg &histObservable) const override
Return p.d.f.
~RooFFTConvPdf() override
Destructor.
RooFit::OwningPtr< RooArgSet > actualParameters(const RooArgSet &nset) const override
Return the parameters on which the cache depends given normalization set nset.
RooRealProxy _pdf2
Second input p.d.f.
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooNumConvPdf to more intuitively reflect the contents of the p...
void setBufferStrategy(BufStrat bs)
Change strategy to fill the overflow buffer on either side of the convolution observable range.
PdfCacheElem * createCache(const RooArgSet *nset) const override
Return specialized cache subclass for FFT calculations.
RooAbsGenContext * genContext(const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const override
Create appropriate generator context for this convolution.
RooFit::OwningPtr< RooArgSet > actualObservables(const RooArgSet &nset) const override
Return the observables to be cached given the normalization set nset.
void fillCacheObject(PdfCacheElem &cache) const override
Fill the contents of the cache the FFT convolution output.
RooSetProxy _cacheObs
Non-convolution observables that are also cached.
Class RooGenContext implement a universal generator context for all RooAbsPdf classes that do not hav...
RooLinearVar is the most general form of a derived real-valued object that can be used by RooRealInte...
void sterilize() override
Clear the cache payload but retain slot mapping w.r.t to normalization and integration sets.
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:40
bool hasBinning(const char *name) const override
Returns true if variable has a binning named 'name'.
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
void setBinning(const RooAbsBinning &binning, const char *name=nullptr)
Add given binning under name 'name' with this variable.
void setRange(const char *name, double min, double max)
Set a fit or plotting range.
const T & arg() const
Return reference to object held in proxy.
RooUniformBinning is an implementation of RooAbsBinning that provides a uniform binning in 'n' bins b...
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
TString fName
Definition TNamed.h:32
Basic string class.
Definition TString.h:139
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
RooConstVar & RooConst(double val)
const Int_t n
Definition legend1.C:16
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:43