doc/master/rf709__BarlowBeeston_8py.html

import math


import ROOT


def create_constraints(gammas, samples, global_observables):

    constraints = ROOT.RooArgSet()

    for i, gamma in enumerate(gammas):

        nom_name = "nom_" + gamma.GetName()

        tau_name = gamma.GetName() + "_tau"

        mean_name = gamma.GetName() + "_mean"

        poisson_name = gamma.GetName() + "_constraint"


        nom_val = sum(func.dataHist().weight(i) for func in samples)


        # Like in HistFactory, the gamma of an empty bin is fixed and not

        # constrained. Otherwise, it would be a flat direction of the likelihood.

        # The other gammas get the relative Poisson uncertainty as initial error.

        if nom_val == 0.0:

            gamma.setConstant(True)

            continue

        gamma.setError(1.0 / math.sqrt(nom_val))


        nominal = ROOT.RooRealVar(nom_name, nom_name, nom_val, 0.0, ROOT.RooNumber.infinity())

        nominal.setConstant(True)

        global_observables.add(nominal)


        tau = ROOT.RooConstVar(tau_name, tau_name, nominal.getVal())


        mean = ROOT.RooProduct(mean_name, mean_name, gamma, tau)


        poisson = ROOT.RooPoisson(poisson_name, poisson_name, nominal, mean)


        # The constraint pdf takes ownership of the objects it depends on, and

        # the returned set takes ownership of the constraint pdf.

        for helper in (nominal, tau, mean):

            ROOT.SetOwnership(helper, False)

        poisson.addOwnedComponents(ROOT.RooArgSet(nominal, tau, mean))


        constraints.addOwned(poisson)

    return constraints


ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.NumIntegration)


# First, construct a likelihood model with a Gaussian signal on top of a uniform background

x = ROOT.RooRealVar("x", "x", -20, 20)

x.setBins(25)


meanG = ROOT.RooRealVar("meanG", "meanG", 1, -10, 10)

sigG = ROOT.RooRealVar("sigG", "sigG", 1.5, -10, 10)

g = ROOT.RooGaussian("g", "Gauss", x, meanG, sigG)

u = ROOT.RooUniform("u", "Uniform", x)


# Generate the data to be fitted

sigData = g.generate(x, 50)

bkgData = u.generate(x, 1000)


sumData = ROOT.RooDataSet("sumData", "Gauss + Uniform", x)

sumData.append(sigData)

sumData.append(bkgData)


# Make histogram templates for signal and background.

# Let's take a signal distribution with low statistics and a more accurate

# background distribution.

# Normally, these come from Monte Carlo simulations, but we will just generate them.

dh_sig = g.generateBinned(x, 50)

dh_bkg = u.generateBinned(x, 10000)


# ***** Case 0 - 'Rigid templates' *****


# Construct histogram shapes for signal and background

p_h_sig = ROOT.RooHistFunc("p_h_sig", "p_h_sig", x, dh_sig)

p_h_bkg = ROOT.RooHistFunc("p_h_bkg", "p_h_bkg", x, dh_bkg)


# Construct scale factors for adding the two distributions

Asig0 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000)

Abkg0 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)


# Construct the sum model

model0 = ROOT.RooRealSumPdf("model0", "model0", [p_h_sig, p_h_bkg], [Asig0, Abkg0], True)


# ***** Case 1 - 'Barlow Beeston' *****


# Construct parameterized histogram shapes for signal and background

gammas_sig1 = ROOT.ParamHistFunc.createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10)

gammas_bkg1 = ROOT.ParamHistFunc.createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10)


p_phf_sig1 = ROOT.ParamHistFunc("p_phf_sig1", "p_phf_sig1", [x], gammas_sig1)

p_phf_bkg1 = ROOT.ParamHistFunc("p_phf_bkg1", "p_phf_bkg1", [x], gammas_bkg1)


p_p_sig1 = ROOT.RooProduct("p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1)

p_p_bkg1 = ROOT.RooProduct("p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1)


Asig1 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000)

Abkg1 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)


# Construct the sum of these

model1_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig1, p_p_bkg1], [Asig1, Abkg1], True)


# Construct the subsidiary poisson measurements constraining the histogram parameters

# These ensure that the bin contents of the histograms are only allowed to vary within

# the statistical uncertainty of the Monte Carlo.

globs1 = ROOT.RooArgSet()

constraints_sig1 = create_constraints(gammas_sig1, [p_h_sig], globs1)

constraints_bkg1 = create_constraints(gammas_bkg1, [p_h_bkg], globs1)


# Construct the joint model with template PDFs and constraints

model1 = ROOT.RooProdPdf("model1", "model1", [model1_tmp, *constraints_sig1, *constraints_bkg1])


# ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) *****


# Construct the histogram shapes, using the same parameters for signal and background.

# There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to

# fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background

# histogram, so the bins fluctuate synchronously. This reduces the number of parameters.

gammas2 = ROOT.ParamHistFunc.createParamSet("p_ph_gamma", x.numBins(), 0, 10)


p_phf2 = ROOT.ParamHistFunc("p_phf2", "p_phf2", [x], gammas2)


p_p_sig2 = ROOT.RooProduct("p_p_sig2", "p_p_sig2", p_h_sig, p_phf2)

p_p_bkg2 = ROOT.RooProduct("p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2)


Asig2 = ROOT.RooRealVar("Asig", "Asig", 1, 0.01, 5000)

Abkg2 = ROOT.RooRealVar("Abkg", "Abkg", 1, 0.01, 5000)


# As before, construct the sum of signal2 and background2

model2_tmp = ROOT.RooRealSumPdf("sp_ph", "sp_ph", [p_p_sig2, p_p_bkg2], [Asig2, Abkg2], True)


# Construct the subsidiary poisson measurements constraining the statistical fluctuations

globs2 = ROOT.RooArgSet()

constraints2 = create_constraints(gammas2, [p_h_sig, p_h_bkg], globs2)


# Construct the joint model

model2 = ROOT.RooProdPdf("model2", "model2", [model2_tmp, *constraints2])


# ************ Fit all models to data and plot *********************


result0 = model0.fitTo(sumData, PrintLevel=0, Save=True)

result1 = model1.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs1)

result2 = model2.fitTo(sumData, PrintLevel=0, Save=True, GlobalObservables=globs2)


can = ROOT.TCanvas("can", "", 1500, 600)

can.Divide(3, 1)


pt = ROOT.TPaveText(-19.5, 1, -2, 25)

pt.SetFillStyle(0)

pt.SetBorderSize(0)


can.cd(1)

frame = x.frame(Title="No template uncertainties")

# Plot data to enable automatic determination of model0 normalisation:

sumData.plotOn(frame)

model0.plotOn(frame, LineColor="b", VisualizeError=result0)

# Plot data again to show it on top of model0 error bands:

sumData.plotOn(frame)

# Plot model components

model0.plotOn(frame, LineColor="b")

model0.plotOn(frame, Components=[p_h_sig], LineColor="kAzure")

model0.plotOn(frame, Components=[p_h_bkg], LineColor="r")

model0.paramOn(frame)


sigData.plotOn(frame, MarkerColor="b")

frame.Draw()


for text in [

    "No template uncertainties",

    "are taken into account.",

    "This leads to low errors",

    "for the parameters A, since",

    "the only source of errors",

    "are the data statistics.",

]:

    pt.AddText(text)

pt.DrawClone()


can.cd(2)

frame = x.frame(Title="Barlow Beeston for Sig & Bkg separately")

sumData.plotOn(frame)

model1.plotOn(frame, LineColor="b", VisualizeError=result1)

# Plot data again to show it on top of error bands:

sumData.plotOn(frame)

model1.plotOn(frame, LineColor="b")

model1.plotOn(frame, Components=[p_p_sig1], LineColor="kAzure")

model1.plotOn(frame, Components=[p_p_bkg1], LineColor="r")

model1.paramOn(frame, Parameters=[Asig1, Abkg1])


sigData.plotOn(frame, MarkerColor="b")

frame.Draw()


pt.Clear()

for text in [

    "With gamma parameters, the",

    "signal & background templates",

    "can adapt to the data.",

    "Note how the blue signal",

    "template changes its shape.",

    "This leads to higher errors",

    "of the scale parameters A.",

]:

    pt.AddText(text)

pt.DrawClone()


can.cd(3)

frame = x.frame(Title="Barlow Beeston light for (Sig+Bkg)")

sumData.plotOn(frame)

model2.plotOn(frame, LineColor="b", VisualizeError=result2)

# Plot data again to show it on top of model0 error bands:

sumData.plotOn(frame)

model2.plotOn(frame, LineColor="b")

model2.plotOn(frame, Components=[p_p_sig2], LineColor="kAzure")

model2.plotOn(frame, Components=[p_p_bkg2], LineColor="r")

model2.paramOn(frame, Parameters=[Asig2, Abkg2])


sigData.plotOn(frame, MarkerColor="b")

frame.Draw()


pt.Clear()

for text in [

    "When signal and background",

    "template share one gamma para-",

    "meter per bin, they adapt less.",

    "The errors of the A parameters",

    "also shrink slightly.",

]:

    pt.AddText(text)

pt.DrawClone()


print(f"Asig [normal ] = {Asig0.getVal()} +/- {Asig0.getError()}")

print(f"Asig [BB     ] = {Asig1.getVal()} +/- {Asig1.getError()}")

print(f"Asig [BBlight] = {Asig2.getVal()} +/- {Asig2.getError()}")


can.SaveAs("rf709_BarlowBeeston.png")

sum
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2335