ToyMC study for simultaneous fit

Roofit and RooStats
1

Dear Experts,

I am trying to generate 1000 toyMCs for simultaneous fit using one PDF and perform fit using another PDF, here is my sample code;

toy = simPdf_toy.generate(ROOT.RooArgList(Mbc,deltaE, sample_toy), totnum) 
##generate toys from simulateneous fit

here, simPdf_toy is the simultaneous fit PDF. Mbc and deltaE are the fitting variables and sample_toy is the RooCategory, defined as;

sample_toy = ROOT.RooCategory("sample_toy", "sample_toy")
sample_toy.defineType("mumu_toy")
sample_toy.defineType("ee_toy")

The sample_toy is included in the toyMC generation as instructed here, RooMCStudy with simultaneous fit - #3 by dennis. But when I perform the fit to the generated “toy” using another PDF, let’s say simPdf;

result = simPdf.fitTo(toy)

I get below errors;

[#0] ERROR:InputArguments -- RooTreeData::split(hmaster) ERROR category sample doesn't depend on any variable in this dataset
[#0] ERROR:Fitting -- RooAbsTestStatistic::initSimMode(nll_simPdf_hmaster) ERROR: index category of simultaneous pdf is missing in dataset, aborting
  Traceback (most recent call last):
   result = simPdf.fitTo(toy)
  TypeError: none of the 2 overloaded methods succeeded. Full details:
  RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, const RooLinkedList& cmdList) =>
  TypeError: takes at least 2 arguments (1 given)
  RooFitResult* RooAbsPdf::fitTo(RooAbsData& data, const RooCmdArg& arg1 =  RooCmdArg::none(), const RooCmdArg& arg2 = RooCmdArg::none(), const RooCmdArg& arg3 = RooCmdArg::none(), const RooCmdArg& arg4 = RooCmdArg::none(), const RooCmdArg& arg5 = RooCmdArg::none(), const RooCmdArg& arg6 = RooCmdArg::none(), const RooCmdArg& arg7 = RooCmdArg::none(), const RooCmdArg& arg8 = RooCmdArg::none()) =>
  runtime_error: RooAbsTestStatistic::initSimMode() ERROR, index category of simultaneous pdf is missing in the dataset, aborting

I didn’t find any example showing toyMC generation for simultaneous fit, so maybe the way I implemented the RooCategory in toy generation has some issues. Could anyone suggests to me, how should I generate the toys for simultaneous fit?

2

Hi,

Can you please post your full code (including the model) showing these problems, so we can investigate it.
Same ting for your other post on fitting the pull distribution

Lorenzo

3

Hi @moneta,

Thank you for your response. Please find my sample code below;

import ROOT


def toyMC():
    Mbc = ROOT.RooRealVar("Mbc", "M_{bc} [GeV/c^{2}]", 5.2, 5.29)
    deltaE = ROOT.RooRealVar("deltaE", "#DeltaE [GeV]", -0.15, 0.1)
    signum = 0
    nbkg_mumu = 553
    nbkg_ee = 387
    totnum = signum + nbkg_mumu + nbkg_ee

    # Mbc backgrund fit mumu channel
    a0_toy = ROOT.RooRealVar("a0_toy ", "a0_toy ", 0)
    bkg_mbc_mumu_toy = ROOT.RooArgusBG(
        "bkg_mbc_mumu_toy ", "bkg_mbc_mumu_toy ", Mbc, ROOT.RooFit.RooConst(5.289), a0_toy
    )
    # Mbc signal fit mumu channel
    mean_mbc_mumu_toy = ROOT.RooRealVar("mean_mbc_mumu_toy ", "mean_mbc_mumu_toy ", 5.27945)
    sigma_mbc_mumu_toy = ROOT.RooRealVar("sigma_mbc_mumu_toy ", "sigma_mbc_mumu_toy ", 0.00268)
    a_mbc_mumu_toy = ROOT.RooRealVar("a_mbc_mumu_toy ", "a_mbc_mumu_toy ", 2.02)
    n_mbc_mumu_toy = ROOT.RooRealVar("n_mbc_mumu_toy ", "n_mbc_mumu_toy ", 10.0)
    sig_mbc_mumu_toy = ROOT.RooCBShape(
        "sig_mbc_mumu_toy ",
        "sig_mbc_mumu_toy ",
        Mbc,
        mean_mbc_mumu_toy,
        sigma_mbc_mumu_toy,
        a_mbc_mumu_toy,
        n_mbc_mumu_toy,
    )
    # Delta E background fit mumu channel
    pol1_toy = ROOT.RooRealVar("pol1_toy ", "pol1_toy", -0.609)
    bkg_de_mumu_toy = ROOT.RooChebychev("bkg_de_mumu_toy", "bkg_de_mumu_toy", deltaE, ROOT.RooArgList(pol1_toy))
    # Delta E signal fit mumu channel
    mean_de_mumu_toy = ROOT.RooRealVar("mean_de_mumu_toy", "mean_de_mumu_toy", 0.012)
    lmbda_de_mumu_toy = ROOT.RooRealVar("lmbda_de_mumu_toy", "lmbda_de_mumu_toy", 0.0311)
    gamma_de_mumu_toy = ROOT.RooRealVar("gamma_de_mumu_toy", "gamma_de_mumu_toy", 0.96)
    delta_de_mumu_toy = ROOT.RooRealVar("delta_de_mumu_toy", "delta_de_mumu_toy", 1.41)
    johnson_mumu_toy = ROOT.RooJohnson(
        "johnson_mumu_toy",
        "johnson_mumu_toy",
        deltaE,
        mean_de_mumu_toy,
        lmbda_de_mumu_toy,
        gamma_de_mumu_toy,
        delta_de_mumu_toy,
    )
    sigma_de_mumu_toy = ROOT.RooRealVar("sigma_de_mumu_toy", "sigma_de_mumu_toy_toy", 0.0212)
    gauss_de_mumu_toy = ROOT.RooGaussian(
        "gauss_de_mumu_toy", "gauss_de_mumu_toy", deltaE, mean_de_mumu_toy, sigma_de_mumu_toy
    )
    frac_de_mumu_toy = ROOT.RooRealVar("frac_de_mumu_toy", "frac_de_mumu_toy", 0.75)
    sig_de_mumu_toy = ROOT.RooAddPdf(
        "sig_de_mumu_toy",
        "sig_de_mumu_toy",
        ROOT.RooArgList(johnson_mumu_toy, gauss_de_mumu_toy),
        ROOT.RooArgList(frac_de_mumu_toy),
    )
    sig_mumu_toy = ROOT.RooProdPdf("sig_mumu_toy", "sig_mumu_toy", ROOT.RooArgList(sig_mbc_mumu_toy, sig_de_mumu_toy))
    bkg_mumu_toy = ROOT.RooProdPdf("bkg_mumu_toy", "bkg_mumu_toy", ROOT.RooArgList(bkg_mbc_mumu_toy, bkg_de_mumu_toy))

    # Mbc backgrund fit ee channel
    bkg_mbc_ee_toy = ROOT.RooArgusBG("bkg_mbc_ee_toy", "bkg_mbc_ee_toy", Mbc, ROOT.RooFit.RooConst(5.289), a0_toy)
    # Mbc signal fit ee channel
    mean_mbc_ee_toy = ROOT.RooRealVar("mean_mbc_ee_toy", "mean_mbc_ee_toy", 5.27949)
    sigma_mbc_ee_toy = ROOT.RooRealVar("sigma_mbc_ee_toy", "sigma_mbc_ee_toy", 0.00271)
    a_mbc_ee_toy = ROOT.RooRealVar("a_mbc_ee_toy", "a_mbc_ee_toy", 1.71)
    n_mbc_ee_toy = ROOT.RooRealVar("n_mbc_ee_toy", "n_mbc_ee_toy", 10.0)
    sig_mbc_ee_toy = ROOT.RooCBShape(
        "sig_mbc_ee_toy", "sig_mbc_ee_toy", Mbc, mean_mbc_ee_toy, sigma_mbc_ee_toy, a_mbc_ee_toy, n_mbc_ee_toy
    )
    # Delta E background fit ee channel
    bkg_de_ee_toy = ROOT.RooChebychev("bkg_de_ee_toy", "bkg_de_ee_toy", deltaE, ROOT.RooArgList(pol1_toy))
    # Delta E signal fit ee channel
    mean_de_ee_toy = ROOT.RooRealVar("mean_de_ee_toy", "mean_de_ee_toy", 0.0171)
    lmbda_de_ee_toy = ROOT.RooRealVar("lmbda_de_ee_toy", "lmbda_de_ee_toy", 0.0336)
    gamma_de_ee_toy = ROOT.RooRealVar("gamma_de_ee_toy", "gamma_de_ee_toy", 1.09)
    delta_de_ee_toy = ROOT.RooRealVar("delta_de_ee_toy", "delta_de_ee_toy", 1.29)
    sig_de_ee_toy = ROOT.RooJohnson(
        "sig_de_ee_toy_toy", "sig_de_ee_toy", deltaE, mean_de_ee_toy, lmbda_de_ee_toy, gamma_de_ee_toy, delta_de_ee_toy
    )
    sig_mumu_toy = ROOT.RooProdPdf("sig_mumu_toy", "sig_mumu_toy", ROOT.RooArgList(sig_mbc_mumu_toy, sig_de_mumu_toy))
    bkg_mumu_toy = ROOT.RooProdPdf("bkg_mumu_toy", "bkg_mumu_toy", ROOT.RooArgList(bkg_mbc_mumu_toy, bkg_de_mumu_toy))
    sig_ee_toy = ROOT.RooProdPdf("sig_ee_toy", "sig_ee_toy", ROOT.RooArgList(sig_mbc_ee_toy, sig_de_ee_toy))
    bkg_ee_toy = ROOT.RooProdPdf("bkg_ee_toy", "bkg_ee_toy", ROOT.RooArgList(bkg_mbc_ee_toy, bkg_de_ee_toy))

    Nsig_toy = ROOT.RooRealVar("Nsig_toy", "Nsig_toy", signum)
    Nbkg_mumu_toy = ROOT.RooRealVar("Nbkg_mumu_toy", "Nbkg_mumu_toy", nbkg_mumu)
    Nbkg_ee_toy = ROOT.RooRealVar("Nbkg_ee_toy", "Nbkg_ee_toy", nbkg_ee)
    model_mumu_toy = ROOT.RooAddPdf(
        "model_mumu_toy",
        "model_mumu_toy",
        ROOT.RooArgList(sig_mumu_toy, bkg_mumu_toy),
        ROOT.RooArgList(Nsig_toy, Nbkg_mumu_toy),
    )
    model_ee_toy = ROOT.RooAddPdf(
        "model_ee_toy", "model_ee_toy", ROOT.RooArgList(sig_ee_toy, bkg_ee_toy), ROOT.RooArgList(Nsig_toy, Nbkg_ee_toy)
    )
    # toy generation
    sample_toy = ROOT.RooCategory("sample_toy", "sample_toy")
    sample_toy.defineType("mumu_toy")
    sample_toy.defineType("ee_toy")

    simPdf_toy = ROOT.RooSimultaneous("simPdf_toy", "simultaneous pdf", sample_toy)
    simPdf_toy.addPdf(model_mumu_toy, "mumu_toy")
    simPdf_toy.addPdf(model_ee_toy, "ee_toy")

    toy = simPdf_toy.generate(ROOT.RooArgList(Mbc, deltaE, sample_toy), totnum)

    # toyMC fitting
    # Mbc backgrund fit
    a0 = ROOT.RooRealVar("a0", "a0", -5, -100, 100)
    bkg_mbc_mumu = ROOT.RooArgusBG("bkg_mbc_mumu", "bkg_mbc_mumu", Mbc, ROOT.RooFit.RooConst(5.289), a0)
    # Mbc signal fit
    mean_mbc_mumu = ROOT.RooRealVar("mean_mbc_mumu", "mean_mbc_mumu", 5.27945)
    sigma_mbc_mumu = ROOT.RooRealVar("sigma_mbc_mumu", "sigma_mbc_mumu", 0.00268)
    a_mbc_mumu = ROOT.RooRealVar("a_mbc_mumu", "a_mbc_mumu", 2.02)
    n_mbc_mumu = ROOT.RooRealVar("n_mbc_mumu", "n_mbc_mumu", 10.0)
    sig_mbc_mumu = ROOT.RooCBShape(
        "sig_mbc_mumu", "sig_mbc_mumu", Mbc, mean_mbc_mumu, sigma_mbc_mumu, a_mbc_mumu, n_mbc_mumu
    )
    # Delta E background fit
    pol1 = ROOT.RooRealVar("pol1", "pol1", -0.23, -1, 1)
    bkg_de_mumu = ROOT.RooChebychev("bkg_de_mumu", "bkg_de_mumu", deltaE, ROOT.RooArgList(pol1))
    # Delta E signal fit
    mean_de_mumu = ROOT.RooRealVar("mean_de_mumu", "mean_de_mumu", 0.012)
    lmbda_de_mumu = ROOT.RooRealVar("lmbda_de_mumu", "lmbda_de_mumu", 0.0311)
    gamma_de_mumu = ROOT.RooRealVar("gamma_de_mumu", "gamma_de_mumu", 0.96)
    delta_de_mumu = ROOT.RooRealVar("delta_de_mumu", "delta_de_mumu", 1.41)
    johnson_mumu = ROOT.RooJohnson(
        "johnson_mumu", "johnson_mumu", deltaE, mean_de_mumu, lmbda_de_mumu, gamma_de_mumu, delta_de_mumu
    )
    sigma_de_mumu = ROOT.RooRealVar("sigma_de_mumu", "sigma_de_mumu", 0.0212)
    gauss_de_mumu = ROOT.RooGaussian("gauss_de_mumu", "gauss_de_mumu", deltaE, mean_de_mumu, sigma_de_mumu)
    frac_de_mumu = ROOT.RooRealVar("frac_de_mumu", "frac_de_mumu", 0.75)
    sig_de_mumu = ROOT.RooAddPdf(
        "sig_de_mumu", "sig_de_mumu", ROOT.RooArgList(johnson_mumu, gauss_de_mumu), ROOT.RooArgList(frac_de_mumu)
    )

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

    # model
    sig_mumu = ROOT.RooProdPdf("sig_mumu", "sig_mumu", ROOT.RooArgList(sig_mbc_mumu, sig_de_mumu))
    bkg_mumu = ROOT.RooProdPdf("bkg_mumu", "bkg_mumu", ROOT.RooArgList(bkg_mbc_mumu, bkg_de_mumu))

    sig_ee = ROOT.RooProdPdf("sig_ee", "sig_ee", ROOT.RooArgList(sig_mbc_ee, sig_de_ee))
    bkg_ee = ROOT.RooProdPdf("bkg_ee", "bkg_ee", ROOT.RooArgList(bkg_mbc_ee, bkg_de_ee))

    Nsig = ROOT.RooRealVar("Nsig", "Nsig", 2, 0, 7)
    Nbkg_mumu = ROOT.RooRealVar("Nbkg_mumu", "Nbkg_mumu", nbkg_mumu, 0, 1000)
    Nbkg_ee = ROOT.RooRealVar("Nbkg_ee", "Nbkg_ee", nbkg_ee, 0, 1000)
    model_mumu = ROOT.RooAddPdf(
        "model_mumu", "model_mumu", ROOT.RooArgList(sig_mumu, bkg_mumu), ROOT.RooArgList(Nsig, Nbkg_mumu)
    )
    model_ee = ROOT.RooAddPdf("model_ee", "model_ee", ROOT.RooArgList(sig_ee, bkg_ee), ROOT.RooArgList(Nsig, Nbkg_ee))

    sample = ROOT.RooCategory("sample", "sample")
    sample.defineType("mumu")
    sample.defineType("ee")

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

    # fitting
    result = simPdf.fitTo(toy)
    sta = result.status()
    pull = (Nsig.getVal() - signum) / Nsig.getError()
    nll = ROOT.RooNLLVar("nll", "nll", model, toy)
    print(Nsig, Nsig.getError(), Nbkg, Nbkg.getError(), pull, nll, sta)
    return [Nsig.getVal(), Nsig.getError(), Nbkg.getVal(), Nbkg.getError(), pull, nll.getVal(), sta]


if __name__ == "__main__":
    file = open("csvs/1000_toyStudy_0.csv", "w")
    for jj in range(1, 1001):
        print(red, jj, end)
        res = toyMC()
        print(res[0], res[1], res[2], res[3], res[4], res[5], res[6])
        file.write("%s %s %s %s %s %s %s \n" % (res[0], res[1], res[2], res[3], res[4], res[5], res[6]))
    file.close()
4

Hi!

The names of the observables in the toy dataset must match the names of the observables in the pdf that you are fitting to it, which includes the index category.

So instead of using this index category with the the _toy suffix in the name:

     sample_toy = ROOT.RooCategory("sample_toy", "sample_toy")
     sample_toy.defineType("mumu_toy")
     sample_toy.defineType("ee_toy")

Just use the same sample category that you also use in your fitting pdf:

     sample = ROOT.RooCategory("sample", "sample")
     sample.defineType("mumu")
     sample.defineType("ee")

Or was there any particular reason why you thought they should be different?

I hope this helps!
Jonas

5

Hi @jonas,

Thank you for your response. But, it doesn’t work when I use the same RooCategory for toy generation and fitting.
Here is the parameter correlation coefficient;

NO. GLOBAL 1 2 3 4 5 6
1 0.00000 1.000 0.000 0.000 0.000 0.000 0.000
2 0.00000 0.000 1.000 0.000 0.000 0.000 0.000
3 0.00000 0.000 0.000 1.000 0.000 0.000 0.000
4 0.00000 0.000 0.000 0.000 1.000 0.000 0.000
5 0.00000 0.000 0.000 0.000 0.000 1.000 0.000
6 0.00000 0.000 0.000 0.000 0.000 0.000 1.000

The fit didn’t converge, maybe the toy generation doesn’t work well with RooCategory!

6

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