Convolution and sampling of multidimensional p.d.f.s

Dear experts,

I am trying to do convolution of a multidimensional p.d.f. as described in the RooFit manual:
M(x,y) = T(x,y) (X) R(x)

I am successful in convolving a 2D p.d.f. with a 1D p.d.f. However, my case calls for convolution of a 3D p.d.f. with 1D p.d.f., like:
M(x,y,z) = T(x,y,z) (X) R(x)

When I do this convolution I don’t have any problems until I try to generate some events. Please see the code snippet at the bottom of the message. My error occurs at the last line and at this point I get an error which says:

[#0] ERROR:InputArguments -- RooDataHist::weight(myPdfConv_myPDF_CONV_gauss_CACHEHIST_Obs[Y,hs,m]_BufFrac0.1_BufStrat0) interpolation in 3 dimensions not yet implemented

My question is if this is the expected behavior. Is that something that has not been implemented yet? Meaning, I can do 2D (X) 1D, but not 3D (X) 1D. Or have I simply misconfigured something to receive this error.

Thanks!

...
...
        // S e t u p   c o m p o n e n t   p d f s 
        // ---------------------------------------
        
        // Construct observable
        RooRealVar Y("Y","Y",-2.5,2.5) ;
        RooRealVar m("m","m",71.,111.) ;
        RooRealVar hs("hs","hs",-1.,1.) ;
        
        // Construct custom PDF
        RooRealVar par1("par1","par1",0.25,0.,1.) ;
        RooRealVar par2("par2","par2",2.4357,0.1,10) ;
        RooRealVar par3("par3","par3",2.1,0.1,10.) ;
        RooZxsec myPdf("myPDF","myPDF", Y, m, hs, par1, par2, par3);
        
        par2.setConstant( kTRUE );
        par3.setConstant( kTRUE );
        
        // Construct gauss(t,mg,sg) 
        RooRealVar mg("mg","mg",0) ;
        RooRealVar sg("sg","sg",2,0.1,10) ;
        RooGaussian gauss("gauss","gauss",m,mg,sg) ;

        mg.setConstant( kTRUE );
        sg.setConstant( kTRUE );

        // C o n s t r u c t   c o n v o l u t i o n   p d f 
        // ---------------------------------------
        
        // Set #bins to be used for FFT sampling to 10000
        m.setBins(1000,"cache") ; 
        
        // Construct landau (x) gauss
        RooFFTConvPdf myPdfConv("myPdfConv","Zxsec (X) gauss", m, myPdf, gauss) ;
        
        myPdfConv.setCacheObservables( RooArgSet(Y,hs) );
        
        // S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f 
        // ----------------------------------------------------------------------
        
        // Sample 1000 events in x from gxlx
        RooDataSet* data = myPdfConv.generate(RooArgSet(Y,m,hs),1000) ;
...
...

Hi,

A RooFFTConvPdf stores its full output (in all observables) internally in a RooHistPdf as the FFT convolution is calculated entire at once, in contrast to a regular pdf, where is each point in the pdf is evaluated at a time.

The default configuration of RooHistPdf that store this output is to have 2nd order interpolation enabled to promote smoothness of the output pdf, however RooHistPdf does not implement interpolation algortithms for >2 dimensions. This is not a (mathematically) fundamental limitation, but I have no generic N-dimensional interpolation algorithm at hand at the moment.

In the absense of that, there are two possibilities to work around this

a) request that no interpolation is performed in the caching RooHistPdf. This can be done in the ctor

Hi,

A RooFFTConvPdf stores its full output (in all observables) internally in a RooHistPdf as the FFT convolution is calculated entire at once, in contrast to a regular pdf, where is each point in the pdf is evaluated at a time.

The default configuration of RooHistPdf that store this output is to have 2nd order interpolation enabled to promote smoothness of the output pdf, however RooHistPdf does not implement interpolation algortithms for >2 dimensions. This is not a (mathematically) fundamental limitation, but I have no generic N-dimensional interpolation algorithm at hand at the moment.

In the absense of that, there are two possibilities to work around this

a) request that no interpolation is performed in the caching RooHistPdf. This can be done in the ctor:

RooFFTConvPdf(const char *name, const char *title, RooRealVar& convVar, RooAbsPdf& pdf1, RooAbsPdf& pdf2, Int_t ipOrder=2);

i.e. set the last parameter to zero.

b) request that the result be cached in a lower dimensional histogram. To do so call

 void RooFFTConfPdf::setCacheObservables(const RooArgSet& obs) 

 with the variables that you'd like to see cached. One option is to only cache the convolution
 observable. This will however trigger recalculation of the cache whenever one of your
 other (y,z) observables changes value, so it is potentially very time consuming, but not necessarily:
 it depends on how many events you have in your likelihood. If it is a small number (compared to
 the number of bins in the (y,z) space for caching) this might still work out OK.

Wouter

1 Like

Hi Wouter,

Thanks for the advice. I chose the solution to turn off interpolation and it seems to work as expected.

As I have gone on, I have run across another issue. I would like to do my convolution in one range of the convoluted observable, and then do a fit in another range. Using the example from my original posting, I did something like this:

...
...
	// C o n s t r u c t   c o n v o l u t i o n   p d f 
	// ---------------------------------------
	
	// Construct landau (x) gauss
	RooFFTConvPdf myPdfConv("myPdfConv","Zxsec (X) gauss", m, myPdf, myCBGaus, 0) ;

	myPdfConv.setCacheObservables( RooArgSet(m,Y,hs) );
	myPdfConv.setBufferFraction(0.7);
	
	// Set #bins to be used for FFT sampling to 10000
	m.setBins(1000,"cache") ; 
	Y.setBins(1000,"cache") ; 
	hs.setBins(1000,"cache") ; 
	
	// S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f 
	// ----------------------------------------------------------------------
	
	// Read in data from text file
	RooDataSet* data = RooDataSet::read("data/dat.txt", RooArgList(Y,m,hs));

	// define fit range
	std::cout << "define fit range..." << std::endl;
	m.setRange("fitm", 81., 101.);
	Y.setRange("fitm", -2.1,2.1);
	hs.setRange("fitm", -1.,1.);
	
	std::cout << "beginning fit..." << std::endl;
	myPdfConv.fitTo(*data, Range("fitm") ) ;

When I run this, I get some strange results. If I set it so, in the code above, the two ranges are the same, I get a different fit result than when I do not set it explicity (define no fit range). Also, the errors of the fit are always about the same, regardless of the range, which is contrary to expectations. The output from the fitter is:

define fit range...
[#1] INFO:Eval -- RooRealVar::setRange(m) new range named 'fitm' created with bounds [86,96]
[#1] INFO:Eval -- RooRealVar::setRange(Y) new range named 'fitm' created with bounds [-2.1,2.1]
[#1] INFO:Eval -- RooRealVar::setRange(hs) new range named 'fitm' created with bounds [-1,1]
beginning fit...
[#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_myPdfConv_dataset) constructing test statistic for sub-range named fitm
[#1] INFO:Eval -- RooRealVar::setRange(m) new range named 'NormalizationRangeForfitm' created with bounds [66,116]
[#1] INFO:Eval -- RooRealVar::setRange(Y) new range named 'NormalizationRangeForfitm' created with bounds [-2.1,2.1]
[#1] INFO:Eval -- RooRealVar::setRange(hs) new range named 'NormalizationRangeForfitm' created with bounds [-1,1]
[#1] INFO:Eval -- RooRealVar::setRange(m) new range named 'fit_nll_myPdfConv_dataset' created with bounds [86,96]
[#1] INFO:Eval -- RooRealVar::setRange(Y) new range named 'fit_nll_myPdfConv_dataset' created with bounds [-2.1,2.1]
[#1] INFO:Eval -- RooRealVar::setRange(hs) new range named 'fit_nll_myPdfConv_dataset' created with bounds [-1,1]
[#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_myPdfConv_dataset) fixing interpretation of coefficients of any RooAddPdf to full domain of observables 
[#1] INFO:Eval -- RooRealVar::setRange(m) new range named 'refrange_fft_myPdfConv' created with bounds [66,116]
[#1] INFO:NumericIntegration -- RooRealIntegral::init(myPDF_Int[Y,hs,m]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(Y,m)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(myPdfConv) creating new cache 0x2daf8f0 with pdf myPDF_CONV_myCBGaus_CACHE_Obs[Y,m,hs] for nset (Y,m,hs,blindState) with code 0

Do you know why this is happening, have I configured it in the right way? Is there something I have over-looked? Or else it should be straightforward, right?

Thanks again,
Nhan