In [1]:
import csv
import numpy as np
import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)
%jsroot on

Welcome to JupyROOT 6.07/07

Getting Kaggle Higgs Machine Learning Data

The data has three categories of entries: Derived variables (those that have a physics meaning), principle variables (those that are used to create the derived variables), and event info comprised of event weights, and classification (signal or background).

In [2]:
inputfile = 'input/training.csv'
all = list(csv.reader(open(inputfile,"rb"), delimiter=','))
headers = all[0]
print 'using the derived variables:',headers[1:14]

print ''

print 'or the principle variables:',headers[14:-2]

print ''

print 'with weights and classification at:',headers[-2:]

using the derived variables: ['DER_mass_MMC', 'DER_mass_transverse_met_lep', 'DER_mass_vis', 'DER_pt_h', 'DER_deltaeta_jet_jet', 'DER_mass_jet_jet', 'DER_prodeta_jet_jet', 'DER_deltar_tau_lep', 'DER_pt_tot', 'DER_sum_pt', 'DER_pt_ratio_lep_tau', 'DER_met_phi_centrality', 'DER_lep_eta_centrality']

or the principle variables: ['PRI_tau_pt', 'PRI_tau_eta', 'PRI_tau_phi', 'PRI_lep_pt', 'PRI_lep_eta', 'PRI_lep_phi', 'PRI_met', 'PRI_met_phi', 'PRI_met_sumet', 'PRI_jet_num', 'PRI_jet_leading_pt', 'PRI_jet_leading_eta', 'PRI_jet_leading_phi', 'PRI_jet_subleading_pt', 'PRI_jet_subleading_eta', 'PRI_jet_subleading_phi', 'PRI_jet_all_pt']

with weights and classification at: ['Weight', 'Label']

In this example I choose to look at the MMC mass since it gives the best separation of signal and background.

I also scale the signal by 10 such that it is visible in the plots.

In [3]:
all = list(csv.reader(open(inputfile,"rb"), delimiter=','))
VBF = [variables for variables in all[1:] if (float(variables[6]) > 200 and float(variables[5]) > 2) ]
VBF_Loose = [variables for variables in all[1:] if (float(variables[6]) > 180 and float(variables[5]) > 1.8) ]
VBF_Tight = [variables for variables in all[1:] if (float(variables[6]) > 220 and float(variables[5]) > 2.2) ]

Creating a workspace

In order to perform statistical tests on this code using the RooStats suite of libraries, first we need to create a workspace. A workspace should contain the statistical model (including any defining parameters such as the normalisations), a pdf of the expected distribution, a parameter of interest (as included in the model), and data to compare this model to.

The first step is to move these distributions from vectors to histograms that can be read by RooFit.

Higgs Enhanced SR

Let's now define a signal region which is just this MMC mass between 100 and 150 GeV

In [4]:
SR = [variables for variables in VBF[1:] if (float(variables[1]) > 100 and float(variables[1]) < 150) ]
SR_Loose = [variables for variables in VBF_Loose[1:] if (float(variables[1]) > 100 and float(variables[1]) < 150) ]
SR_Tight = [variables for variables in VBF_Tight[1:] if (float(variables[1]) > 100 and float(variables[1]) < 150) ]
In [5]:
def get_var(region, region_name=None, h_params=[5,100,150]):
    signal = [float(row[1]) for row in region[1:] if row[-1] == 's' ]
    sigweights = [float(row[-2])*10 for row in region[1:] if row[-1] == 's' ] 
    background = [float(row[1]) for row in region[1:] if row[-1] == 'b' ]
    bkgweights = [float(row[-2]) for row in region[1:] if row[-1] == 'b' ]
    h3 = ROOT.TH1D("{} data".format(region_name),"data histogram",h_params[0],h_params[1],h_params[2])
    h1 = ROOT.TH1D("{} signal".format(region_name),"signal histogram",h_params[0],h_params[1],h_params[2])
    for i in range(len(signal)):     
        h1.Fill(signal[i],sigweights[i])
        h3.Fill(signal[i],sigweights[i])
    h2 = ROOT.TH1D("{} background".format(region_name),"background histogram",h_params[0],h_params[1],h_params[2])
    for i in range(len(background)):     
        h2.Fill(background[i],bkgweights[i])
        h3.Fill(background[i],bkgweights[i])
    return h3, h1, h2, sum(sigweights), sum(bkgweights)
In [6]:
data_hist, sig_hist, bkg_hist, ns, nb = get_var(SR,"SR", [5,100,150])
data_up_hist, sig_up_hist, bkg_up_hist, nsup, nbup = get_var(SR_Loose,"SR_Loose", [5,100,150])
data_down_hist, sig_down_hist, bkg_down_hist, nsdown, nbdown = get_var(SR_Tight,"SR_Tight", [5,100,150])

print 'n signal {0}(+{1}%,{2}%)'.format(ns,nsup/ns*100-100,nsdown/ns*100-100)
print 'n background {0}(+{1}%,{2}%)'.format(nb,nbup/nb*100-100,nbdown/nb*100-100)

n signal 421.876963784(+7.59966405052%,-7.06385384704%)
n background 1819.96369349(+20.7978888549%,-16.2194878628%)

Create Histfactory Measurement

First we set the Parameter of interest, and several constant parameters.

In [7]:
meas = ROOT.RooStats.HistFactory.Measurement("meas", "meas")
meas.SetPOI( "SigXsecOverSM" )

We don't include a specific Luminocity here, but since HistFactory gives it to us for free, we set it constant with a dummy value of 2% (as fairly standard for events with ATLAS)

In [8]:
meas.SetLumi( 1.0 )
meas.SetLumiRelErr( 0.02 )
meas.AddConstantParam("Lumi")

Create a channel and set the measured value of data (no extenal histogram is necessarry for cut-and-count)

Here we use an, 'asimov' dataset which is the nominal signal+background model.

In [9]:
chan = ROOT.RooStats.HistFactory.Channel( "SR" )
chan.SetStatErrorConfig(0.05, "Poisson")
chan.SetData( data_hist )

Create the signal sample and set its histogram

In [10]:
signal = ROOT.RooStats.HistFactory.Sample( "signal" )
signal.SetNormalizeByTheory( False )
signal.SetHisto( sig_hist )

Add the parmaeter of interest and a systematic and try to make intelligent choice of upper bound

In [11]:
signal.AddNormFactor( "SigXsecOverSM", 1, 0, 3)

Since we have a signal normalisation uncertainty, we can add it too

In [12]:
uncertainty_up   = nsup/ns
uncertainty_down = nsdown/ns
signal.AddOverallSys( "signal_norm_uncertainty",  uncertainty_down, uncertainty_up )

and a signal shape uncertainty

In [13]:
signal_shape = ROOT.RooStats.HistFactory.HistoSys("signal_shape")
signal_shape.SetHistoHigh( sig_up_hist )
signal_shape.SetHistoLow( sig_down_hist )
signal.AddHistoSys( signal_shape )

Finally, we add this sample to the channel

In [14]:
chan.AddSample( signal )

Create the background sample and set its histogram

In [15]:
background = ROOT.RooStats.HistFactory.Sample( "background" )
background.SetNormalizeByTheory( False )
background.SetHisto( bkg_hist )
In [16]:
uncertainty_up   = nbup/nb
uncertainty_down = nbdown/nb
background.AddOverallSys( "background_uncertainty",  uncertainty_down, uncertainty_up )
background_shape = ROOT.RooStats.HistFactory.HistoSys("background_shape")
background_shape.SetHistoHigh( bkg_up_hist )
background_shape.SetHistoLow( bkg_down_hist )
background.AddHistoSys( background_shape )
In [17]:
chan.AddSample( background )
In [18]:
meas.AddChannel(chan)

Make workspace!

In [19]:
hist2workspace = ROOT.RooStats.HistFactory.HistoToWorkspaceFactoryFast(meas)
In [20]:
workspace = hist2workspace.MakeSingleChannelModel(meas, chan)


-------------------
Starting to process SR channel with 1 observables
lumi str = [1,0,10]
lumi Error str = nominalLumi[1,0,1.2],0.02
making normFactor: SigXsecOverSM
Gaussian::alpha_signal_shapeConstraint(alpha_signal_shape,nom_alpha_signal_shape[0.,-10,10],1.)
Gaussian::alpha_background_shapeConstraint(alpha_background_shape,nom_alpha_background_shape[0.,-10,10],1.)
-----------------------------------------
import model into workspace
RooDataSet::AsimovData[obs_x_SR,weight:binWeightAsimov] = 5 entries (2241.84 weighted)

RooWorkspace(SR) SR workspace contents

variables
---------
(Lumi,SigXsecOverSM,alpha_background_shape,alpha_background_uncertainty,alpha_signal_norm_uncertainty,alpha_signal_shape,binWidth_obs_x_SR_0,binWidth_obs_x_SR_1,nom_alpha_background_shape,nom_alpha_background_uncertainty,nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nominalLumi,obs_x_SR,weightVar)

p.d.f.s
-------
RooRealSumPdf::SR_model[ binWidth_obs_x_SR_0 * L_x_signal_SR_overallSyst_x_HistSyst + binWidth_obs_x_SR_1 * L_x_background_SR_overallSyst_x_HistSyst ] = 26.5578/2241.84
RooGaussian::alpha_background_shapeConstraint[ x=alpha_background_shape mean=nom_alpha_background_shape sigma=1 ] = 1
RooGaussian::alpha_background_uncertaintyConstraint[ x=alpha_background_uncertainty mean=nom_alpha_background_uncertainty sigma=1 ] = 1
RooGaussian::alpha_signal_norm_uncertaintyConstraint[ x=alpha_signal_norm_uncertainty mean=nom_alpha_signal_norm_uncertainty sigma=1 ] = 1
RooGaussian::alpha_signal_shapeConstraint[ x=alpha_signal_shape mean=nom_alpha_signal_shape sigma=1 ] = 1
RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.02 ] = 1
RooProdPdf::model_SR[ lumiConstraint * alpha_signal_norm_uncertaintyConstraint * alpha_signal_shapeConstraint * alpha_background_uncertaintyConstraint * alpha_background_shapeConstraint * SR_model(obs_x_SR) ] = 26.5578

functions
--------
RooProduct::L_x_background_SR_overallSyst_x_HistSyst[ 1 * background_SR_overallSyst_x_HistSyst ] = 219.314
RooProduct::L_x_signal_SR_overallSyst_x_HistSyst[ 1 * signal_SR_overallSyst_x_HistSyst ] = 46.264
PiecewiseInterpolation::background_SR_Hist_alpha[ ] = 219.314
RooHistFunc::background_SR_Hist_alpha_0high[ depList=(obs_x_SR) ] = 271.717
RooHistFunc::background_SR_Hist_alpha_0low[ depList=(obs_x_SR) ] = 189.227
RooHistFunc::background_SR_Hist_alphanominal[ depList=(obs_x_SR) ] = 219.314
RooStats::HistFactory::FlexibleInterpVar::background_SR_epsilon[ paramList=(alpha_background_uncertainty) ] = 1
RooProduct::background_SR_overallSyst_x_HistSyst[ background_SR_Hist_alpha * background_SR_epsilon ] = 219.314
PiecewiseInterpolation::signal_SR_Hist_alpha[ ] = 46.264
RooHistFunc::signal_SR_Hist_alpha_0high[ depList=(obs_x_SR) ] = 49.4179
RooHistFunc::signal_SR_Hist_alpha_0low[ depList=(obs_x_SR) ] = 43.9572
RooHistFunc::signal_SR_Hist_alphanominal[ depList=(obs_x_SR) ] = 46.264
RooStats::HistFactory::FlexibleInterpVar::signal_SR_epsilon[ paramList=(alpha_signal_norm_uncertainty) ] = 1
RooProduct::signal_SR_overallNorm_x_sigma_epsilon[ SigXsecOverSM * signal_SR_epsilon ] = 1
RooProduct::signal_SR_overallSyst_x_HistSyst[ signal_SR_Hist_alpha * signal_SR_overallNorm_x_sigma_epsilon ] = 46.264

datasets
--------
RooDataSet::asimovData(obs_x_SR)
RooDataSet::obsData(obs_x_SR)

embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_SR_Hist_alphanominalDHist(obs_x_SR)
RooDataHist::signal_SR_Hist_alpha_0lowDHist(obs_x_SR)
RooDataHist::signal_SR_Hist_alpha_0highDHist(obs_x_SR)
RooDataHist::background_SR_Hist_alphanominalDHist(obs_x_SR)
RooDataHist::background_SR_Hist_alpha_0lowDHist(obs_x_SR)
RooDataHist::background_SR_Hist_alpha_0highDHist(obs_x_SR)

named sets
----------
ModelConfig_GlobalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nom_alpha_background_uncertainty,nom_alpha_background_shape)
ModelConfig_Observables:(obs_x_SR)
coefList:(binWidth_obs_x_SR_0,binWidth_obs_x_SR_1)
constraintTerms:(lumiConstraint,alpha_signal_norm_uncertaintyConstraint,alpha_signal_shapeConstraint,alpha_background_uncertaintyConstraint,alpha_background_shapeConstraint)
globalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nom_alpha_background_uncertainty,nom_alpha_background_shape)
likelihoodTerms:(SR_model)
obsAndWeight:(weightVar,obs_x_SR)
observables:(obs_x_SR)
observablesSet:(obs_x_SR)
shapeList:(L_x_signal_SR_overallSyst_x_HistSyst,L_x_background_SR_overallSyst_x_HistSyst)

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

Setting Parameter(s) of Interest as: SigXsecOverSM
In [21]:
workspace.SetName('HistFactoryWorkspace')
workspace.writeToFile("output/HistFactoryWorkspace.root")
Out[21]:
False

Adding a Control Region

Here the CR is the region around the Z peak. 80<mmc<100

In [22]:
CR = [variables for variables in VBF[1:] if (float(variables[1]) < 100 and float(variables[1]) > 80) ]
CR_Loose = [variables for variables in VBF_Loose[1:] if (float(variables[1]) < 100 and float(variables[1]) > 80) ]
CR_Tight = [variables for variables in VBF_Tight[1:] if (float(variables[1]) < 100 and float(variables[1]) > 80) ]
In [23]:
h_params = [1,80,100]
data_cr, sig_cr, bkg_cr, ns, nb = get_var(CR,"CR",h_params)
data_up_cr, sig_up_cr, bkg_up_cr, nsup, nbup = get_var(CR_Loose,"CR_Loose",h_params)
data_down_cr, sig_down_cr, bkg_down_cr, nsdown, nbdown = get_var(CR_Tight,"CR_Tight",h_params)

print 'n signal {0}(+{1}%,{2}%)'.format(ns,nsup/ns*100-100,nsdown/ns*100-100)
print 'n background {0}(+{1}%,{2}%)'.format(nb,nbup/nb*100-100,nbdown/nb*100-100)

n signal 49.7659558473(+10.2652888066%,-7.27659576988%)
n background 1348.93556723(+23.4177300037%,-18.9804492244%)
In [24]:
chan = ROOT.RooStats.HistFactory.Channel( "CR" )
chan.SetStatErrorConfig(0.05, "Poisson")
chan.SetData( data_cr )
signal = ROOT.RooStats.HistFactory.Sample( "signal" )
signal.SetNormalizeByTheory( False )
signal.SetHisto( sig_cr )
uncertainty_up   = nsup/ns
uncertainty_down = nsdown/ns
signal.AddOverallSys( "signal_norm_uncertainty",  uncertainty_down, uncertainty_up )
signal_shape = ROOT.RooStats.HistFactory.HistoSys("background_shape")
signal_shape.SetHistoHigh( sig_up_cr )
signal_shape.SetHistoLow( sig_down_cr )
background.AddHistoSys( background_shape )
chan.AddSample( signal )

background = ROOT.RooStats.HistFactory.Sample( "background" )
background.SetNormalizeByTheory( False )
background.SetHisto( bkg_cr )
uncertainty_up   = nbup/nb
uncertainty_down = nbdown/nb
background.AddOverallSys( "background_uncertainty",  uncertainty_down, uncertainty_up )
background_shape = ROOT.RooStats.HistFactory.HistoSys("background_shape")
background_shape.SetHistoHigh( bkg_up_cr )
background_shape.SetHistoLow( bkg_down_cr )
background.AddHistoSys( background_shape )
chan.AddSample( background )

meas.AddChannel(chan)
In [25]:
hist2workspace = ROOT.RooStats.HistFactory.HistoToWorkspaceFactoryFast(meas)
workspace = hist2workspace.MakeCombinedModel(meas)


-------------------
Starting to process SR channel with 1 observables
lumi str = [1,0,10]
lumi Error str = nominalLumi[1,0,1.2],0.02
making normFactor: SigXsecOverSM
Gaussian::alpha_signal_shapeConstraint(alpha_signal_shape,nom_alpha_signal_shape[0.,-10,10],1.)
Gaussian::alpha_background_shapeConstraint(alpha_background_shape,nom_alpha_background_shape[0.,-10,10],1.)
-----------------------------------------
import model into workspace
RooDataSet::AsimovData[obs_x_SR,weight:binWeightAsimov] = 5 entries (2241.84 weighted)

RooWorkspace(SR) SR workspace contents

variables
---------
(Lumi,SigXsecOverSM,alpha_background_shape,alpha_background_uncertainty,alpha_signal_norm_uncertainty,alpha_signal_shape,binWidth_obs_x_SR_0,binWidth_obs_x_SR_1,nom_alpha_background_shape,nom_alpha_background_uncertainty,nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nominalLumi,obs_x_SR,weightVar)

p.d.f.s
-------
RooRealSumPdf::SR_model[ binWidth_obs_x_SR_0 * L_x_signal_SR_overallSyst_x_HistSyst + binWidth_obs_x_SR_1 * L_x_background_SR_overallSyst_x_HistSyst ] = 26.5578/2241.84
RooGaussian::alpha_background_shapeConstraint[ x=alpha_background_shape mean=nom_alpha_background_shape sigma=1 ] = 1
RooGaussian::alpha_background_uncertaintyConstraint[ x=alpha_background_uncertainty mean=nom_alpha_background_uncertainty sigma=1 ] = 1
RooGaussian::alpha_signal_norm_uncertaintyConstraint[ x=alpha_signal_norm_uncertainty mean=nom_alpha_signal_norm_uncertainty sigma=1 ] = 1
RooGaussian::alpha_signal_shapeConstraint[ x=alpha_signal_shape mean=nom_alpha_signal_shape sigma=1 ] = 1
RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.02 ] = 1
RooProdPdf::model_SR[ lumiConstraint * alpha_signal_norm_uncertaintyConstraint * alpha_signal_shapeConstraint * alpha_background_uncertaintyConstraint * alpha_background_shapeConstraint * SR_model(obs_x_SR) ] = 26.5578

functions
--------
RooProduct::L_x_background_SR_overallSyst_x_HistSyst[ 1 * background_SR_overallSyst_x_HistSyst ] = 219.314
RooProduct::L_x_signal_SR_overallSyst_x_HistSyst[ 1 * signal_SR_overallSyst_x_HistSyst ] = 46.264
PiecewiseInterpolation::background_SR_Hist_alpha[ ] = 219.314
RooHistFunc::background_SR_Hist_alpha_0high[ depList=(obs_x_SR) ] = 271.717
RooHistFunc::background_SR_Hist_alpha_0low[ depList=(obs_x_SR) ] = 189.227
RooHistFunc::background_SR_Hist_alphanominal[ depList=(obs_x_SR) ] = 219.314
RooStats::HistFactory::FlexibleInterpVar::background_SR_epsilon[ paramList=(alpha_background_uncertainty) ] = 1
RooProduct::background_SR_overallSyst_x_HistSyst[ background_SR_Hist_alpha * background_SR_epsilon ] = 219.314
PiecewiseInterpolation::signal_SR_Hist_alpha[ ] = 46.264
RooHistFunc::signal_SR_Hist_alpha_0high[ depList=(obs_x_SR) ] = 49.4179
RooHistFunc::signal_SR_Hist_alpha_0low[ depList=(obs_x_SR) ] = 43.9572
RooHistFunc::signal_SR_Hist_alphanominal[ depList=(obs_x_SR) ] = 46.264
RooStats::HistFactory::FlexibleInterpVar::signal_SR_epsilon[ paramList=(alpha_signal_norm_uncertainty) ] = 1
RooProduct::signal_SR_overallNorm_x_sigma_epsilon[ SigXsecOverSM * signal_SR_epsilon ] = 1
RooProduct::signal_SR_overallSyst_x_HistSyst[ signal_SR_Hist_alpha * signal_SR_overallNorm_x_sigma_epsilon ] = 46.264

datasets
--------
RooDataSet::asimovData(obs_x_SR)
RooDataSet::obsData(obs_x_SR)

embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_SR_Hist_alphanominalDHist(obs_x_SR)
RooDataHist::signal_SR_Hist_alpha_0lowDHist(obs_x_SR)
RooDataHist::signal_SR_Hist_alpha_0highDHist(obs_x_SR)
RooDataHist::background_SR_Hist_alphanominalDHist(obs_x_SR)
RooDataHist::background_SR_Hist_alpha_0lowDHist(obs_x_SR)
RooDataHist::background_SR_Hist_alpha_0highDHist(obs_x_SR)

named sets
----------
ModelConfig_GlobalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nom_alpha_background_uncertainty,nom_alpha_background_shape)
ModelConfig_Observables:(obs_x_SR)
coefList:(binWidth_obs_x_SR_0,binWidth_obs_x_SR_1)
constraintTerms:(lumiConstraint,alpha_signal_norm_uncertaintyConstraint,alpha_signal_shapeConstraint,alpha_background_uncertaintyConstraint,alpha_background_shapeConstraint)
globalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nom_alpha_background_uncertainty,nom_alpha_background_shape)
likelihoodTerms:(SR_model)
obsAndWeight:(weightVar,obs_x_SR)
observables:(obs_x_SR)
observablesSet:(obs_x_SR)
shapeList:(L_x_signal_SR_overallSyst_x_HistSyst,L_x_background_SR_overallSyst_x_HistSyst)

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

Setting Parameter(s) of Interest as: SigXsecOverSM 


-------------------
Starting to process CR channel with 1 observables
lumi str = [1,0,10]
lumi Error str = nominalLumi[1,0,1.2],0.02
signal_CR has no variation histograms 
processing hist CR signal
Gaussian::alpha_background_shapeConstraint(alpha_background_shape,nom_alpha_background_shape[0.,-10,10],1.)
-----------------------------------------
import model into workspace
RooDataSet::AsimovData[obs_x_CR,weight:binWeightAsimov] = 1 entries (1398.7 weighted)

RooWorkspace(CR) CR workspace contents

variables
---------
(Lumi,alpha_background_shape,alpha_background_uncertainty,alpha_signal_norm_uncertainty,binWidth_obs_x_CR_0,binWidth_obs_x_CR_1,nom_alpha_background_shape,nom_alpha_background_uncertainty,nom_alpha_signal_norm_uncertainty,nominalLumi,obs_x_CR,weightVar)

p.d.f.s
-------
RooRealSumPdf::CR_model[ binWidth_obs_x_CR_0 * L_x_signal_CR_overallSyst_x_Exp + binWidth_obs_x_CR_1 * L_x_background_CR_overallSyst_x_HistSyst ] = 69.9351/1398.7
RooGaussian::alpha_background_shapeConstraint[ x=alpha_background_shape mean=nom_alpha_background_shape sigma=1 ] = 1
RooGaussian::alpha_background_uncertaintyConstraint[ x=alpha_background_uncertainty mean=nom_alpha_background_uncertainty sigma=1 ] = 1
RooGaussian::alpha_signal_norm_uncertaintyConstraint[ x=alpha_signal_norm_uncertainty mean=nom_alpha_signal_norm_uncertainty sigma=1 ] = 1
RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.02 ] = 1
RooProdPdf::model_CR[ lumiConstraint * alpha_signal_norm_uncertaintyConstraint * alpha_background_uncertaintyConstraint * alpha_background_shapeConstraint * CR_model(obs_x_CR) ] = 69.9351

functions
--------
RooProduct::L_x_background_CR_overallSyst_x_HistSyst[ 1 * background_CR_overallSyst_x_HistSyst ] = 1348.94
RooProduct::L_x_signal_CR_overallSyst_x_Exp[ 1 * signal_CR_overallSyst_x_Exp ] = 49.766
PiecewiseInterpolation::background_CR_Hist_alpha[ ] = 1348.94
RooHistFunc::background_CR_Hist_alpha_0high[ depList=(obs_x_CR) ] = 1664.83
RooHistFunc::background_CR_Hist_alpha_0low[ depList=(obs_x_CR) ] = 1092.9
RooHistFunc::background_CR_Hist_alphanominal[ depList=(obs_x_CR) ] = 1348.94
RooStats::HistFactory::FlexibleInterpVar::background_CR_epsilon[ paramList=(alpha_background_uncertainty) ] = 1
RooProduct::background_CR_overallSyst_x_HistSyst[ background_CR_Hist_alpha * background_CR_epsilon ] = 1348.94
RooStats::HistFactory::FlexibleInterpVar::signal_CR_epsilon[ paramList=(alpha_signal_norm_uncertainty) ] = 1
RooHistFunc::signal_CR_nominal[ depList=(obs_x_CR) ] = 49.766
RooProduct::signal_CR_overallSyst_x_Exp[ signal_CR_nominal * signal_CR_epsilon ] = 49.766

datasets
--------
RooDataSet::asimovData(obs_x_CR)
RooDataSet::obsData(obs_x_CR)

embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_CRnominalDHist(obs_x_CR)
RooDataHist::background_CR_Hist_alphanominalDHist(obs_x_CR)
RooDataHist::background_CR_Hist_alpha_0lowDHist(obs_x_CR)
RooDataHist::background_CR_Hist_alpha_0highDHist(obs_x_CR)

named sets
----------
ModelConfig_GlobalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_background_uncertainty,nom_alpha_background_shape)
ModelConfig_Observables:(obs_x_CR)
coefList:(binWidth_obs_x_CR_0,binWidth_obs_x_CR_1)
constraintTerms:(lumiConstraint,alpha_signal_norm_uncertaintyConstraint,alpha_background_uncertaintyConstraint,alpha_background_shapeConstraint)
globalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_background_uncertainty,nom_alpha_background_shape)
likelihoodTerms:(CR_model)
obsAndWeight:(weightVar,obs_x_CR)
observables:(obs_x_CR)
observablesSet:(obs_x_CR)
shapeList:(L_x_signal_CR_overallSyst_x_Exp,L_x_background_CR_overallSyst_x_HistSyst)

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

Setting Parameter(s) of Interest as: SigXsecOverSM 
WARNING: Can't find parameter of interest: SigXsecOverSM in Workspace. Not setting in ModelConfig.
full list of observables:
RooArgList:: = (obs_x_SR,obs_x_CR)


------------------
 Entering combination
-----------------------------------------
create toy data for SR,CR
RooDataSet::AsimovData0[obs_x_SR,channelCat,weight:binWeightAsimov] = 5 entries (2241.84 weighted)
RooDataSet::AsimovData1[obs_x_CR,channelCat,weight:binWeightAsimov] = 1 entries (1398.7 weighted)
Merging data for channel SR
Merging data for channel CR

RooWorkspace(combined) combined contents

variables
---------
(channelCat,nom_alpha_background_shape,nom_alpha_background_uncertainty,nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,obs_x_CR,obs_x_SR,weightVar)

datasets
--------
RooDataSet::asimovData(obs_x_SR,obs_x_CR,weightVar,channelCat)
RooDataSet::obsData(channelCat,obs_x_SR,obs_x_CR)

named sets
----------
ModelConfig_GlobalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nom_alpha_background_uncertainty,nom_alpha_background_shape)
ModelConfig_Observables:(obs_x_SR,obs_x_CR,weightVar,channelCat)
globalObservables:(nom_alpha_signal_norm_uncertainty,nom_alpha_signal_shape,nom_alpha_background_uncertainty,nom_alpha_background_shape)
observables:(obs_x_SR,obs_x_CR,weightVar,channelCat)



----------------
 Importing combined model
setting Lumi constant
Setting Parameter(s) of Interest as: SigXsecOverSM
In [26]:
workspace.SetName('HistFactoryWorkspace')
workspace.writeToFile("output/HistFactoryCombinedWorkspace.root")
Out[26]:
False

A quick test!

In [27]:
infile = ROOT.TFile.Open("output/HistFactoryCombinedWorkspace.root")
w = infile.Get("HistFactoryWorkspace")
mc = w.obj("ModelConfig")
data = w.data("obsData")
x = w.var("obs_x_SR")
In [28]:
pl = ROOT.RooStats.ProfileLikelihoodCalculator(data,mc)
pl.SetConfidenceLevel(0.99)
In [29]:
interval = pl.GetInterval()
In [30]:
plot = ROOT.RooStats.LikelihoodIntervalPlot(interval)
plot.SetNPoints(50)
plot.SetMaximum(5)
c = ROOT.TCanvas()
plot.Draw()
c.Draw()
In [ ]: