How to Preserve Uncertainties When Creating TH1?

[feickert]

In RooStats what is the proper way to save a RooDataHist to a TH1 such that the uncertainties on the RooDataHist are preserved? When I use RooDataHist::createHistogram() the resulting TH1 has unreasonable uncertainties on it.

The issue “RooFit: saving RooDataHist as ROOT TH1 histogram?” touches on how to get a TH1, but doesn’t give a solution to my problem.

Below is the example PyROOT code demonstrating the problem with the plots showing the problem attached.

import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)

# A one bin model with flat signal and background samples
n_events = 1000
n_bins = 1
range_low = 0
range_high = 10

frac_sig = 0.1
frac_bkg = 1 - frac_sig

# create input histograms with default Poisson uncertainties
signal_hist = ROOT.TH1F("signal_hist", "signal_hist",
                        n_bins, range_low, range_high)
signal_hist.SetBinContent(1, frac_sig * n_events)

background_hist = ROOT.TH1F(
    "background_hist", "background_hist", n_bins, range_low, range_high)
background_hist.SetBinContent(1, frac_bkg * n_events)

data_hist = ROOT.TH1F("data_hist", "data_hist", n_bins, range_low, range_high)
data_hist.SetBinContent(1, n_events)

# Create the measurement
meas = ROOT.RooStats.HistFactory.Measurement("meas", "meas")
ROOT.SetOwnership(meas, False) #c.f. https://root-forum.cern.ch/t/ownership-in-python/9573

meas.SetPOI("SigXsecOverSM")
meas.AddConstantParam("Lumi")
meas.SetLumi(1.0)
meas.SetLumiRelErr(0.10)
meas.SetExportOnly(False)

# Create a channel
SR = ROOT.RooStats.HistFactory.Channel("SR")
ROOT.SetOwnership(SR, False)
SR.SetData(data_hist)

# Create the signal sample
signal = ROOT.RooStats.HistFactory.Sample("signal")
ROOT.SetOwnership(signal, False)
signal.SetNormalizeByTheory(False)
signal.SetHisto(signal_hist)
signal.AddNormFactor("SigXsecOverSM", 1, 0, 3)
# Add normalisation systematic uncertainty of 2%
signal.AddOverallSys("signal_uncertainty", 1.0 - 0.02, 1.0 + 0.02)

# Create the background sample
background = ROOT.RooStats.HistFactory.Sample("background")
ROOT.SetOwnership(background, False)
background.SetNormalizeByTheory(False)
background.SetHisto(background_hist)
# Add normalisation systematic uncertainty of 2%
background.AddOverallSys("background_uncertainty", 1.0 - 0.02, 1.0 + 0.02)

# Add the samples to the channel
SR.AddSample(signal)
SR.AddSample(background)

# Add the channel to the measurement
meas.AddChannel(SR)

# make the factory
factory = ROOT.RooStats.HistFactory.HistoToWorkspaceFactoryFast()
ws = factory.MakeSingleChannelModel(meas, SR)
ws.SetName('HistFactoryWorkspace')

model = ws.pdf('simPdf')  # 'simPdf' is the default name assigned

c1 = ROOT.TCanvas()
frame = ws.var("obs_x_SR").frame()
print(ws.var("obs_x_SR"))
ws.data("obsData").plotOn(frame)
frame.Draw()
c1.Draw()
c1.SaveAs("data_RooPlot.pdf")
print(ws.data("obsData"))

roo_data_hist = ROOT.RooDataHist("h", "h", ROOT.RooArgSet(
    ws.var("obs_x_SR")), ws.data("obsData"))
h_chan_data = roo_data_hist.createHistogram("h_chan_data", ws.var("obs_x_SR"))
print(roo_data_hist)
print(h_chan_data)
roo_data_hist.Print("v")
c2 = ROOT.TCanvas()
h_chan_data.Draw()
c2.Draw()
c2.SaveAs("data_TH1.pdf")

example.py (2.6 KB)

[lukasheinrich]

The issue seems to come from trying to use fillHistogram on a RooDataSet that is weighted (which is how HF creates its data objects). I’m attaching a pure RooFit minimal working example:

debug.py (911 Bytes)

as a workaround one can call Sumw2(0) on the histogram after the fill to rebuild the errors.

[feickert]

As an example of how Lukas’s suggestion fixes things, I’ve included the fixed version of the code and the resultant plot:

import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)

# A one bin model with flat signal and background samples
n_events = 1000
n_bins = 1
range_low = 0
range_high = 10

frac_sig = 0.1
frac_bkg = 1 - frac_sig

# create input histograms with default Poisson uncertainties
signal_hist = ROOT.TH1F("signal_hist", "signal_hist",
                        n_bins, range_low, range_high)
signal_hist.SetBinContent(1, frac_sig * n_events)

background_hist = ROOT.TH1F(
    "background_hist", "background_hist", n_bins, range_low, range_high)
background_hist.SetBinContent(1, frac_bkg * n_events)

data_hist = ROOT.TH1F("data_hist", "data_hist", n_bins, range_low, range_high)
data_hist.SetBinContent(1, n_events)

# Create the measurement
meas = ROOT.RooStats.HistFactory.Measurement("meas", "meas")
ROOT.SetOwnership(meas, False)

meas.SetPOI("SigXsecOverSM")
meas.AddConstantParam("Lumi")
meas.SetLumi(1.0)
meas.SetLumiRelErr(0.10)
meas.SetExportOnly(False)

# Create a channel
SR = ROOT.RooStats.HistFactory.Channel("SR")
ROOT.SetOwnership(SR, False)
SR.SetData(data_hist)

# Create the signal sample
signal = ROOT.RooStats.HistFactory.Sample("signal")
ROOT.SetOwnership(signal, False)
signal.SetNormalizeByTheory(False)
signal.SetHisto(signal_hist)
signal.AddNormFactor("SigXsecOverSM", 1, 0, 3)
# Add normalisation systematic uncertainty of 2%
signal.AddOverallSys("signal_uncertainty", 1.0 - 0.02, 1.0 + 0.02)

# Create the background sample
background = ROOT.RooStats.HistFactory.Sample("background")
ROOT.SetOwnership(background, False)
background.SetNormalizeByTheory(False)
background.SetHisto(background_hist)
# Add normalisation systematic uncertainty of 2%
background.AddOverallSys("background_uncertainty", 1.0 - 0.02, 1.0 + 0.02)

# Add the samples to the channel
SR.AddSample(signal)
SR.AddSample(background)

# Add the channel to the measurement
meas.AddChannel(SR)

# make the factory
factory = ROOT.RooStats.HistFactory.HistoToWorkspaceFactoryFast()
ws = factory.MakeSingleChannelModel(meas, SR)
ws.SetName('HistFactoryWorkspace')

model = ws.pdf('simPdf')  # 'simPdf' is the default name assigned

c1 = ROOT.TCanvas()
frame = ws.var("obs_x_SR").frame()
print(ws.var("obs_x_SR"))
ws.data("obsData").plotOn(frame)
frame.Draw()
c1.Draw()
c1.SaveAs("data_RooPlot.pdf")
print(ws.data("obsData"))

roo_data_hist = ROOT.RooDataHist("h", "h", ROOT.RooArgSet(
    ws.var("obs_x_SR")), ws.data("obsData"))
h_chan_data = roo_data_hist.createHistogram("h_chan_data", ws.var("obs_x_SR"))
# The fix is here
h_chan_data.Sumw2(0) # clear the uncertainties
h_chan_data.Sumw2() # rebuild the uncertainty
h_chan_data.SetMinimum(0)
h_chan_data.SetMaximum(h_chan_data.GetBinContent(1) * 1.1)
#
print(roo_data_hist)
print(h_chan_data)
roo_data_hist.Print("v")
c2 = ROOT.TCanvas()
h_chan_data.Draw()
c2.Draw()
c2.SaveAs("data_TH1.pdf")

fixed_example.py (2.7 KB)