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>