Problem using convolution with RooFFTConvPDF


we are performing an unbinned shape analysis where the signal is modelled by a Breit-Wigner convolved with a resolution model. Using a Gaussian for this as implemnted in RooVoigtian works fine. However, we are trying to use a crystal ball for the resolution and encounter problems. For this we are using RooFFTConvPDF to convolve Breit-Wigner and Crystal ball. We try to set limits on a search for a narrow resonance on the invariant mass distribution. The code runs, but we are seeing limits that are a factor three stronger with respect to the Voigtian, while we would expect slightly weaker limits because of the enlarged width. This behaviour can be seen using both the MCMCCalculator and the ProfileLikelihoodCalculator.

The signal models are build this way:

[code] if doCB:
print " – build convolution"
wspace.factory(“BreitWigner::bw(mass_scaled, peak, width)”)
wspace.factory(“RooCBShape::cb(mass_scaled, mean[0.0], sigma, alpha[1.43], n[3])”)

   bw = wspace.pdf("bw")
   cb = wspace.pdf("cb")

   mass.setMax("cache",6000);  ## need to be adjusted to be higher than limit setting

   sigpdf = RooFFTConvPdf("sigpdf","sigpdf",mass,bw,cb)

wspace.factory(“Voigtian::sigpdf(mass_scaled, peak, width, sigma)”)

As nothing else changes between the two approaches, we suspect the convolution to create problems here. Any ideas what is going wrong?

A minimal example of the limit setting code looks like this:

[code]#include “RooStats/SimpleInterval.h”
#include “RooStats/BayesianCalculator.h”
#include “RooStats/MCMCCalculator.h”
#include “RooStats/MCMCIntervalPlot.h”
#include “RooStats/SequentialProposal.h”
#include “RooStats/ProposalHelper.h”
#include “RooStats/ProfileLikelihoodCalculator.h”
#include “RooStats/LikelihoodInterval.h”
#include “RooStats/LikelihoodIntervalPlot.h”
#include “RooStats/HypoTestResult.h”

using namespace std;
using namespace RooFit;
using namespace RooStats;

void cbLimit_Minimal(){

RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();

TFile *f = new TFile(“dimuon_13TeV_2016.root”);

RooWorkspace * _myWS = (RooWorkspace*) f->Get(“dimuon_13TeV_2016”);



ModelConfig* mc = (ModelConfig*) _myWS->obj(“mc”);

RooDataSet* prunedData = (RooDataSet*)_myWS->data(“data”);

_myWS->var(“ratio”)->setRange(0, 1.08e-5);

//getting estimate of limit from PLC
ProfileLikelihoodCalculator plc(*prunedData, *mc);

mpPlrInt = plc.GetInterval();
RooRealVar * _poi = (RooRealVar *)mc->GetParametersOfInterest()->first();
double upper_limit = mpPlrInt->UpperLimit( *_poi );

cout << "Upper limit from PCL: " << upper_limit << endl;

cout << “Now testing MCMC” << endl;

MCMCInterval * mcInt;

// FIXME: testing: this proposal function seems fairly robust
SequentialProposal sp(0.5);

MCMCCalculator mcmc( *prunedData, *mc );
mcmc.SetNumIters(10000); // Metropolis-Hastings algorithm iterations


mcmc.SetNumBurnInSteps(1000); // first N steps to be ignored as burn-in

mcInt = mcmc.GetInterval();

RooRealVar * p_first_poi = (RooRealVar*) mc->GetParametersOfInterest()->first();
double poi_limit = mcInt->UpperLimit(*p_first_poi);

cout << "limti from MCMC: " << poi_limit << endl;


Attached the macro to run over the two workspaces we use.

Jan and Santiago
dimuon_13TeV_2016_CB.root (166 KB)
dimuon_13TeV_2016.root (68.8 KB)
cbLimit_Minimal.C (2.03 KB)


After investigating carefully, I have not spotted any particular problem with the FFT model. There could be some numerical issues who can cause this effect.
I have noticed that for example the fit fails to get the right value when injecting a fake signal.

I would try to re-parametrize the parameter of interest “ratio”, to be in a more sensible range (e.g. [0,10] ) by rescaling it by a 10^-5 or 6 factor.