Source code for theory.model

.. module:: model
   :synopsis: Create model from file

.. moduleauthor:: Alicia Wongel <>

import pyslha
from import logger
from import GeV
from smodels.theory import lheReader, crossSection
from smodels.theory.particle import Particle,MultiParticle,ParticleList
from smodels.theory.exceptions import SModelSTheoryError as SModelSError

[docs]class Model(object): """ An instance of this class holds all the relevant information from the input model. This class contains the input file, the particles of the model and their production cross-sections. """ def __init__(self, BSMparticles, SMparticles, label=None): """ Initializes the model :parameter BSMparticles: list with BSM particle objects :parameter SMparticles: list with SM particle objects :parameter label: Optional string to label the model """ self.inputFile = None #allBSMparticles stores all possible particles for a given model: self.allBSMparticles = [particle.copy() for particle in BSMparticles] #BSMparticles is used to store the model particles defined by the inputFile self.BSMparticles = self.allBSMparticles[:] #at initialization assume all particles will be used self.SMparticles = SMparticles[:] self.label = label #Check if for each PDG there is a unique particle object defined allPDGs = self.getValuesFor('pdg') for pdg in allPDGs: p = self.getParticlesWith(pdg=pdg) if len(p) > 1: raise SModelSError("PDG = %s has been defined for multiple particles (%s). Check your model definitions." %(pdg,p)) def __str__(self): if not self.label: return 'Model: %s' %self.inputFile else: return 'Model: %s' %self.label def __repr__(self): return str(self) def __eq__(self,other): """ Simple model comparison. Check if all the particles are the same (including their ordering) and if the model label is the same """ if self.label != other.label: return False if len(self.SMparticles) != len(other.SMparticles): return False if len(self.BSMparticles) != len(other.BSMparticles): return False if self.SMparticles != other.SMparticles: return False if self.BSMparticles != other.BSMparticles: return False return True
[docs] def getParticlesWith(self,**kwargs): """ Return the particle objects with the listed attributes. MultiParticle objects are added if any of its particles matches the listed attributes. In order to avoid double counting, MultiParticle objects are only included if they do not contain any of the particles already in the list. For instance, if MP.particles = [e-,e+] and [e-] already appears in the list, MP will not be added. :returns: List of particle objects """ particleList = [] for particle in self.BSMparticles + self.SMparticles: pList = [particle] if isinstance(particle,MultiParticle): pList += particle.particles for p in pList: if any(not hasattr(p,attr) for attr in kwargs.keys()): continue if any(getattr(p,attr) != value for attr,value in kwargs.items()): continue particleList.append(particle) #Remove repeated entries. If the list contains a Particle and a MultiParticle #which contains the Particle, remove the MultiParticle cleanList = [] #First include all Particle objects: for particle in particleList: if isinstance(particle,MultiParticle): continue if any(particle is ptc for ptc in cleanList): continue cleanList.append(particle) #Now only include MultiParticle objects, if they do not contain #any of the particles already in the list: for particle in particleList: if not isinstance(particle,MultiParticle): continue if any(particle is ptc for ptc in cleanList): continue if any((particle.contains(p) and not particle is p) for p in particleList): continue cleanList.append(particle) if not cleanList: logger.warning("Particle with attributes %s not found in %s" %(str(kwargs),self)) return cleanList
[docs] def getValuesFor(self,attributeStr): """ Returns a list with all the values for attribute appearing in the model. :param attributeStr: String for the desired attribute :returns: list of values """ valueList = [] for particle in self.BSMparticles+self.SMparticles: if not hasattr(particle,attributeStr): continue value = getattr(particle,attributeStr) if isinstance(particle,MultiParticle) and isinstance(value,list): for v in value: if not v in valueList: valueList.append(v) else: valueList.append(value) return valueList
[docs] def getModelDataFrom(self,inputFile): """ Reads the input file (LHE or SLHA) and extract the relevant information (masses, widths, BRs and cross-sections). If a http address is given, it will attempt to download the file. :param inputFile: input file (SLHA or LHE) :return: dictionary with masses, dictionary with decays and XSectionList object """ #Download input file, if requested if inputFile.startswith("http") or inputFile.startswith("ftp"):"Asked for remote slhafile %s. Fetching it." % inputFile) import requests import os.path r = requests.get(inputFile) if r.status_code != 200: logger.error("Could not retrieve remote file %d: %s" %(r.status_code, r.reason)) raise SModelSError() basename = os.path.basename(inputFile) f = open(basename, "w") f.write(r.text) f.close() inputFile = basename #Trick to suppress pyslha error messages: import sys storeErr = sys.stderr #Try to read file assuming it is an SLHA file: try: sys.stderr = None res = pyslha.readSLHAFile(inputFile) massDict = res.blocks['MASS'] #Make sure both PDG signs appear in massDict for pdg,mass in massDict.items(): if not -pdg in massDict: massDict[-pdg] = abs(mass) decaysDict = res.decays xsections = crossSection.getXsecFromSLHAFile(inputFile) #If fails assume it is an LHE file: except (IOError,AttributeError,KeyError): massDict,decaysDict = lheReader.getDictionariesFrom(inputFile) xsections = crossSection.getXsecFromLHEFile(inputFile)"Using LHE input. All unstable particles will be assumed to have prompt decays.")"Using LHE input. All particles not appearing in the events will be removed from the model.") finally: sys.stderr = storeErr return massDict,decaysDict,xsections
[docs] def getEvenOddList(self): """ Get the list of even and odd particles, according to the Z2 parity value defined for each particle. :return: list with PDGs of even particles, list with PDGs of odd particles """ allPDGs = list(set(self.getValuesFor('pdg'))) evenPDGs = [] oddPDGs = [] for pdg in allPDGs: if all(p.Z2parity == 1 for p in self.getParticlesWith(pdg=pdg)): evenPDGs.append(pdg) elif all(p.Z2parity == -1 for p in self.getParticlesWith(pdg=pdg)): oddPDGs.append(pdg) return evenPDGs,oddPDGs
[docs] def filterCrossSections(self): """ Remove cross-sections for even particles or particles which do not belong to the model. Valid cross-sections are stored in self.xsections :return: Number of cross-sections after filter. """ evenPDGs,oddPDGs = self.getEvenOddList() modelPDGs = [particle.pdg for particle in self.BSMparticles if isinstance(particle.pdg,int)] for xsec in self.xsections.xSections[:]: for pid in if not pid in modelPDGs: logger.debug("Cross-section for %s includes particles not belonging to model and will be ignored" %str( self.xsections.delete(xsec) if not pid in oddPDGs: logger.debug("Cross-section for %s includes even particles and will be ignored" %str( self.xsections.delete(xsec) return len(self.xsections.xSections)
[docs] def setMasses(self,massDict,roundMasses,minMass): """ Define particle masses using massDict. :param roundMasses: If set, it will round the masses to this number of digits (int) :param massDict: dictionary with PDGs as keys and masses as values. :param minMass: Minimal mass for BSM particles in the model. Any particle with mass below minMass will have its mass set to minMass. """ for particle in self.BSMparticles: if isinstance(particle,MultiParticle): continue if not hasattr(particle,'pdg') or not hasattr(particle,'Z2parity'): raise SModelSError("PDG and/or Z2-parity for particle %s has not been defined" %particle.label) pdg = particle.pdg if abs(pdg) in massDict.keys(): if roundMasses and int(roundMasses) > 0: particle.mass = round(abs(massDict[abs(pdg)]),int(roundMasses))*GeV else: particle.mass = abs(massDict[abs(pdg)])*GeV if minMass > 0*GeV: particle.mass = max(particle.mass,minMass) else: raise SModelSError("No mass found for %i in input file %s." %(particle.pdg,self.inputFile))
[docs] def setDecays(self,decaysDict,promptWidth,stableWidth,erasePrompt): allPDGs = list(set(self.getValuesFor('pdg'))) evenPDGs,oddPDGs = self.getEvenOddList() for particle in self.BSMparticles: if isinstance(particle,MultiParticle): continue if not hasattr(particle,'pdg') or not hasattr(particle,'Z2parity'): raise SModelSError("PDG and/or Z2-parity for particle %s has not been defined" %particle.label) pdg = particle.pdg particle.decays = [] if pdg in decaysDict: particleData = decaysDict[pdg] chargeConj = 1 elif -pdg in decaysDict: particleData = decaysDict[-pdg] chargeConj = -1 else: logger.error("Decay information for particle %i could not be found" %pdg) raise SModelSError() particle.totalwidth = abs(particleData.totalwidth)*GeV if particle.totalwidth < stableWidth: particle.totalwidth = 0.*GeV #Treat particle as stable logger.debug("Particle %s has width below the threshold and will be assumed as stable" %particle.pdg) continue if particle.totalwidth > promptWidth: particle.totalwidth = float('inf')*GeV #Treat particle as prompt logger.debug("Particle %s has width above the threshold and will be assumed as prompt." %particle.pdg) if erasePrompt and particle.Z2parity == -1: logger.debug("Erasing quantum numbers of (prompt) particle %s." %particle.pdg) for attr in erasePrompt: if hasattr ( particle, attr ): delattr(particle,attr) else: particle.decays.append(None) #Include possibility for particle being long-lived (non-prompt) for decay in particleData.decays: pids = decay.ids missingIDs = set(pids).difference(set(allPDGs)) if missingIDs:"Particle(s) %s is not defined within model. Decay %s will be ignored" %(missingIDs,decay)) continue oddPids = [pid for pid in decay.ids if abs(pid) in oddPDGs] evenPids = [pid for pid in decay.ids if abs(pid) in evenPDGs] if len(oddPids) != 1 or len(evenPids+oddPids) != len(decay.ids): logger.debug("Decay %i -> %s is not of the form Z2-odd -> Z2-odd + [Z2-even particles] and will be ignored" %(pdg,pids)) continue #Conjugated decays if needed #(if pid*chargeConj is not in model, assume the particle is its own anti-particle) decayIDs = [pid*chargeConj if pid*chargeConj in allPDGs else pid for pid in decay.ids] newDecay = pyslha.Decay(,nda=decay.nda,parentid=decay.parentid,ids=decayIDs) #Convert PDGs to particle objects: daughters = [] for pdg in newDecay.ids: daughter = self.getParticlesWith(pdg=pdg) if not daughter: raise SModelSError("Particle with PDG = %i was not found in model. Check the model definitions." %pdg) elif len(daughter) > 1: raise SModelSError("Multiple particles defined with PDG = %i. PDG ids must be unique." %pdg) else: daughter = daughter[0] daughters.append(daughter) oddParticles = [p for p in daughters if p.Z2parity == -1] evenParticles = ParticleList([p for p in daughters if p.Z2parity == 1]) newDecay.oddParticles = oddParticles newDecay.evenParticles = evenParticles particle.decays.append(newDecay) #Check if all unstable particles have decay channels defined: for p in self.BSMparticles: if p.totalwidth < stableWidth: continue ndecays = len([dec for dec in p.decays if dec is not None]) if ndecays == 0: if p.Z2parity == -1: logger.warning("No valid decay found for %s. It will be considered stable." %p) p.totalwidth = 0.*GeV
[docs] def updateParticles(self, inputFile, promptWidth = None, stableWidth = None, roundMasses = 1, erasePrompt=['spin','eCharge','colordim'], minMass=1.0*GeV): """ Update mass, total width and branches of allParticles particles from input SLHA or LHE file. :param inputFile: input file (SLHA or LHE) :param promptWidth: Maximum width for considering particles as decaying prompt. If None, it will be set 1e-11 GeV. :param stableWidth: Minimum width for considering particles as stable. If None, it will be set 1e-25 GeV. :param roundMasses: If set, it will round the masses to this number of digits (int) :param erasePrompt: If set, all particles with prompt decays (totalwidth > promptWidth) will have the corresponding properties (quantum numbers). So all particles with the same mass and Z2parity will be considered as equal when combining elements. :param minMass: Minimal mass for BSM particles in the model. Any particle with mass below minMass will have its mass set to minMass. """ self.inputFile = inputFile if promptWidth is None: promptWidth = 1e-11*GeV if stableWidth is None: stableWidth = 1e-25*GeV massDict,decaysDict,xsections = self.getModelDataFrom(self.inputFile) self.xsections = xsections #Restric BSM particles to the overlap between allBSMparticles and massDict: self.BSMparticles = [] for particle in self.allBSMparticles: if not particle.pdg in massDict and not -particle.pdg in massDict: continue self.BSMparticles.append(particle) if len(self.BSMparticles) == len(self.allBSMparticles):"Loaded %i BSM particles" %(len(self.BSMparticles))) else:"Loaded %i BSM particles (%i particles not found in %s)" %(len(self.BSMparticles),len(self.allBSMparticles)-len(self.BSMparticles), self.inputFile)) #Remove cross-sections for even particles: nXsecs = self.filterCrossSections() if nXsecs == 0: msg = "No cross-sections found in %s for the Z2 odd BSM particles. "%self.inputFile msg += "Check if the model is compatible with the input file." logger.error(msg) raise SModelSError(msg) #Set particle masses: self.setMasses(massDict,roundMasses,minMass) #Set particle decays self.setDecays(decaysDict,promptWidth,stableWidth,erasePrompt) #Reset particle equality from all particles: for p in Particle.getinstances(): p._comp = {p._id : 0} if isinstance(p,MultiParticle): for ptc in p.particles: p._comp[ptc._id] = 0 #Reset particle equality from all particle lists: for pL in ParticleList.getinstances(): pL._comp = {pL._id : 0}