Hi,
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.setBins(10000,"cache");
mass.setMin("cache",0);
mass.setMax("cache",6000); ## need to be adjusted to be higher than limit setting
sigpdf = RooFFTConvPdf("sigpdf","sigpdf",mass,bw,cb)
getattr(wspace,'import')(sigpdf,ROOT.RooCmdArg())
else:
wspace.factory(“Voigtian::sigpdf(mass_scaled, peak, width, sigma)”)
[/code]
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();
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
TFile *f = new TFile(“dimuon_13TeV_2016.root”);
RooWorkspace * _myWS = (RooWorkspace*) f->Get(“dimuon_13TeV_2016”);
_myWS->Print();
_myWS->var(“peak”)->setVal(1000);
ModelConfig* mc = (ModelConfig*) _myWS->obj(“mc”);
RooDataSet* prunedData = (RooDataSet*)_myWS->data(“data”);
_myWS->var(“ratio”)->setRange(0, 1.08e-5);
_myWS->var(“ratio”)->setVal(1e-5);
//getting estimate of limit from PLC
ProfileLikelihoodCalculator plc(*prunedData, *mc);
plc.SetConfidenceLevel(0.95);
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.SetConfidenceLevel(0.95);
mcmc.SetNumIters(10000); // Metropolis-Hastings algorithm iterations
mcmc.SetProposalFunction(sp);
mcmc.SetNumBurnInSteps(1000); // first N steps to be ignored as burn-in
mcmc.SetLeftSideTailFraction(0.0);
mcmc.SetNumBins(100);
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;
}[/code]
Attached the macro to run over the two workspaces we use.
Cheers,
Jan and Santiago
dimuon_13TeV_2016_CB.root (166 KB)
dimuon_13TeV_2016.root (68.8 KB)
cbLimit_Minimal.C (2.03 KB)