Why the FitGauss(true) doesn't work for plotParam in RooMCStudy

Dear Experts,

I am using RooMCStudy to generate toyMC for simultaneous fit, as below;

 RooMCStudy *mcstudy = new RooMCStudy(simPdf, RooArgList(Mbc,deltaE, sample), Silence(), Extended(), FitOptions(Save(true), PrintEvalErrors(0)));
 mcstudy->generateAndFit(1001);
 RooPlot *frame1 = mcstudy->plotParam(Nsig, Bins(40), FitGauss(true));
 RooPlot *frame2 = mcstudy->plotError(Nsig, Bins(40));
 RooPlot *frame3 = mcstudy->plotPull(Nsig, Bins(40), FitGauss(true));

I plotted the Nsig, its error, and its pull. The Gaussian fit works well for pull distribution. But, it doesn’t perform fit for Nsig parameter. Could you please let me know how should I solve this issue?

Thanks,

Hello @CS_p ,

could you please share a minimal, self-contained reproducer for the issue that we can run in order to debug the problem?

@jonas might be able to help further.

Cheers,
Enrico

Hi @eguiraud,

Thank you for your response. Here is a working example;

using namespace RooFit;
using namespace RooStats;
void toyMC(){


 const int nbkg_mumu = 553;
 const int nbkg_ee = 389;
 const int nsig =0;
 RooRealVar Mbc("Mbc","Mbc",5.2,5.29);
 RooRealVar deltaE("deltaE","deltaE", -0.15, 0.1);

 //Mbc backgrund fit for mumu channel
 RooRealVar a0("a0","a0",0);
 RooArgusBG bkg_mbc_mumu("bkg_mbc_mumu","bkg_mbc_mumu",Mbc, RooConst(5.289), a0);
//Mbc signal fit for mumu channel
RooRealVar mean_mbc_mumu("mean_mbc_mumu","mean_mbc_mumu",5.27945);
RooRealVar sigma_mbc_mumu("sigma_mbc_mumu","sigma_mbc_mumu",0.00268);
RooRealVar a_mbc_mumu("a_mbc_mumu", "a_mbc_mumu", 2.02);
RooRealVar n_mbc_mumu("n_mbc_mumu", "n_mbc_mumu", 10.0);
RooCBShape sig_mbc_mumu("sig_mbc_mumu", "sig_mbc_mumu", Mbc, mean_mbc_mumu, sigma_mbc_mumu, a_mbc_mumu, n_mbc_mumu);
//Delta E background fit for mumu channel
RooRealVar pol1("pol1","pol1",-0.609);
RooChebychev bkg_de_mumu("bkg_de_mumu","bkg_de_mumu", deltaE,RooArgList(pol1));
//Delta E signal fit for mumu channel
RooRealVar mean_de_mumu("mean_de_mumu","mean_de_mumu",0.012);
RooRealVar lmbda_de_mumu("lmbda_de_mumu","lmbda_de_mumu",0.0311);
RooRealVar gamma_de_mumu("gamma_de_mumu", "gamma_de_mumu", 0.96);
RooRealVar delta_de_mumu ("delta_de_mumu", "delta_de_mumu", 1.41);
RooJohnson johnson_mumu("johnson_mumu", "johnson_mumu", deltaE, mean_de_mumu, lmbda_de_mumu, gamma_de_mumu, delta_de_mumu);
RooRealVar sigma_de_mumu ("sigma_de_mumu","sigma_de_mumu",0.0212);
RooGaussian gauss_de_mumu("gauss_de_mumu","gauss_de_mumu",deltaE,mean_de_mumu,sigma_de_mumu);
RooRealVar frac_de_mumu("frac_de_mumu","frac_de_mumu",0.75);
RooAddPdf sig_de_mumu("sig_de_mumu","sig_de_mumu",      RooArgList(johnson_mumu,gauss_de_mumu), RooArgList(frac_de_mumu));

//Mbc backgrund fit for ee channel
RooArgusBG bkg_mbc_ee ("bkg_mbc_ee","bkg_mbc_ee",Mbc,RooConst(5.289), a0);
//Mbc signal fit
RooRealVar mean_mbc_ee("mean_mbc_ee","mean_mbc_ee",5.27949);
RooRealVar  sigma_mbc_ee("sigma_mbc_ee","sigma_mbc_ee",0.00271);
RooRealVar a_mbc_ee("a_mbc_ee", "a_mbc_ee",1.71);
RooRealVar n_mbc_ee("n_mbc_ee", "n_mbc_ee", 10.0);
RooCBShape sig_mbc_ee("sig_mbc_ee", "sig_mbc_ee", Mbc, mean_mbc_ee, sigma_mbc_ee, a_mbc_ee, n_mbc_ee);
//Delta E background fit for ee channel
RooChebychev bkg_de_ee("bkg_de_ee","bkg_de_ee", deltaE, RooArgList(pol1));
//Delta E signal fit for ee channel
RooRealVar mean_de_ee("mean_de_ee","mean_de_ee",0.0171);
RooRealVar lmbda_de_ee("lmbda_de_ee","lmbda_de_ee",0.0336);
RooRealVar gamma_de_ee("gamma_de_ee", "gamma_de_ee",1.09);
RooRealVar delta_de_ee("delta_de_ee", "delta_de_ee",1.29);
RooJohnson sig_de_ee("sig_de_ee", "sig_de_ee", deltaE, mean_de_ee, lmbda_de_ee, gamma_de_ee, delta_de_ee);

//model
RooProdPdf sig_mumu("sig_mumu","sig_mumu", RooArgList(sig_mbc_mumu,sig_de_mumu));
RooProdPdf bkg_mumu("bkg_mumu","bkg_mumu",RooArgList(bkg_mbc_mumu,bkg_de_mumu));

RooProdPdf sig_ee("sig_ee","sig_ee", RooArgList(sig_mbc_ee,sig_de_ee));
RooProdPdf bkg_ee("bkg_ee","bkg_ee",RooArgList(bkg_mbc_ee,bkg_de_ee));

RooRealVar Nsig("Nsig","Nsig",nsig,-60,50);
RooRealVar Nbkg_mumu("Nbkg_mumu","Nbkg_mumu",nbkg_mumu);
RooRealVar Nbkg_ee("Nbkg_ee","Nbkg_ee",nbkg_ee);
RooAddPdf model_mumu("model_mumu","model_mumu",RooArgList(sig_mumu,bkg_mumu),RooArgList(Nsig,Nbkg_mumu));
RooAddPdf model_ee("model_ee","model_ee",RooArgList(sig_ee,bkg_ee),RooArgList(Nsig,Nbkg_ee));
  
RooCategory sample("sample", "sample");
sample.defineType("mumu");
sample.defineType("ee");

RooSimultaneous simPdf("simPdf", "simultaneous pdf", sample);
simPdf.addPdf(model_mumu,"mumu");
simPdf.addPdf(model_ee,"ee");

 RooMCStudy *mcstudy = new RooMCStudy(simPdf, RooArgList(Mbc,deltaE, sample), Silence(), Extended(), FitOptions(Save(true), PrintEvalErrors(0)));
 mcstudy->generateAndFit(1001);

 // E x p l o r e   r e s u l t s   o f   s t u d y
 // ------------------------------------------------

 // Make plots of the distributions of mean, the error on mean and the pull of mean
 RooPlot *frame1 = mcstudy->plotParam(Nsig, Bins(40), FitGauss(true));
 RooPlot *frame2 = mcstudy->plotError(Nsig, Bins(40));
 RooPlot *frame3 = mcstudy->plotPull(Nsig, Bins(40), FitGauss(true));

 // Plot distribution of minimized likelihood
 RooPlot *frame4 = mcstudy->plotNLL(Bins(40));

 // Draw all plots on a canvas
 gStyle->SetOptStat(0);
 TCanvas *c = new TCanvas(" ", " ", 900, 900);
 c->Divide(2, 2);
 c->cd(1);
 gPad->SetLeftMargin(0.15);
 frame1->GetYaxis()->SetTitleOffset(1.4);
 frame1->Draw();
 c->cd(2);
 gPad->SetLeftMargin(0.15);
 frame2->GetYaxis()->SetTitleOffset(1.4);
 frame2->Draw();
 c->cd(3);
 gPad->SetLeftMargin(0.15);
 frame3->GetYaxis()->SetTitleOffset(1.4);
 frame3->Draw();
 c->cd(4);
 gPad->SetLeftMargin(0.15);
 frame4->GetYaxis()->SetTitleOffset(1.4);
 frame4->Draw();
}

Hi @CS_p!

It doesn’t work because it’s not implemented. In the documentation for RooMCStudy::plotParam(), the optional FitGauss() argument does not appear.

What you can do is to write a function to fit a Gaussian based on the code that the RooMCStudy uses to fit the Gaussian in plotPull(), only more general:

void fitGaussToFrame(RooPlot& frame, RooDataSet const& fitParData)
{
   RooWorkspace ws;
   auto plotVar = frame.getPlotVar();
   const std::string plotVarName = plotVar->GetName();
   ws.import(*plotVar);
   std::stringstream expr;
   const double xMin = frame.GetXaxis()->GetXmin();
   const double xMax = frame.GetXaxis()->GetXmax();
   const double interval = xMax - xMin;
   expr << "Gaussian::frameGauss(" << plotVarName << ", frameMean[" << (0.5 * (xMax + xMin)) << ", " << xMin << ", "
        << xMax << "], frameSigma[" << (0.05 * interval) << ", " << (0.005 * interval) << ", " << (0.25 * interval)
        << "])";
   ws.factory(expr.str().c_str());

   RooRealVar& frameMean = *ws.var("frameMean");
   RooRealVar& frameSigma = *ws.var("frameSigma");
   RooAbsPdf& frameGauss = *ws.pdf("frameGauss");

   frameGauss.fitTo(const_cast<RooDataSet&>(fitParData), RooFit::Minos(0), RooFit::PrintLevel(-1)) ;
   frameGauss.plotOn(&frame) ;

   // Instead of using paramOn() without command arguments to plot the fit
   // parameters, we are building the parameter label ourselves for more
   // flexibility and pass this together with an appropriate layout
   // parametrization to paramOn().
   const int sigDigits = 2;
   const char * options = "ELU";
   std::stringstream ss;
   ss << "Fit parameters:\n"
      << "#mu: " << *std::unique_ptr<TString>{frameMean.format(sigDigits, options)}
      << "\n#sigma: " << *std::unique_ptr<TString>{frameSigma.format(sigDigits, options)};
   // We set the parameters constant to disable the default label. Still, we
   // use param() on as a wrapper for the text box generation.
   frameMean.setConstant(true);
   frameSigma.setConstant(true);
   frameGauss.paramOn(&frame, RooFit::Label(ss.str().c_str()), RooFit::Layout(0.60, 0.9, 0.9));
}

} // namespace

Then, you can include all your Gaussian fits like this:

 // Make plots of the distributions of mean, the error on mean and the pull of mean
 RooPlot *frame1 = mcstudy->plotParam(Nsig, Bins(40));
 fitGaussToFrame(*frame1, mcstudy->fitParDataSet());
 RooPlot *frame2 = mcstudy->plotError(Nsig, Bins(40));
 fitGaussToFrame(*frame2, mcstudy->fitParDataSet());
 RooPlot *frame3 = mcstudy->plotPull(Nsig, Bins(40), FitGauss(true));

 // Plot distribution of minimized likelihood
 RooPlot *frame4 = mcstudy->plotNLL(Bins(40));
 fitGaussToFrame(*frame4, mcstudy->fitParDataSet());

Here you have the full updated script:

I hope this helps!

Cheers,
Jonas

PS: If you think that it’s worth it if the other functions like plotParam() and plotError() are also supporting the FitGauss() argument, please feel free to open a ROOT improvement issue on GitHub, including the link to this forum post. Thanks!

Thanks @Jonas! This is working fine now.

I think it is worth including the gaussian fit for the parameter of interest (here, Nsig) in addition to its pull because we have to check the stability of the fitting by performing a linearity test i.e., Nsig (generated) vs Nsig (obtained). The Nsig (obtained) should be from Gaussian fit.

I will open the issue on GitHub.

Thanks

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.