# Set up the path to SModelS installation folder
import sys; sys.path.append("."); import smodels_paths
from imp import reload
from smodels.tools import runtime
from smodels import particlesLoader
from smodels.theory import slhaDecomposer,lheDecomposer
from smodels.tools.physicsUnits import fb, GeV, TeV
from smodels.theory.theoryPrediction import theoryPredictionsFor
from smodels.experiment.databaseObj import Database
from smodels.tools import coverage
from smodels.tools.smodelsLogging import setLogLevel
setLogLevel("info")
#Define your model (list of rEven and rOdd particles)
runtime.modelFile = 'smodels.share.models.mssm'
reload(particlesLoader) #Make sure all the model particles are up-to-date
# Path to input file (either a SLHA or LHE file)
slhafile = 'inputFiles/slha/lightEWinos.slha'
# Set main options for decomposition
sigmacut = 0.01 * fb
mingap = 5. * GeV
maxcond = 0.2
# Decompose model (use slhaDecomposer for SLHA input or lheDecomposer for LHE input)
toplist = slhaDecomposer.decompose(slhafile, sigmacut, doCompress=True, doInvisible=True, minmassgap=mingap)
# Access basic information from decomposition, using the topology list and topology objects:
print( "\n Decomposition Results: " )
print( "\t Total number of topologies: %i " %len(toplist) )
nel = sum([len(top.elementList) for top in toplist])
print( "\t Total number of elements = %i " %nel )
# Set the path to the database
database = Database("official")
# Load the experimental results to be used.
# In this case, all results are employed.
listOfExpRes = database.getExpResults()
# Print basic information about the results loaded.
# Count the number of loaded UL and EM experimental results:
nUL, nEM = 0, 0
for exp in listOfExpRes:
expType = exp.getValuesFor('dataType')[0]
if expType == 'upperLimit':
nUL += 1
elif expType == 'efficiencyMap':
nEM += 1
print("\n Loaded Database with %i UL results and %i EM results " %(nUL,nEM))
# Compute the theory predictions for each experimental result and print them:
print("\n Theory Predictions and Constraints:")
rmax = 0.
bestResult = None
for expResult in listOfExpRes:
predictions = theoryPredictionsFor(expResult, toplist)
if not predictions: continue # Skip if there are no constraints from this result
print('\n %s (%i TeV)' %(expResult.globalInfo.id,expResult.globalInfo.sqrts.asNumber(TeV)))
for theoryPrediction in predictions:
dataset = theoryPrediction.dataset
datasetID = theoryPrediction.dataId()
mass = theoryPrediction.mass
txnames = [str(txname) for txname in theoryPrediction.txnames]
PIDs = theoryPrediction.PIDs
print( "------------------------" )
print( "TxNames = ",txnames )
print( "Theory Prediction = ",theoryPrediction.xsection.value ) #Signal cross section
# Get the corresponding upper limit:
print( "UL for theory prediction = ",theoryPrediction.upperLimit )
# Compute the r-value
r = theoryPrediction.getRValue()
print( "r = ",r )
#Compute likelihhod and chi^2 for EM-type results:
if theoryPrediction.dataType() == 'efficiencyMap':
theoryPrediction.computeStatistics()
print( 'Chi2, likelihood=', theoryPrediction.chi2, theoryPrediction.likelihood )
# Check condition violation
exceedsMaxCond = False
CondViolation = theoryPrediction.getmaxCondition()
if CondViolation == 'N/A' or CondViolation == None:
print( "no condition violation" )
elif CondViolation <= maxcond:
print( "Condition violation = ", CondViolation, " (OK)" )
else:
print( "Condition violation ", CondViolation, " exceeds chosen bound of ", maxcond )
exceedsMaxCond=True
if r > rmax and exceedsMaxCond == False:
rmax = r
bestResult = expResult.globalInfo.id
# Print the most constraining experimental result
print( "\nThe largest r-value (theory/upper limit ratio) is ",rmax )
if rmax > 1.:
print( "(The input model is likely excluded by %s)" %bestResult )
else:
print( "(The input model is not excluded by the simplified model results)" )
#Find out missing topologies for sqrts=8*TeV:
uncovered = coverage.Uncovered(toplist,sqrts=13.*TeV)
#Print uncovered cross-sections:
print( "\nTotal missing topology cross section (fb): %10.3E\n" %(uncovered.getMissingXsec()) )
print( "Total cross section where we are outside the mass grid (fb): %10.3E\n" %(uncovered.getOutOfGridXsec()) )
print( "Total cross section in long cascade decays (fb): %10.3E\n" %(uncovered.getLongCascadeXsec()) )
print( "Total cross section in decays with asymmetric branches (fb): %10.3E\n" %(uncovered.getAsymmetricXsec()) )
#Print some of the missing topologies:
print( 'Missing topologies (up to 3):' )
for topo in uncovered.missingTopos.topos[:3]:
print( 'Topology:',topo.topo )
print( 'Contributing elements (up to 2):' )
for el in topo.contributingElements[:2]:
print( el,'cross-section (fb):', el.missingX )
#Print elements with long cascade decay:
print( '\nElements outside the grid (up to 2):' )
for topo in uncovered.outsideGrid.topos[:2]:
print( 'Topology:',topo.topo )
print( 'Contributing elements (up to 4):' )
for el in topo.contributingElements[:4]:
print( el,'cross-section (fb):', el.missingX )
print( '\tmass:',el.getMasses() )