PyRoot: Workspace embedded data sets disappear

Hello,

I have been trying to conver some code to pyroot in order to better display what it is doing.

However, I have an issue that when a workspace is read in in pyroot using:

tfile=ROOT.TFile("workspace_error.root","READ")
workspace=ROOT.RooWorkspace(tfile.Get("w"))
tfile.Close()

It seems all the embedded data sets ( that are there in root terminal browser) when it is created are lost when I read it back in. This means that fits fail. Further, if I save the histogram again the embedded datasets are lost and cannot be found in root terminal browser.

I cannot share my original program but have attached an example the reproduces the issue, however I cannot seem to upload it so apologies for the block of code:

import ROOT
import numpy as np
import numpy.random as rnd
from scipy.stats import norm
from scipy.stats import expon as exp





def create_mock_workspace():
    #need to create and save some mock workspace with:
    #-signal model
    #-fake background with it's own fake data
    #-fake data
    #-fit the total model to data
    #generate a pdf
    #need to inverse background pdf
    #grab 1000 random samples from the background 
    bkg_data=norm.ppf(rnd.rand(1000),loc=1)
    #create a datahist
    background_hist=ROOT.TH1F("bkg_hist","histogram of background data",100,0,2)
    for b in bkg_data:
        background_hist.Fill(b)
    #create signal samples
    sig_data=exp.ppf(rnd.rand(1000))
    sig_hist=ROOT.TH1F("sig_hist","histogram of sig data",100,0,2)
    for s in sig_data:
        sig_hist.Fill(s)
    #----Now start the workspace
    workspace=ROOT.RooWorkspace("w")
    x=ROOT.RooRealVar("x","x[cm]",0,2)
    x.setBins(100)
    workspace.Import(x)
    sig_amplitude=ROOT.RooRealVar("sig_amp","signal amplitude", 1)
    bkg_amplitude=ROOT.RooRealVar("bkg_amp","background amplitude", 0.8)
    workspace.defineSet("obs", "x")
    workspace.Import(sig_amplitude)
    workspace.Import(bkg_amplitude)
    #--- need to add the signal models--
    bkg_data_hist=ROOT.RooDataHist("data_hist_bkg","bkg RooDataHist",workspace.set("obs"),background_hist)
    sig_data_hist=ROOT.RooDataHist("data_hist_bkg","bkg RooDataHist",workspace.set("obs"),background_hist)
    bkg_pdf=ROOT.RooHistPdf("pdf_bkg","background pdf",workspace.set("obs"),bkg_data_hist,1)
    sig_pdf=ROOT.RooHistPdf("pdf_sig","signal pdf",workspace.set("obs"),sig_data_hist,1)
    workspace.Import(bkg_pdf)
    workspace.Import(sig_pdf)
    #----create the event model----
    workspace.factory("SUM::EventModel(sig_amp*pdf_sig,bkg_amp*pdf_bkg)")
    workspace.factory("Gaussian::constraint_bkg(bkg_amp,0.8,0.2)")
    workspace.factory("PROD::TotalEventModel(EventModel,constraint_bkg)")
    #-----create signal model--------
    workspace.defineSet("poi", "sig_amp")
    workspace.defineSet("nuis", "bkg_amp")
    sig_bg_model = ROOT.RooStats.ModelConfig("SB_model")
    sig_bg_model.SetName("SB_model")
    sig_bg_model.SetWorkspace(workspace)
    sig_bg_model.SetPdf(workspace.pdf("TotalEventModel"))
    sig_bg_model.SetObservables(workspace.set("obs"))
    sig_bg_model.SetParametersOfInterest(workspace.set("poi"))
    sig_bg_model.SetNuisanceParameters(workspace.set("nuis"))
    workspace.Import(sig_bg_model)
    #save the workspace
    tfile=ROOT.TFile("workspace_error.root","RECREATE")
    workspace.Write()
    print("===================workspace after being created==================================")
    workspace.Print()
    tfile.Close()

def perform_fit():
    #create pseudo data from bkg 
    pseudo_bkg_data=norm.ppf(rnd.rand(100),loc=1)
    pseudo_hist=ROOT.TH1F("obs_data","pseudo observed data",100,0,2)
    for p in pseudo_bkg_data:
        pseudo_hist.Fill(p)
    tfile=ROOT.TFile("workspace_error.root","READ")
    workspace=ROOT.RooWorkspace(tfile.Get("w"))
    tfile.Close()
    print("===================workspace after being read in==================================")
    workspace.Print()
    #get the model
    signal_background_model=ROOT.RooStats.ModelConfig(workspace.obj("SB_model"))
    #redefine the data
    obs_data=ROOT.RooDataHist("data_hist_pseudo", "pseudo RooDataHist", workspace.set("obs"), pseudo_hist)
    #set up arguments
    model_nuisance_parameters=signal_background_model.GetNuisanceParameters()
    arglist = ROOT.RooLinkedList()
    arglist.Add(ROOT.RooFit.Extended(False))
    arglist.Add(ROOT.RooFit.Constrain(model_nuisance_parameters)) #I've used a seperate variable
    #for this rather than for getting the prefit values -> not sure if necessary!
    arglist.Add( ROOT.RooFit.Offset(True))
    arglist.Add( ROOT.RooFit.Strategy(1))
    arglist.Add( ROOT.RooFit.Save(True))
    arglist.Add( ROOT.RooFit.SumW2Error(False))
    arglist.Add(ROOT.RooFit.Verbose(False))
    arglist.Add(ROOT.RooFit.PrintLevel(-1))
    arglist.Add(ROOT.RooFit.Warnings(True))
    signal_background_model.GetPdf().fitTo(obs_data,arglist)
    tfile=ROOT.TFile("fitted_workspace_error.root","RECREATE")
    workspace.Write()
    print("===================workspace after being fitted to fake data==================================")
    workspace.Print()
    tfile.Close()


def main():
    print("Running Main")
    create_mock_workspace()
    perform_fit()
    tfile=ROOT.TFile("fitted_workspace_error.root","READ")
    workspace=ROOT.RooWorkspace(tfile.Get("w"))
    print("===================workspace after being fitted and reread==================================")
    workspace.Print()
    tfile.Close()

    return 0

main()

which yields terminal output:

user> python reproduce_workspace_error.py 
Running Main

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::x
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::sig_amp
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::bkg_amp
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing dataset data_hist_bkg
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooHistPdf::pdf_bkg
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooHistPdf::pdf_sig
===================workspace after being created==================================

RooWorkspace(w) w contents

variables
---------
(bkg_amp,sig_amp,x)

p.d.f.s
-------
RooAddPdf::EventModel[ sig_amp * pdf_sig + bkg_amp * pdf_bkg ] = 575
RooProdPdf::TotalEventModel[ EventModel * constraint_bkg ] = 575
RooGaussian::constraint_bkg[ x=bkg_amp mean=0.8 sigma=0.2 ] = 1
RooHistPdf::pdf_bkg[ pdfObs=(x) ] = 575
RooHistPdf::pdf_sig[ pdfObs=(x) ] = 575

embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::data_hist_bkg(x)

named sets
----------
SB_model_NuisParams:(bkg_amp)
SB_model_Observables:(x)
SB_model_POI:(sig_amp)
nuis:(bkg_amp)
obs:(x)
poi:(sig_amp)

generic objects
---------------
RooStats::ModelConfig::SB_model

===================workspace after being read in==================================

RooWorkspace(w) w contents

variables
---------
(bkg_amp,sig_amp,x)

p.d.f.s
-------
RooAddPdf::EventModel[ sig_amp * pdf_sig + bkg_amp * pdf_bkg ] = 575
RooProdPdf::TotalEventModel[ EventModel * constraint_bkg ] = 575
RooGaussian::constraint_bkg[ x=bkg_amp mean=0.8 sigma=0.2 ] = 1
RooHistPdf::pdf_bkg[ pdfObs=(x) ] = 575
RooHistPdf::pdf_sig[ pdfObs=(x) ] = 575

named sets
----------
SB_model_NuisParams:(bkg_amp)
SB_model_Observables:(x)
SB_model_POI:(sig_amp)
nuis:(bkg_amp)
obs:(x)
poi:(sig_amp)

generic objects
---------------
RooStats::ModelConfig::SB_model

[#1] INFO:Minization --  Including the following constraint terms in minimization: (constraint_bkg)
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_TotalEventModel_data_hist_pseudo_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --  The following expressions have been identified as constant and will be precalculated and cached: (TotalEventModel)
Error in <ROOT::Math::Fitter::SetFCN>: FCN function has zero parameters 
[#0] WARNING:Minization -- RooMinimizer::hesse: Error, run Migrad before Hesse!
[#0] WARNING:Minization -- RooMinimizer::save: Error, run minimization before!
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
===================workspace after being fitted to fake data==================================

RooWorkspace(w) w contents

variables
---------
(bkg_amp,sig_amp,x)

p.d.f.s
-------
RooAddPdf::EventModel[ sig_amp * pdf_sig + bkg_amp * pdf_bkg ] = 575
RooProdPdf::TotalEventModel[ EventModel * constraint_bkg ] = 575
RooGaussian::constraint_bkg[ x=bkg_amp mean=0.8 sigma=0.2 ] = 1
RooHistPdf::pdf_bkg[ pdfObs=(x) ] = 575
RooHistPdf::pdf_sig[ pdfObs=(x) ] = 575

named sets
----------
SB_model_NuisParams:(bkg_amp)
SB_model_Observables:(x)
SB_model_POI:(sig_amp)
nuis:(bkg_amp)
obs:(x)
poi:(sig_amp)

generic objects
---------------
RooStats::ModelConfig::SB_model

===================workspace after being fitted and reread==================================

RooWorkspace(w) w contents

variables
---------
(bkg_amp,sig_amp,x)

p.d.f.s
-------
RooAddPdf::EventModel[ sig_amp * pdf_sig + bkg_amp * pdf_bkg ] = 575
RooProdPdf::TotalEventModel[ EventModel * constraint_bkg ] = 575
RooGaussian::constraint_bkg[ x=bkg_amp mean=0.8 sigma=0.2 ] = 1
RooHistPdf::pdf_bkg[ pdfObs=(x) ] = 575
RooHistPdf::pdf_sig[ pdfObs=(x) ] = 575

named sets
----------
SB_model_NuisParams:(bkg_amp)
SB_model_Observables:(x)
SB_model_POI:(sig_amp)
nuis:(bkg_amp)
obs:(x)
poi:(sig_amp)

generic objects
---------------
RooStats::ModelConfig::SB_model

user> 

Hi @Joshua_Green,

It seems that the ownership of the workspace (and its child objects) is transferred to gROOT in the first tfile.Close() in perform_fit(), but embedded datasets are not written to the newly opened TFile. Probably @jonas can help you better with your issue.

Cheers,
J.

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