Dear experts,
I am trying to perform a simultaneous fit over a mass spectrum, where the mass variable has a different lower bound for each category. I have created a simple toy macro to demonstrate, which is attached, along with the plots it generates. I have also broken down the logic below:
I start by initialising the RooRealVar that will be my fitting variables, deltaMass, and defining the 2 different ranges:
int delta_low = 0;
int delta_high = 210;
int delta_nbins = (delta_high - delta_low) / 2;
RooRealVar deltaMass("Delta_M", "m[D^{*0}] - m[D^{0}]", delta_low, delta_high,
"MeV/c^{2}");
deltaMass.setRange("range_pi0", 135, 210);
deltaMass.setRange("range_gamma", 50, 210);
deltaMass.setBins(delta_nbins);
I have 2 neutral categories: gamma and pi0 (these represent the 2 neutral particles in the D*0 -> D0 pi0/gamma decay). I need the fit to the D*0 - D0 mass variable to be performed over the range 135-210 MeV for the pi0 category, and 50-210 MeV for the gamma category.
I then define the Categories and PDFs for the 2 neutrals:
// DEFINE INDEX CATEGORIES
RooCategory neutral("neutral", "neutral");
neutral.defineType("pi0");
neutral.defineType("gamma");
// DEFINE pi0 PDFs
// Signal
RooRealVar meanDeltaSigPi0("meanDeltaSigPi0", "Mean of Sig m[Delta]", 145,
140, 150);
RooRealVar sigmaDeltaSigPi0("sigmaDeltaSigPi0", "Sigma of pi Delta Gaussian",
3, 0, 5);
RooGaussian pdfDeltaSigPi0("pdfDeltaSigPi0", "Sig pi Delta PDF", deltaMass,
meanDeltaSigPi0, sigmaDeltaSigPi0);
// Background
RooRealVar thDeltaBkgPi0("thDeltaBkgPi0", "", 1.3415e+02); //, 30, 50);
RooRealVar aDeltaBkgPi0("aDeltaBkgPi0", "", 5.3768e-01); //, -1, 1);
RooRealVar bDeltaBkgPi0("bDeltaBkgPi0", "", -4.4782e-01); //, -1, 1);
RooRealVar cDeltaBkgPi0("cDeltaBkgPi0", "", 5.8018e+00); //, 0, 35);
RooDstD0BG pdfDeltaBkgPi0("pdfDeltaBkgPi0", "", deltaMass, thDeltaBkgPi0,
cDeltaBkgPi0, aDeltaBkgPi0, bDeltaBkgPi0);
// ROOADDPDF
RooRealVar yieldSigPi0("yieldSigPi0", "", 1000, 0, 5000);
RooRealVar yieldBkgPi0("yieldBkgPi0", "", 10000, 0, 50000);
RooArgSet yieldsPi0;
yieldsPi0.add(yieldSigPi0);
yieldsPi0.add(yieldBkgPi0);
RooArgSet functionsPi0;
functionsPi0.add(pdfDeltaSigPi0);
functionsPi0.add(pdfDeltaBkgPi0);
RooAddPdf pdfPi0("pdfPi0", "", functionsPi0, yieldsPi0);
// DEFINE gamma PDFs
// Signal Gamma
RooRealVar meanDeltaSigGamma1("meanDeltaSigGamma1", "Mean of Sig m[Delta]",
145, 140, 150);
RooRealVar sigmaDeltaSigGamma1("sigmaDeltaSigGamma1",
"Sigma of pi Delta Gaussian", 10, 0, 20);
RooGaussian pdfDeltaSigGamma1("pdfDeltaSigGamma1", "Sig pi Delta PDF",
deltaMass, meanDeltaSigGamma1,
sigmaDeltaSigGamma1);
// Signal Pi0 in Gamma mode
RooRealVar meanDeltaSigGamma2("meanDeltaSigGamma2", "Mean of Sig m[Delta]",
85, 80, 90);
RooRealVar sigmaDeltaSigGamma2("sigmaDeltaSigGamma2",
"Sigma of pi Delta Gaussian", 10, 0, 20);
RooGaussian pdfDeltaSigGamma2("pdfDeltaSigGamma2", "Sig pi Delta PDF",
deltaMass, meanDeltaSigGamma2,
sigmaDeltaSigGamma2);
// Background
RooRealVar thDeltaBkgGamma("thDeltaBkgGamma", "", 5.2160e+01); //, 30, 50);
RooRealVar aDeltaBkgGamma("aDeltaBkgGamma", "", 1.0341e+00); //, -1, 1);
RooRealVar bDeltaBkgGamma("bDeltaBkgGamma", "", -1.2284e+00); //, -1, 1);
RooRealVar cDeltaBkgGamma("cDeltaBkgGamma", "", 4.2450e+01); //, 0, 35);
RooDstD0BG pdfDeltaBkgGamma("pdfDeltaBkgGamma", "", deltaMass,
thDeltaBkgGamma, cDeltaBkgGamma, aDeltaBkgGamma,
bDeltaBkgGamma);
// ROOADDPDF
RooRealVar yieldSigGamma1("yieldSigGamma1", "", 2000, 0, 5000);
RooRealVar crossFeed("crossFeed", "", 1.5, 0, 5);
RooFormulaVar yieldSigGamma2("yieldSigGamma2", "", "@0*@1",
RooArgList(yieldSigPi0, crossFeed));
RooRealVar yieldBkgGamma("yieldBkgGamma", "", 20000, 0, 50000);
RooArgSet yieldsGamma;
yieldsGamma.add(yieldSigGamma1);
yieldsGamma.add(yieldSigGamma2);
yieldsGamma.add(yieldBkgGamma);
RooArgSet functionsGamma;
functionsGamma.add(pdfDeltaSigGamma1);
functionsGamma.add(pdfDeltaSigGamma2);
functionsGamma.add(pdfDeltaBkgGamma);
RooAddPdf pdfGamma("pdfGamma", "", functionsGamma, yieldsGamma);
I then generate the datasets for each category and cut them to be in the associated range:
// Generate datasets
auto dataPi0_1 =
std::unique_ptr<RooDataSet>(pdfPi0.generate(RooArgSet(deltaMass), 11000));
std::unique_ptr<RooDataSet> dataPi0;
dataPi0 =
std::unique_ptr<RooDataSet>(dynamic_cast<RooDataSet *>(dataPi0_1->reduce(
RooFit::CutRange("range==range::range_pi0"))));
if (dataPi0.get() == nullptr) {
throw std::runtime_error("Could not reduce pi0 input dataset.");
}
auto dataGamma_1 =
std::unique_ptr<RooDataSet>(pdfGamma.generate(RooArgSet(deltaMass), 11000));
std::unique_ptr<RooDataSet> dataGamma;
dataGamma =
std::unique_ptr<RooDataSet>(dynamic_cast<RooDataSet *>(dataGamma_1->reduce(
RooFit::CutRange("range==range::range_gamma"))));
if (dataGamma.get() == nullptr) {
throw std::runtime_error("Could not reduce gamma input dataset.");
}
Next, I make the combined dataset and simultaneous PDF:
RooDataSet combData("combData", "", RooArgSet(deltaMass),
RooFit::Index(neutral),
RooFit::Import("pi0", *dataPi0.get()),
RooFit::Import("gamma", *dataGamma.get()));
// CONSTRUCT SIMULTANEOUS PDF
RooSimultaneous simPdf("simPdf", "", neutral);
simPdf.addPdf(pdfPi0, "pi0");
simPdf.addPdf(pdfGamma, "gamma");
I then fit the simultaneous PDF to the combined dataset, with the Range, SumCoefRange and SplitRange functionality specified:
std::unique_ptr<RooFitResult> result = std::unique_ptr<RooFitResult>(
simPdf.fitTo(combData, RooFit::Extended(true),
RooFit::Range("range"),
RooFit::SumCoefRange("range"),
RooFit::SplitRange(true),
RooFit::Save(),
RooFit::Strategy(2), RooFit::Minimizer("Minuit2"),
RooFit::Offset(true)));
From the documentation, I understand that SplitRange should match a range to categories using rangeName_{INDEX}, where INDEX is the category label. I have therefore labelled my categories and ranges accordingly. However, this option doesn’t seem to be doing anything. The fit is still being performed over the entire range (0-210 MeV). I have also tried commenting out the Range and SumCoefRange options, but it makes no difference.
I can fit in one specific range by removing the SplitRange option and passing Range and SumCoefRange either “range_gamma” or “range_pi0”, but I cannot get the fit to work simultaneously over the 2 ranges.
I need the fit to be simultaneous because the signal shapes in the pi0 and gamma mode share some physics parameters (the magenta shapes in the plots). The reason I need the different ranges is is because the background PDFs, DstD0BG shapes, cause the fit to become unstable if you extend the fit range below the kinematic threshold.
Your help with this would be appreciated as I’ve searched the tutorial macros and can’t find any example of this functionality.
Best wishes,
Alexandra
SimToy.cpp (8.1 KB)
SimToy.pdf (35.7 KB)