Dear experts,
I am trying to perform template morphing using RooMomentMorph with RooHistPdf. The morphed function seems to produce inconsistent binning when drawn, which may cause issues when fitting this function as part of a composite signal + background model.
Here is the setup:
-
I create a morphed function using shifted-mean Gaussians. When I draw this morphed PDF, the binning does not match the original histogram binning.
-
I then construct background models using PiecewiseInterpolation for systematic evaluation, build the full signal+background model, generate a toy dataset, and perform a fit.
-
The fit only works with
EvalBackend("Legacy"), and drawing the fit results shows that the binning inconsistency from the morphing distorts the visualization and (maybe) could impact the fit.
-
How can I ensure that RooMomentMorph respects the original histogram binning when drawn?
-
Is this behavior expected, and are there recommended workarounds to ensure consistent binning for fitting and plotting in a composite model?
Here attached is the full running code I am using.
Thanks in advance for any guidance!
Root Version: 6.32.02
#include <RooRealVar.h>
#include <RooGaussian.h>
#include <RooDataHist.h>
#include <RooMomentMorph.h>
#include <RooFitResult.h>
#include <RooPlot.h>
#include <TCanvas.h>
#include <vector>
#include <iostream>
using namespace std;
using namespace RooFit;
RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
void analytical_morph()
{
// Observable
RooRealVar x("x", "Observable", -10, 10);
int nbins = 50;
x.setBins(nbins);
// Continuous mean parameter you want to interpolate over
RooRealVar mean("mean", "Mean parameter", 0.0, -4, 4);
RooRealVar alpha("alpha", "Alpha parameter", 0.0, -1.0, 1.0);
// Fixed sigma for all templates
RooConstVar sigma("sigma", "Width of Gaussians", 1);
// Discrete means for template histograms
std::vector<double> meanPoints = {-4, -3, -2, -1, 0, 1, 2, 3, 4};
// Storage for objects
std::vector<RooDataHist *> dataHists;
std::vector<RooHistPdf *> histPdfs;
RooArgList pdfList;
TVectorD mean_pts(meanPoints.size());
// Generate templates with consistent binning
int nEvents = 100000; // Number of events for template generation
for (size_t i = 0; i < meanPoints.size(); ++i)
{
mean_pts[i] = meanPoints[i];
// Create Gaussian PDF
RooConstVar *mu = new RooConstVar(Form("mu_%zu", i), "Mean", meanPoints[i]);
RooGaussian *gauss = new RooGaussian(Form("gauss_%zu", i), "Gaussian", x, *mu, sigma);
// Generate dataset and create histogram with consistent binning
RooDataSet *dataset = gauss->generate(x, nEvents);
TH1 *hist = dataset->createHistogram(Form("hist_%zu", i), x, Binning(nbins));
hist->Scale(1.0 / hist->Integral()); // Normalize
RooDataHist *dh = new RooDataHist(Form("dh_%zu", i), "DataHist", x, hist);
RooHistPdf *histPdf = new RooHistPdf(Form("hist_%zu", i), "HistPdf", x, *dh);
dataHists.push_back(dh);
histPdfs.push_back(histPdf);
pdfList.add(*histPdf);
// Clean up temporary objects
delete dataset;
delete hist;
}
// Create morphing PDF using template histograms
RooMomentMorph morph("morph", "Morphed PDF",
mean,
RooArgList(x),
pdfList,
mean_pts,
RooMomentMorph::Linear);
// Draw morph for several mean values
TCanvas *c1 = new TCanvas("c1", "Morphing Test", 800, 600);
RooPlot *frame = x.frame(Bins(nbins));
double step = 0.7;
std::vector<double> test_means;
std::vector<Color_t> colors;
for (double val = -4.0; val <= 4.0 + 1e-6; val += step)
{
test_means.push_back(val);
colors.push_back(static_cast<EColor>(static_cast<double>(kRed) +
(val + 4.0) / 8.0 * (static_cast<double>(kBlue) - static_cast<double>(kRed))));
}
for (size_t i = 0; i < test_means.size(); ++i)
{
mean.setVal(test_means[i]);
morph.plotOn(frame, LineColor(colors[i]), LineWidth(2),
Name(Form("morph_mean_%.1f", test_means[i])));
}
frame->SetTitle("Morph for different mean values");
frame->Draw();
// Create background PDFs with consistent binning
// Gaussian background shapes
RooGaussian bg_pdf("bg_pdf", "Background PDF", x, RooConst(0), RooConst(5));
RooGaussian bg_pdf_m("bg_pdf_m", "Background PDF -1", x, RooConst(-1), RooConst(5));
RooGaussian bg_pdf_p("bg_pdf_p", "Background PDF +1", x, RooConst(1), RooConst(5));
// Generate datasets and create histograms with consistent binning
RooDataSet *dset_bg = bg_pdf.generate(x, nEvents);
RooDataSet *dset_bg_m = bg_pdf_m.generate(x, nEvents);
RooDataSet *dset_bg_p = bg_pdf_p.generate(x, nEvents);
TH1 *h_bg = dset_bg->createHistogram("h_bg", x, Binning(nbins));
TH1 *h_bg_m = dset_bg_m->createHistogram("h_bg_m", x, Binning(nbins));
TH1 *h_bg_p = dset_bg_p->createHistogram("h_bg_p", x, Binning(nbins));
h_bg->Scale(1.0 / h_bg->Integral());
h_bg_m->Scale(1.0 / h_bg_m->Integral());
h_bg_p->Scale(1.0 / h_bg_p->Integral());
// Create RooDataHist objects
RooDataHist bg_dh("bg_dh", "Background DataHist", x, h_bg);
RooDataHist bg_dh_m("bg_dh_m", "Background DataHist -1", x, h_bg_m);
RooDataHist bg_dh_p("bg_dh_p", "Background DataHist +1", x, h_bg_p);
// Create HistPdf objects
RooHistPdf bg_histPdf("bg_histPdf", "Background HistPdf", x, bg_dh);
RooHistPdf bg_histPdf_m("bg_histPdf_m", "Background HistPdf -1", x, bg_dh_m);
RooHistPdf bg_histPdf_p("bg_histPdf_p", "Background HistPdf +1", x, bg_dh_p);
// Create interpolated background
PiecewiseInterpolation bg_pdf_temp("bg_pdf_temp", "Background PDF alpha",
bg_histPdf, bg_histPdf_m, bg_histPdf_p, alpha);
RooWrapperPdf bg_pdf_alpha("bg_pdf_alpha", "Background PDF alpha wrapper",
bg_pdf_temp, true);
// Exponential background shapes (same approach)
RooExponential expo_pdf("expo_pdf", "Exponential PDF", x, RooConst(-0.2));
RooExponential expo_pdf_m("expo_pdf_m", "Exponential PDF -0.1", x, RooConst(-0.3));
RooExponential expo_pdf_p("expo_pdf_p", "Exponential PDF +0.1", x, RooConst(-0.1));
RooDataSet *dset_expo = expo_pdf.generate(x, nEvents);
RooDataSet *dset_expo_m = expo_pdf_m.generate(x, nEvents);
RooDataSet *dset_expo_p = expo_pdf_p.generate(x, nEvents);
TH1 *h_expo = dset_expo->createHistogram("h_expo", x, Binning(nbins));
TH1 *h_expo_m = dset_expo_m->createHistogram("h_expo_m", x, Binning(nbins));
TH1 *h_expo_p = dset_expo_p->createHistogram("h_expo_p", x, Binning(nbins));
h_expo->Scale(1.0 / h_expo->Integral());
h_expo_m->Scale(1.0 / h_expo_m->Integral());
h_expo_p->Scale(1.0 / h_expo_p->Integral());
RooDataHist expo_dh("expo_dh", "Expo DataHist", x, h_expo);
RooDataHist expo_dh_m("expo_dh_m", "Expo DataHist -0.1", x, h_expo_m);
RooDataHist expo_dh_p("expo_dh_p", "Expo DataHist +0.1", x, h_expo_p);
RooHistPdf expo_histPdf("expo_histPdf", "Expo HistPdf", x, expo_dh);
RooHistPdf expo_histPdf_m("expo_histPdf_m", "Expo HistPdf -0.1", x, expo_dh_m);
RooHistPdf expo_histPdf_p("expo_histPdf_p", "Expo HistPdf +0.1", x, expo_dh_p);
PiecewiseInterpolation expo_pdf_temp("expo_pdf_temp", "Expo PDF alpha",
expo_histPdf, expo_histPdf_m, expo_histPdf_p, alpha);
RooWrapperPdf expo_pdf_alpha("expo_pdf_alpha", "Expo PDF alpha wrapper",
expo_pdf_temp, true);
// Full PDF
RooRealVar nsig("nsig", "Signal yield", 100, 0, 1000);
RooRealVar nbkg("nbkg", "Background yield", 1000, 0, 4000);
RooRealVar nbkg2("nbkg2", "Background yield 2", 500, 0, 2000);
RooAddPdf full_pdf("full_pdf", "Full PDF",
RooArgList(morph, bg_pdf_alpha, expo_pdf_alpha),
RooArgList(nsig, nbkg, nbkg2));
// Generate toy data
double mean_val = 2.3;
double alpha_val = 0.4;
alpha.setVal(alpha_val);
mean.setVal(mean_val);
RooDataHist *toyData = full_pdf.generateBinned(x, Extended());
// Reset for fit
mean.setVal(0.0);
alpha.setVal(0.0);
nsig.setVal(0);
nbkg.setVal(1600);
nbkg2.setVal(800);
std::cout << "Starting analytical fit..." << std::endl;
RooFitResult *fitResult = full_pdf.fitTo(*toyData,
Save(true),
Extended(),
Strategy(1),
PrintLevel(1),
EvalBackend("Legacy"));
fitResult->Print();
std::cout << "Fitted mean: " << mean.getVal() << " +/- " << mean.getError() << std::endl;
std::cout << "True mean: " << mean_val << std::endl;
std::cout << "Fitted alpha: " << alpha.getVal() << " +/- " << alpha.getError() << std::endl;
std::cout << "True alpha: " << alpha_val << std::endl;
// Plot fit result
TCanvas *c2 = new TCanvas("c2", "Analytical Morph Fit", 800, 600);
RooPlot *frame2 = x.frame(Bins(nbins));
toyData->plotOn(frame2);
full_pdf.plotOn(frame2, LineColor(kRed), LineWidth(3));
full_pdf.plotOn(frame2, Components("morph"), LineColor(kGreen), LineStyle(kDashed));
full_pdf.plotOn(frame2, Components("bg_pdf_alpha"), LineColor(kBlue), LineStyle(kDashed));
full_pdf.plotOn(frame2, Components("expo_pdf_alpha"), LineColor(kMagenta), LineStyle(kDashed));
frame2->SetTitle(Form("Analytical Fit - mean: %.3f #pm %.3f (true: %.3f)",
mean.getVal(), mean.getError(), mean_val));
frame2->Draw();
/*
// Cleanup
for (auto* ptr : dataHists) delete ptr;
for (auto* ptr : histPdfs) delete ptr;
delete dset_bg; delete dset_bg_m; delete dset_bg_p;
delete dset_expo; delete dset_expo_m; delete dset_expo_p;
delete h_bg; delete h_bg_m; delete h_bg_p;
delete h_expo; delete h_expo_m; delete h_expo_p;
delete toyData;
delete fitResult;
*/
}