134#ifndef ROOFIT_MATH_FFTW3
141 static void (*
doFFT)(
int,
double *,
double *,
double *) =
nullptr;
149 std::stringstream
ss;
150 ss <<
"RooFFTConvPdf evaluation Failed! The interpreter could not find the fftw3.h header.\n";
151 ss <<
"You have three options to fix this problem:\n";
152 ss <<
" 1) Install fftw3 on your system so that the interpreter can find it\n";
153 ss <<
" 2) In case fftw3.h is installed somewhere else,\n"
154 <<
" tell ROOT with gInterpreter->AddIncludePath() where to find it\n";
155 ss <<
" 3) Compile ROOT with the -Dfftw3=ON in the CMake configuration,\n"
156 <<
" such that ROOT comes with built-in fftw3 convolution routines\n";
157 oocoutE(
nullptr, Eval) <<
ss.str() << std::endl;
158 throw std::runtime_error(
"RooFFTConvPdf evaluation Failed! The interpreter could not find the fftw3.h header");
162void RooFFTConvPdf_doFFT(int n, double *input1, double *input2, double *output)
164 auto fftr2c1_Out = reinterpret_cast<fftw_complex *>(fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1)));
165 auto fftr2c2_Out = reinterpret_cast<fftw_complex *>(fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1)));
166 auto fftc2r_In = reinterpret_cast<fftw_complex *>(fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1)));
168 fftw_plan fftr2c1_plan = fftw_plan_dft_r2c(1, &n, input1, fftr2c1_Out, FFTW_ESTIMATE);
169 fftw_plan fftr2c2_plan = fftw_plan_dft_r2c(1, &n, input2, fftr2c2_Out, FFTW_ESTIMATE);
170 fftw_plan fftc2r_plan = fftw_plan_dft_c2r(1, &n, fftc2r_In, output, FFTW_ESTIMATE);
172 // Real->Complex FFT Transform on p.d.f. samplings
173 fftw_execute(fftr2c1_plan);
174 fftw_execute(fftr2c2_plan);
176 // Loop over first half +1 of complex output results, multiply
177 // and set as input of reverse transform
178 for (Int_t i = 0; i < n / 2 + 1; i++) {
179 double re1 = fftr2c1_Out[i][0];
180 double re2 = fftr2c2_Out[i][0];
181 double im1 = fftr2c1_Out[i][1];
182 double im2 = fftr2c2_Out[i][1];
183 fftc2r_In[i][0] = re1 * re2 - im1 * im2;
184 fftc2r_In[i][1] = re1 * im2 + re2 * im1;
187 // Reverse Complex->Real FFT transform product
188 fftw_execute(fftc2r_plan);
190 fftw_destroy_plan(fftr2c1_plan);
191 fftw_destroy_plan(fftr2c2_plan);
192 fftw_destroy_plan(fftc2r_plan);
194 fftw_free(fftr2c1_Out);
195 fftw_free(fftr2c2_Out);
196 fftw_free(fftc2r_In);
200 doFFT = reinterpret_cast<void(*)(
int,
double*,
double*,
double*)
>(
gInterpreter->ProcessLine(
"RooFFTConvPdf_doFFT;"));
208using std::endl, std::string, std::ostream;
226 _x(
"!x",
"Convolution Variable",
this, convVar),
227 _xprime(
"!xprime",
"External Convolution Variable",
this,
false),
228 _pdf1(
"!pdf1",
"pdf1",
this, pdf1,
false),
229 _pdf2(
"!pdf2",
"pdf2",
this, pdf2,
false),
230 _params(
"!params",
"effective parameters",
this),
234 _shift2((convVar.getMax(
"cache") + convVar.getMin(
"cache")) / 2),
252 _xprime(
"!xprime",
"External Convolution Variable",
this,
pdfConvVar),
253 _pdf1(
"!pdf1",
"pdf1",
this, pdf1,
false),
254 _pdf2(
"!pdf2",
"pdf2",
this, pdf2,
false),
255 _params(
"!params",
"effective parameters",
this),
259 _shift2((convVar.getMax(
"cache") + convVar.getMin(
"cache")) / 2),
279 _bufFrac(
other._bufFrac),
280 _bufStrat(
other._bufStrat),
281 _shift1(
other._shift1),
282 _shift2(
other._shift2),
283 _cacheObs(
"!cacheObs",
this,
other._cacheObs)
309 coutI(Caching) <<
"Changing internal binning of variable '" << convVar.
GetName()
310 <<
"' in FFT '" <<
fName <<
"'"
312 <<
" to " <<
optimal <<
" to improve the precision of the numerical FFT."
313 <<
" This can be done manually by setting an additional binning named 'cache'." << std::endl;
316 coutE(Caching) <<
"The internal binning of variable " << convVar.
GetName()
317 <<
" is not uniform. The numerical FFT will likely yield wrong results." << std::endl;
331 name.Append(
"_CONV_") ;
368 if (
self._shift1!=0) {
385 if (
self._shift2!=0) {
427 oocoutW(&
self, Eval) <<
"The FFT convolution '" <<
self.GetName() <<
"' will run with " <<
N
428 <<
" bins. A decent accuracy for difficult convolutions is typically only reached with n >= 1000. Suggest to increase the number"
429 " of bins of the observable '" <<
convObs->GetName() <<
"'." << std::endl;
462 ret.add(*pdf1Clone) ;
463 ret.add(*pdf2Clone) ;
464 if (pdf1Clone->ownedComponents()) {
465 ret.add(*pdf1Clone->ownedComponents()) ;
467 if (pdf2Clone->ownedComponents()) {
468 ret.add(*pdf2Clone->ownedComponents()) ;
512 std::vector<RooAbsLValue*>
obsLV(
n);
586#ifndef ROOFIT_MATH_FFTW3
589 std::vector<double>
output(N2);
604 if (
aux.fftr2c1 ==
nullptr ||
aux.fftr2c2 ==
nullptr ||
aux.fftc2r ==
nullptr) {
605 coutF(Eval) <<
"RooFFTConvPdf::fillCacheSlice(" <<
GetName() <<
"Cannot get a handle to fftw. Maybe ROOT was built without it?" << std::endl;
606 throw std::runtime_error(
"Cannot get a handle to fftw.");
612 aux.fftr2c1->Transform();
616 aux.fftr2c2->Transform();
620 for (
Int_t i=0 ; i<N2/2+1 ; i++) {
625 aux.fftr2c1->GetPointComplex(i,
re1,
im1) ;
626 aux.fftr2c2->GetPointComplex(i,
re2,
im2) ;
630 aux.fftc2r->SetPointComplex(i,t) ;
634 aux.fftc2r->Transform() ;
642 for (
Int_t i =0 ; i<
N ; i++) {
647 while (
j>=N2)
j-= N2 ;
650#ifndef ROOFIT_MATH_FFTW3
678 std::vector<double> array(N2);
685 if (
histX->getMax()>=0 &&
histX->getMin()<=0) {
687 }
else if (
histX->getMin()>0) {
702 std::vector<double> tmp(N2);
708 for (k=0 ; k<N2 ; k++) {
720 for (k=0 ; k<
Nbuf ; k++) {
723 for (k=0 ; k<
N ; k++) {
729 for (k=0 ; k<
Nbuf ; k++) {
730 tmp[
N+
Nbuf+k] = val ;
738 for (k=0 ; k<
N ; k++) {
742 for (k=1 ; k<=
Nbuf ; k++) {
752 for (
Int_t i=0 ; i<N2 ; i++) {
792 for(
auto const& arg : *
obs1) {
809 for(
auto const& arg : *
obs1) {
878 cxcoutI(Generation) <<
"RooFFTConvPdf::genContext() input p.d.f " <<
_pdf1.
arg().GetName()
879 <<
" has internal generator that is safe to use in current context" << std::endl ;
882 cxcoutI(Generation) <<
"RooFFTConvPdf::genContext() input p.d.f. " <<
_pdf2.
arg().GetName()
883 <<
" has internal generator that is safe to use in current context" << std::endl ;
886 cxcoutI(Generation) <<
"RooFFTConvPdf::genContext() generation requested for observables other than the convolution observable " <<
_x.
arg().GetName() << std::endl ;
893 cxcoutI(Generation) <<
"RooFFTConvPdf::genContext() selecting accept/reject generator context because one or both of the input "
894 <<
"p.d.f.s cannot use internal generator and/or "
895 <<
"observables other than the convolution variable are requested for generation" << std::endl ;
900 cxcoutI(Generation) <<
"RooFFTConvPdf::genContext() selecting specialized convolution generator context as both input "
901 <<
"p.d.fs are safe for internal generator and only "
902 <<
"the convolution observables is requested for generation" << std::endl ;
915 coutE(InputArguments) <<
"RooFFTConvPdf::setBufferFraction(" <<
GetName() <<
") fraction should be greater than or equal to zero" << std::endl ;
949 os <<
_pdf1.
arg().GetName() <<
"(" <<
_x.
arg().GetName() <<
") (*) " <<
_pdf2.
arg().GetName() <<
"(" <<
_x.
arg().GetName() <<
") " ;
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Common abstract base class for objects that represent a value and a "shape" in RooFit.
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
Abstract base class for RooRealVar binning definitions.
Abstract base class for p.d.f.s that need or want to cache their evaluate() output in a RooHistPdf de...
virtual const char * binningName() const
RooObjCacheManager _cacheMgr
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
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.
Abstract base class for generator contexts of RooAbsPdf objects.
Abstract base class for objects that are lvalues, i.e.
Abstract interface for all probability density functions.
Abstract base class for objects that represent a real value and implements functionality common to al...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooAbsArg * absArg() const
Return pointer to contained argument.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
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 ...
Container class to hold N-dimensional binned data.
const RooArgSet * get() const override
Get bin centre of current bin.
Container class to hold unbinned data.
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
PDF for the numerical (FFT) convolution of two PDFs.
friend class RooConvGenContext
RooSetProxy _params
Effective parameters of this p.d.f.
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.
friend class FFTCacheElem
RooSetProxy _cacheObs
Non-convolution observables that are also cached.
Implements a universal generator context for all RooAbsPdf classes that do not have or need a special...
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.
Variable that can be changed from the outside.
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.
const T & arg() const
Return reference to object held in proxy.
const char * GetName() const override
Returns name of object.
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)
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.