import os,sys
# Set up the path to SModelS installation folder
import sys; sys.path.append("."); import smodels_paths
import numpy as np
import matplotlib.pyplot as plt

from smodels.matching.modelTester import getCombiner

Define input (SLHA) file and the parameters file

slhafile = os.path.expanduser("./inputFiles/slha/ewino_example.slha")
# Define parameters file with combineAnas = ATLAS-SUSY-2019-08,ATLAS-SUSY-2019-09:
parfile = os.path.expanduser("./parameters_comb.ini")

Define some basic parameters for plotting the likelihoods

expected = False # whether to plot the observed or expected likelihood
normalize = True # whether to normalize the likelihoods
muvals = np.linspace(0.,5.,200) # Signal strength values for which to evaluate the likelihoods

Run SModelS and get the analysis combination

combiner = getCombiner(slhafile, parfile)

Use the combination to evaluate the likelihoods

llhdDict = combiner.getLlhds(muvals,expected,normalize)

Compute L_SM, L_BSM, L_max and the UL on mu

muhat = combiner.muhat()
lmax = combiner.lmax()
lsm = combiner.lsm()
lbsm = combiner.likelihood(mu=1.0)
muULDict = {'combined' : combiner.getUpperLimitOnMu()}
for theoryPred in combiner.theoryPredictions:
    anaID = theoryPred.analysisId()
    muULDict[anaID] = theoryPred.getUpperLimitOnMu()

Plot the results

fig = plt.figure(figsize=(8,4))
ymin = 0.
ymax = 0.
for anaID,l in llhdDict.items():
    if anaID == 'combined':
        zorder = 100
        linestyle = '--'
    else:
        zorder = None
        linestyle = '-'
        
    p = plt.plot(muvals,l,label=anaID,linewidth=3,zorder=zorder,linestyle=linestyle)
    ymin = min(ymin,min(l))
    ymax = max(ymax,max(l))
    
    plt.vlines(muULDict[anaID],ymin=ymin,ymax=ymax,linestyle='--',linewidth=2,
           label=r'$\mu_{UL}$ (%s)' %anaID,color=p[0].get_color(),alpha=0.7)
    
    
plt.vlines(muhat,ymin=ymin,ymax=ymax,linestyle='--',linewidth=2,
           label=r'$\hat{\mu}$',color='black',alpha=0.7)


# plt.yscale('log')
# plt.ylim(1e-1,1e1)
plt.xlim(muvals.min(),muvals.max())
plt.xlabel(r'$\mu$')
if normalize:
    plt.ylabel('Normalized Likelihood')
else:
    plt.ylabel('Likelihood')
plt.legend(framealpha=1)
plt.title(r'$\hat{\mu} = $ %1.2f, $L_{max} =$ %1.2e, $L_{SM} =$ %1.2e, $L_{BSM} =$ %1.2e' %(muhat,lmax,lsm,lbsm))
plt.tight_layout()
plt.show()