#!/usr/bin/env python3
"""
.. module:: coverage
:synopsis: Definitions of classes used to find, group and format missing topologies
.. moduleauthor:: Ursula Laa <ursula.laa@lpsc.in2p3.fr>
.. moduleauthor:: Andre Lessa <lessa.a.p@gmail.com>
.. moduleauthor:: Alicia Wongel <alicia.wongel@gmail.com>
"""
from smodels.tools.physicsUnits import fb
from smodels.tools.reweighting import reweightFactorFor
from smodels.theory.branch import Branch
from smodels.theory.particle import MultiParticle, ParticleList, Particle
from smodels.share.models.SMparticles import e,mu,ta,taC,eC,muC,W,WC,t,tC,q,c,g,pion,nu
from smodels.theory.auxiliaryFunctions import index_bisect
from smodels.theory.exceptions import SModelSTheoryError as SModelSError
#Default definitions for the uncovered topology categories/groups:
##Element filters for each group:
##(it should be a function which takes an Element object as input
##and returns True if the element belongs to the group and False otherwise)
filtersDefault = {'missing (prompt)' : lambda el: not ('prompt' in el.coveredBy),
'missing (displaced)' : lambda el: not ('displaced' in el.coveredBy),
# 'missing (long cascade)' : lambda el: (not el.coveredBy) and el._getLength() > 3,
'missing (all)' : lambda el: (not el.coveredBy),
'outsideGrid (all)' : lambda el: (el.coveredBy and not el.testedBy)}
##Description for each group (optional and only used for printing):
##(if not defined, the group label will be used instead)
descriptionDefault = {'missing (prompt)' : 'missing topologies with prompt decays',
'missing (displaced)' : 'missing topologies with displaced decays',
# 'missing (long cascade)' : 'missing topologies with long cascade decays',
'missing (all)' : 'missing topologies',
'outsideGrid (all)' : 'topologies outside the grid'}
##Default final BSM states for grouping topologies:
#Used to construct BSM final states:
MET = Particle(label='MET', Z2parity = -1, eCharge = 0, colordim = 1)
HSCPp = Particle(label='HSCP+', Z2parity = -1, eCharge = +1, colordim = 1)
HSCPm = Particle(label='HSCP-', Z2parity = -1, eCharge = -1, colordim = 1)
HSCP = MultiParticle(label='HSCP', particles = [HSCPp,HSCPm])
RHadronG = Particle(label='RHadronG', Z2parity = -1, eCharge = 0, colordim = 8)
RHadronU = Particle(label='RHadronU', Z2parity = -1, eCharge = 2./3., colordim = 3)
RHadronD = Particle(label='RHadronD', Z2parity = -1, eCharge = -1./3., colordim = 3)
RHadronQ = MultiParticle(label='RHadronQ', particles = [RHadronU,RHadronU.chargeConjugate(),
RHadronD,RHadronD.chargeConjugate()])
bsmDefault = [MET,HSCP,RHadronG,RHadronQ]
##Weight factors for each group:
##(it should be a function which takes an Element object as input
##and returns the reweighting factor to be applied to the element weight. It is relevant if only
##the fraction of the weight going into prompt or displaced decays is required)
factorsDefault = {}
for key in filtersDefault:
if 'prompt' in key.lower():
factorsDefault[key] = lambda el: reweightFactorFor(el,'prompt')
elif 'displaced' in key.lower():
factorsDefault[key] = lambda el: reweightFactorFor(el,'displaced')
else:
#If not specified assumed all fractions
#(note that we can not include any long-lived fraction since this is already included in
#the topologies where the meta-stable particle appears as a final state,
#so the total is = (fraction of all decays being prompt)
#+ (fraction of at least one displaced decay and no long-lived decays)
factorsDefault[key] = lambda el: reweightFactorFor(el,'prompt') + reweightFactorFor(el,'displaced')
[docs]class Uncovered(object):
"""
Wrapper object for defining and holding a list of coverage groups (UncoveredGroup objects).
The class builds a series of UncoveredGroup objects and stores them.
"""
def __init__(self,topoList, sqrts=None, sigmacut=0*fb,
groupFilters = filtersDefault,
groupFactors = factorsDefault,
groupdDescriptions = descriptionDefault,
smFinalStates=None,bsmFinalSates=None):
"""
Inititalize the object.
:param topoList: TopologyList object used to select elements from.
:param sqrts: Value (with units) for the center of mass energy used to compute the missing cross sections.
If not specified the largest value available will be used.
:param sigmacut: Minimum cross-section/weight value (after applying the reweight factor)
for an element to be included. The value should in fb (unitless)
:param groupFilters: Dictionary containing the groups' labels and the method for selecting
elements.
:param groupFactors: Dictionary containing the groups' labels and the method for reweighting
cross sections.
:param groupdDescriptions: Dictionary containing the groups' labels and strings describing the group
(used for printout)
:param smFinalStates: List of (inclusive) Particle or MultiParticle objects used for grouping Z2-even
particles when creating GeneralElements.
:param bsmFinalSates: List of (inclusive) Particle or MultiParticle objects used for grouping Z2-odd
particles when creating GeneralElements.
"""
#Sanity checks:
if not isinstance(groupFilters,dict):
raise SModelSError("groupFilters input should be a Dictionary and not %s" %type(groupFilters))
if not isinstance(groupFactors,dict):
raise SModelSError("groupFactors input should be a Dictionary and not %s" %type(groupFactors))
if sorted(groupFilters.keys()) != sorted(groupFactors.keys()):
raise SModelSError("Keys in groupFilters and groupFactors do not match")
if any(not hasattr(gFilter,'__call__') for gFilter in groupFilters.values()):
raise SModelSError("Group filters must be callable methods")
if smFinalStates is None:
#Define multiparticles for conveniently grouping SM final states
taList = MultiParticle('ta', [ta,taC])
lList = MultiParticle('l', [e,mu,eC,muC])
WList = MultiParticle('W', [W,WC])
tList = MultiParticle('t', [t,tC])
jetList = MultiParticle('jet', [q,c,g,pion])
nuList = nu
smFinalStates = [WList, lList, tList, taList, nuList, jetList]
if bsmFinalSates is None:
#Define inclusive BSM states to group/label the last BSM states:
bsmFinalStates = bsmDefault
if sqrts is None:
sqrts = max([xsec.info.sqrts for xsec in topoList.getTotalWeight()])
else:
sqrts = sqrts
#Store the relevant element cross-sections to improve performance:
for el in topoList.getElements():
xsec = el.weight.getXsecsFor(sqrts)
if xsec:
el._totalXsec = xsec[0].value.asNumber(fb)
else:
el._totalXsec = 0.
self.groups = []
#Create each uncovered group and get the topologies from topoList
for gLabel,gFilter in groupFilters.items():
#Initialize the uncovered topology list:
uncoveredTopos = UncoveredGroup(label=gLabel,elementFilter=gFilter,
reweightFactor = groupFactors[gLabel],
smFinalStates=smFinalStates,
bsmFinalStates=bsmFinalStates,
sqrts=sqrts, sigmacut=sigmacut.asNumber(fb))
if groupdDescriptions and gLabel in groupdDescriptions:
uncoveredTopos.description = groupdDescriptions[gLabel]
else:
uncoveredTopos.description = gLabel
#Fill the list with the elements in topoList:
uncoveredTopos.getToposFrom(topoList)
self.groups.append(uncoveredTopos)
[docs] def getGroup(self,label):
"""
Returns the group with the required label. If not found, returns None.
:param label: String corresponding to the specific group label
:return: UncoveredGroup object which matches the label
"""
for group in self.groups:
if group.label == label:
return group
return None
[docs]class UncoveredGroup(object):
"""
Holds information about a single coverage group: criteria for selecting and grouping elements,
function for reweighting cross sections, etc.
"""
def __init__(self, label, elementFilter, reweightFactor,
smFinalStates, bsmFinalStates, sqrts, sigmacut=0.):
"""
:param label: Group label
:param elementFilter: Function which takes an element as argument and returns True (False) if
the element should (not) be selected.
:param reweightFactor: Function which takes an element as argument and returns the reweighting
factor to be applied to the element weight.
:param smFinalStates: List of Particle/MultiParticle objects used to group Z2-even particles appearing
in the final state
:param bsmFinalStates: List of Particle/MultiParticle objects used to group Z2-odd particles appearing
in the final state
:param sqrts: Value (with units) for the center of mass energy used to compute the missing cross sections.
If not specified the largest value available will be used.
:param sigmacut: Minimum cross-section/weight value (after applying the reweight factor)
for an element to be included. The value should in fb (unitless)
"""
self.generalElements = []
self.smFinalStates = smFinalStates
self.bsmFinalStates = bsmFinalStates
self.sqrts = sqrts
self.sigmacut = sigmacut
self.label = label
self.elementFilter = elementFilter
self.reweightFactor = reweightFactor
def __str__(self):
return self.label
def __repr__(self):
return str(self)
[docs] def getToposFrom(self,topoList):
"""
Select the elements from topoList according to self.elementFilter
and build GeneralElements from the selected elements.
The GeneralElement weights corresponds to the missing cross-section
with double counting from compressed elements already accounted for.
"""
#First select all elements according to the filter (type of uncovered/missing topology):
elementList = [el for el in topoList.getElements() if self.elementFilter(el)]
#Get missing xsections including the reweight factor:
missingXandEls = [[self.getMissingX(el)*self.reweightFactor(el),el] for el in elementList]
#Only keep the ones elements sigmacut:
if self.sigmacut:
missingXandEls = [x for x in missingXandEls[:] if x[0] > self.sigmacut]
else:
missingXandEls = [x for x in missingXandEls[:] if x[0] > 0.]
#Sort according to largest missingX, smallest size and largest ID
missingXandEls = sorted(missingXandEls, key = lambda pt: [pt[0],-pt[1]._getLength(),-pt[1].elID], reverse=True)
#Split lists of elements and missingX:
missingXsecs = [pt[0] for pt in missingXandEls]
elementList = [pt[1] for pt in missingXandEls]
#Remove all elements which are related to each other in order to avoid double counting
#(keep always the first appearance in the list, so we always keep the ones with largest missing xsec)
elementListUnique = []
missingXsecsUnique = []
ancestors = set() #Keep track of all the ancestor ids of the elements in the unique list
for i,element in enumerate(elementList):
#Get ancestor object ids
#(safer than element.elID since some elements can share elID = 0 if they
#were never directly inserted in the TopologyList)
ancestorsIDs = set([id(element)] + [id(el) for el in element.getAncestors()])
#If the element has any common ancestor with any of the previous elements
#or if it is an ancestor of any of the previous elements
#skip it to avoid double counting
if ancestors.intersection(ancestorsIDs):
continue
elementListUnique.append(element)
missingXsecsUnique.append(missingXsecs[i])
ancestors = ancestors.union(ancestorsIDs)
#Now that we only have unique elements with their effective missing cross-sections
#we create General Elements out of them
for i,el in enumerate(elementListUnique):
missingX = missingXsecsUnique[i]
if not missingX:
continue
self.addToGeneralElements(el,missingX)
#Finally sort general elements by their missing cross-section:
self.generalElements = sorted(self.generalElements[:], key = lambda genEl: genEl.missingX, reverse=True)
[docs] def getMissingX(self,element):
"""
Calculate total missing cross section of an element, by recursively checking if its
mothers already appear in the list.
:param element: Element object
:returns: missing cross section without units (in fb)
"""
ancestorList = element.getAncestors()
alreadyChecked = [] # keep track of which elements have already been checked
#Get the (pre-loaded) total element weight in fb:
missingX = element._totalXsec
if not missingX:
return 0.
overlapXsec = 0.
for ancestor in ancestorList: #check all ancestors (ancestorList is sorted by generation)
#Skip entries which correspond to the element itself
if ancestor is element:
continue
#Make sure we do not subtract the same mother twice
if any(ancestor is el for el in alreadyChecked):
continue
alreadyChecked.append(ancestor)
#Check if ancestor passes the group filter (if it has been covered/tested or not):
if self.elementFilter(ancestor):
continue
#Subtract the weight of the ancestor and skip checking for all the older family tree of the ancestor
#(avoid subtracting the same weight twice from mother and grandmother for instance)
#(since the ancestorList is sorted by generation, the mother always
#appears before the grandmother in the list)
alreadyChecked += ancestor.getAncestors()
if not hasattr(ancestor,'_totalXsec'):
xsec = ancestor.weight.getXsecsFor(self.sqrts)
ancestor._totalXsec = xsec[0].value.asNumber(fb)
overlapXsec += ancestor._totalXsec
return missingX-overlapXsec
[docs] def getTotalXSec(self, sqrts=None):
"""
Calculate total missing topology cross section at sqrts. If no sqrts is given use self.sqrts
:ivar sqrts: sqrts
"""
xsec = 0.
if not sqrts: sqrts = self.sqrts
for genEl in self.generalElements:
xsec += genEl.missingX
return xsec
[docs] def addToGeneralElements(self, el, missingX):
"""
Adds an element to the list of missing topologies = general elements.
If the element contributes to a missing topology that is already in the list, add element and weight to topology.
:parameter el: element to be added
:parameter missingX: missing cross-section for the element (in fb)
"""
newGenEl = GeneralElement(el,missingX,self.smFinalStates,self.bsmFinalStates)
index = index_bisect(self.generalElements,newGenEl)
if index != len(self.generalElements) and self.generalElements[index] == newGenEl:
self.generalElements[index]._contributingElements.append(el)
self.generalElements[index].missingX += missingX
else:
self.generalElements.insert(index,newGenEl)
[docs]class GeneralElement(object):
"""
This class represents a simplified (general) element which does
only holds information about its even particles and decay type.
The even particles are replaced/grouped by the particles defined in smFinalStates.
"""
def __init__(self,el,missingX,smFinalStates,bsmFinalStates):
self.branches = [Branch() for _ in el.branches]
self.missingX = missingX
self.finalBSMstates = []
for ib,branch in enumerate(el.branches):
newParticles = []
for vertex in branch.evenParticles:
newVertex = vertex[:]
for ip,particle in enumerate(vertex):
for particleList in smFinalStates:
if particleList.contains(particle):
newVertex[ip] = particleList
newParticles.append(ParticleList(newVertex))
self.branches[ib].evenParticles = newParticles
finalBSM = branch.oddParticles[-1]
for bsmFS in bsmFinalStates:
if finalBSM == bsmFS:
finalBSM = bsmFS
break
self.branches[ib].finalBSMstate = finalBSM
self.branches[ib].setInfo()
self.sortBranches()
self.finalBSMstates = [br.finalBSMstate for br in self.branches]
self._contributingElements = [el]
self._outputDescription = str(self)
def __cmp__(self,other):
"""
Compares the element with other for any branch ordering.
The comparison is made based on branches.
OBS: The elements and the branches must be sorted!
:param other: element to be compared (GeneralElement object)
:return: -1 if self < other, 0 if self == other, +1, if self > other.
"""
if not isinstance(other,GeneralElement):
return -1
#Compare branches:
if self.branches != other.branches:
comp = (self.branches > other.branches)
if comp:
return 1
else:
return -1
#Compare decay type
if self.finalBSMstates != other.finalBSMstates:
comp = (self.finalBSMstates > other.finalBSMstates)
if comp:
return 1
else:
return -1
return 0
def __eq__(self,other):
return self.__cmp__(other)==0
def __lt__(self,other):
return self.__cmp__(other)<0
def __str__(self):
"""
Create the element bracket notation string, e.g. [[[jet]],[[jet]],
including the decay type.
:returns: string representation of the element (in bracket notation)
"""
elStr = "["+",".join([str(br) for br in self.branches])+"]"
elStr = elStr.replace(" ", "").replace("'", "")
name = elStr.replace('~','') + ' (%s)'%(','.join([str(p) for p in self.finalBSMstates]))
return name
def __repr__(self):
return self.__str__()
[docs] def sortBranches(self):
"""
Sort branches. The smallest branch is the first one.
If branches are equal, sort accoding to decayType.
"""
#Now sort branches
self.branches = sorted(self.branches, key = lambda br: (br,br.finalBSMstate))