pyslha.py

Fri, 21 Jan 2011 13:50:54 +0000

author
Andy Buckley <andy@insectnation.org>
date
Fri, 21 Jan 2011 13:50:54 +0000
changeset 109
0502c4942fbc
parent 106
513c294bc9dd
child 111
4f426f1f5923
permissions
-rw-r--r--

Fix units issue in slhaplot -- thanks to Martin Spinrath!

     1 #! /usr/bin/env python
     3 """\
     4 A simple but flexible parser of SUSY Les Houches Accord (SLHA) model and decay files.
     6 pyslha is a parser/writer module for particle physics SUSY Les Houches Accord
     7 (SLHA) supersymmetric spectrum/decay files, and a collection of scripts which
     8 use the interface, e.g. for conversion to and from the legacy ISAWIG format, or
     9 to plot the mass spectrum and decay chains.
    11 The current release supports SLHA version 1. Assistance with supporting version
    12 2 will be gladly accepted!
    14 The plotting script provides output in PDF via LaTeX and the TikZ graphics package,
    15 as LaTeX/TikZ source for direct embedding into documents, and in the format used by
    16 the Rivet make-plots system.
    18 TODOs:
    19  * Identify HERWIG decay matrix element to use in ISAWIG
    20  * Split writeSLHA into writeSLHA{Blocks,Decays}
    21  * Handle SLHA2
    22  * Handle RPV SUSY in ISAWIG
    23 """
    25 __author__ = "Andy Buckley <andy.buckley@cern.ch"
    26 __version__ = "1.0.5"
    29 def _autotype(var):
    30     """Automatically convert strings to numerical types if possible."""
    31     if type(var) is not str:
    32         return var
    33     if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
    34         return int(var)
    35     try:
    36         f = float(var)
    37         return f
    38     except ValueError:
    39         return var
    41 def _autostr(var):
    42     """Automatically numerical types to the right sort of string."""
    43     if type(var) is float:
    44         return "%e" % var
    45     return str(var)
    49 class Block(object):
    50     """
    51     Object representation of any BLOCK elements read from the SLHA file.  Blocks
    52     have a name, may have an associated Q value, and then a collection of data
    53     entries, stored as a recursive dictionary. Types in the dictionary are
    54     numeric (int or float) when a cast from the string in the file has been
    55     possible.
    56     """
    57     def __init__(self, name, q=None):
    58         self.name = name
    59         self.entries = {}
    60         self.q = _autotype(q)
    62     def add_entry(self, entry):
    63         #print entry
    64         nextparent = self.entries
    65         if len(entry) < 2:
    66             raise Exception("Block entries must be at least a 2-tuple")
    67         #print "in", entry
    68         entry = map(_autotype, entry)
    69         #print "out", entry
    70         for e in entry[:-2]:
    71             if e is not entry[-1]:
    72                 nextparent = nextparent.setdefault(e, {})
    73         nextparent[entry[-2]] = entry[-1]
    74         #print self.entries
    76     def __cmp__(self, other):
    77         return cmp(self.name, other.name)
    79     def __str__(self):
    80         s = self.name
    81         if self.q is not None:
    82             s += " (Q=%s)" % self.q
    83         s += "\n"
    84         s += str(self.entries)
    85         return s
    88 class Decay(object):
    89     """
    90     Object representing a decay entry on a particle decribed by the SLHA file.
    91     'Decay' objects are not a direct representation of a DECAY block in an SLHA
    92     file... that role, somewhat confusingly, is taken by the Particle class.
    94     Decay objects have three properties: a branching ratio, br, an nda number
    95     (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
    96     decay occurs. The PDG ID of the particle whose decay this represents may
    97     also be stored, but this is normally known via the Particle in which the
    98     decay is stored.
    99     """
   100     def __init__(self, br, nda, ids, parentid=None):
   101         self.parentid = parentid
   102         self.br = br
   103         self.nda = nda
   104         self.ids = ids
   105         assert(self.nda == len(self.ids))
   107     def __cmp__(self, other):
   108         return cmp(other.br, self.br)
   110     def __str__(self):
   111         return "%e %s" % (self.br, self.ids)
   114 class Particle(object):
   115     """
   116     Representation of a single, specific particle, decay block from an SLHA
   117     file.  These objects are not themselves called 'Decay', since that concept
   118     applies more naturally to the various decays found inside this
   119     object. Particle classes store the PDG ID (pid) of the particle being
   120     represented, and optionally the mass (mass) and total decay width
   121     (totalwidth) of that particle in the SLHA scenario. Masses may also be found
   122     via the MASS block, from which the Particle.mass property is filled, if at
   123     all. They also store a list of Decay objects (decays) which are probably the
   124     item of most interest.
   125     """
   126     def __init__(self, pid, totalwidth=None, mass=None):
   127         self.pid = pid
   128         self.totalwidth = totalwidth
   129         self.mass = mass
   130         self.decays = []
   132     def add_decay(self, br, nda, ids):
   133         self.decays.append(Decay(br, nda, ids))
   134         self.decays.sort()
   136     def __cmp__(self, other):
   137         if abs(self.pid) == abs(other.pid):
   138             return cmp(self.pid, other.pid)
   139         return cmp(abs(self.pid), abs(other.pid))
   141     def __str__(self):
   142         s = str(self.pid)
   143         if self.mass is not None:
   144             s += " : mass = %e GeV" % self.mass
   145         if self.totalwidth is not None:
   146             s += " : total width = %e GeV" % self.totalwidth
   147         for d in self.decays:
   148             if d.br > 0.0:
   149                 s += "\n  %s" % d
   150         return s
   153 def readSLHAFile(spcfilename, **kwargs):
   154     """
   155     Read an SLHA file, returning dictionaries of blocks and decays.
   157     Other keyword parameters are passed to readSLHA.
   158     """
   159     f = open(spcfilename, "r")
   160     rtn = readSLHA(f.read(), kwargs)
   161     f.close()
   162     return rtn
   165 def readSLHA(spcstr, ignorenobr=False):
   166     """
   167     Read an SLHA definition from a string, returning dictionaries of blocks and
   168     decays.
   170     If the ignorenobr parameter is True, do not store decay entries with a
   171     branching ratio of zero.
   172     """
   173     blocks = {}
   174     decays = {}
   175     #
   176     import re
   177     currentblock = None
   178     currentdecay = None
   179     for line in spcstr.splitlines():
   180         ## Handle (ignore) comment lines
   181         if line.startswith("#"):
   182             continue
   183         if "#" in line:
   184             line = line[:line.index("#")]
   186         ## Handle BLOCK/DECAY start lines
   187         if line.upper().startswith("BLOCK"):
   188             #print line
   189             match = re.match(r"BLOCK\s+(\w+)(\s+Q=\s*.+)?", line.upper())
   190             if not match:
   191                 continue
   192             blockname = match.group(1)
   193             qstr = match.group(2)
   194             if qstr is not None:
   195                 qstr = qstr[2:].strip()
   196             currentblock = blockname
   197             currentdecay = None
   198             blocks[blockname] = Block(blockname, q=qstr)
   199         elif line.upper().startswith("DECAY"):
   200             match = re.match(r"DECAY\s+(\d+)\s+([\d\.E+-]+).*", line.upper())
   201             if not match:
   202                 continue
   203             pdgid = int(match.group(1))
   204             width = float(match.group(2))
   205             currentblock = "DECAY"
   206             currentdecay = pdgid
   207             decays[pdgid] = Particle(pdgid, width)
   208         else:
   209             ## In-block line
   210             if currentblock is not None:
   211                 items = line.split()
   212                 if len(items) < 1:
   213                     continue
   214                 if currentblock != "DECAY":
   215                     if len(items) < 2:
   216                         ## Treat the ALPHA block differently
   217                         blocks[currentblock].value = _autotype(items[0])
   218                         blocks[currentblock].entries = _autotype(items[0])
   219                     else:
   220                         blocks[currentblock].add_entry(items)
   221                 else:
   222                     br = float(items[0])
   223                     nda = int(items[1])
   224                     ids = map(int, items[2:])
   225                     if br > 0.0 or not ignorenobr:
   226                         decays[currentdecay].add_decay(br, nda, ids)
   228     ## Try to populate Particle masses from the MASS block
   229     # print blocks.keys()
   230     try:
   231         for pid in blocks["MASS"].entries.keys():
   232             if decays.has_key(pid):
   233                 decays[pid].mass = blocks["MASS"].entries[pid]
   234     except:
   235         raise Exception("No MASS block found, from which to populate particle masses")
   237     return blocks, decays
   240 def readISAWIGFile(isafilename, **kwargs):
   241     """
   242     Read a spectrum definition from a file in the ISAWIG format, returning
   243     dictionaries of blocks and decays. While this is not an SLHA format, it is
   244     informally supported as a useful mechanism for converting ISAWIG spectra to
   245     SLHA.
   247     Other keyword parameters are passed to readSLHA.
   248     """
   249     f = open(isafilename, "r")
   250     rtn = readISAWIG(f.read(), kwargs)
   251     f.close()
   252     return rtn
   255 def readISAWIG(isastr, ignorenobr=False):
   256     """
   257     Read a spectrum definition from a string in the ISAWIG format, returning
   258     dictionaries of blocks and decays. While this is not an SLHA format, it is
   259     informally supported as a useful mechanism for converting ISAWIG spectra to
   260     SLHA.
   262     ISAWIG parsing based on the HERWIG SUSY specification format, from
   263     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   265     If the ignorenobr parameter is True, do not store decay entries with a
   266     branching ratio of zero.
   267     """
   269     ## PDG MC ID codes mapped to HERWIG SUSY ID codes, based on
   270     ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   271     HERWIGID2PDGID = {}
   272     HERWIGID2PDGID[203] =  25 ## HIGGSL0
   273     HERWIGID2PDGID[204] =  35 ## HIGGSH0
   274     HERWIGID2PDGID[205] =  36 ## HIGGSA0
   275     HERWIGID2PDGID[206] =  37 ## HIGGS+
   276     HERWIGID2PDGID[207] = -37 ## HIGGS-
   277     HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   278     HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   279     HERWIGID2PDGID[402] =  1000002 ## SSULBR
   280     HERWIGID2PDGID[408] = -1000002 ## SSUL
   281     HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   282     HERWIGID2PDGID[409] = -1000003 ## SSSL
   283     HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   284     HERWIGID2PDGID[410] = -1000004 ## SSCL
   285     HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   286     HERWIGID2PDGID[411] = -1000005 ## SSB1
   287     HERWIGID2PDGID[406] =  1000006 ## SST1BR
   288     HERWIGID2PDGID[412] = -1000006 ## SST1
   289     HERWIGID2PDGID[413] =  2000001 ## SSDR
   290     HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   291     HERWIGID2PDGID[414] =  2000002 ## SSUR
   292     HERWIGID2PDGID[420] = -2000002 ## SSURBR
   293     HERWIGID2PDGID[415] =  2000003 ## SSSR
   294     HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   295     HERWIGID2PDGID[416] =  2000004 ## SSCR
   296     HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   297     HERWIGID2PDGID[417] =  2000005 ## SSB2
   298     HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   299     HERWIGID2PDGID[418] =  2000006 ## SST2
   300     HERWIGID2PDGID[424] = -2000006 ## SST2BR
   301     HERWIGID2PDGID[425] =  1000011 ## SSEL-
   302     HERWIGID2PDGID[431] = -1000011 ## SSEL+
   303     HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   304     HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   305     HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   306     HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   307     HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   308     HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   309     HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   310     HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   311     HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   312     HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   313     HERWIGID2PDGID[437] =  2000011 ## SSEL-
   314     HERWIGID2PDGID[443] = -2000011 ## SSEL+
   315     HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   316     HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   317     HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   318     HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   319     HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   320     HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   321     HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   322     HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   323     HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   324     HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   325     HERWIGID2PDGID[449] =  1000021 ## GLUINO
   326     HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   327     HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   328     HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   329     HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   330     HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   331     HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   332     HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   333     HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   334     HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   336     blocks = {}
   337     decays = {}
   338     LINES = isastr.splitlines()
   340     def getnextvalidline():
   341         while LINES:
   342             s = LINES.pop(0).strip()
   343             ## Return None if EOF reached
   344             if len(s) == 0:
   345                 continue
   346             ## Strip comments
   347             if "#" in s:
   348                 s = s[:s.index("#")].strip()
   349             ## Return if non-empty
   350             if len(s) > 0:
   351                 return s
   353     def getnextvalidlineitems():
   354         return map(_autotype, getnextvalidline().split())
   356     ## Populate MASS block and create decaying particle objects
   357     masses = Block("MASS")
   358     numentries = int(getnextvalidline())
   359     for i in xrange(numentries):
   360         hwid, mass, lifetime = getnextvalidlineitems()
   361         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   362         pdgid = HERWIGID2PDGID.get(hwid, hwid)
   363         masses.add_entry((pdgid, mass))
   364         decays[pdgid] = Particle(pdgid, width, mass)
   365         #print pdgid, mass, width
   366     blocks["MASS"] = masses
   368     ## Populate decays
   369     for n in xrange(numentries):
   370         numdecays = int(getnextvalidline())
   371         for d in xrange(numdecays):
   372             #print n, numentries-1, d, numdecays-1
   373             decayitems = getnextvalidlineitems()
   374             hwid = decayitems[0]
   375             pdgid = HERWIGID2PDGID.get(hwid, hwid)
   376             br = decayitems[1]
   377             nme = decayitems[2]
   378             daughter_hwids = decayitems[3:]
   379             daughter_pdgids = []
   380             for hw in daughter_hwids:
   381                 if hw != 0:
   382                     daughter_pdgids.append(HERWIGID2PDGID.get(hw, hw))
   383             if not decays.has_key(pdgid):
   384                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   385                 decays[pdgid] = Particle(pdgid)
   386             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   389     ## Now the SUSY parameters
   390     TANB, ALPHAH = getnextvalidlineitems()
   391     blocks["MINPAR"] = Block("MINPAR")
   392     blocks["MINPAR"].add_entry((3, TANB))
   393     blocks["ALPHA"] = Block("ALPHA")
   394     blocks["ALPHA"].entries = ALPHAH
   395     #
   396     ## Neutralino mixing matrix
   397     blocks["NMIX"] = Block("NMIX")
   398     for i in xrange(1, 5):
   399         nmix_i = getnextvalidlineitems()
   400         for j, v in enumerate(nmix_i):
   401             blocks["NMIX"].add_entry((i, j+1, v))
   402     #
   403     ## Chargino mixing matrices V and U
   404     blocks["VMIX"] = Block("VMIX")
   405     vmix = getnextvalidlineitems()
   406     blocks["VMIX"].add_entry((1, 1, vmix[0]))
   407     blocks["VMIX"].add_entry((1, 2, vmix[1]))
   408     blocks["VMIX"].add_entry((2, 1, vmix[2]))
   409     blocks["VMIX"].add_entry((2, 2, vmix[3]))
   410     blocks["UMIX"] = Block("UMIX")
   411     umix = getnextvalidlineitems()
   412     blocks["UMIX"].add_entry((1, 1, umix[0]))
   413     blocks["UMIX"].add_entry((1, 2, umix[1]))
   414     blocks["UMIX"].add_entry((2, 1, umix[2]))
   415     blocks["UMIX"].add_entry((2, 2, umix[3]))
   416     #
   417     THETAT, THETAB, THETAL = getnextvalidlineitems()
   418     import math
   419     blocks["STOPMIX"] = Block("STOPMIX")
   420     blocks["STOPMIX"].add_entry((1, 1,  math.cos(THETAT)))
   421     blocks["STOPMIX"].add_entry((1, 2, -math.sin(THETAT)))
   422     blocks["STOPMIX"].add_entry((2, 1,  math.sin(THETAT)))
   423     blocks["STOPMIX"].add_entry((2, 2,  math.cos(THETAT)))
   424     blocks["SBOTMIX"] = Block("SBOTMIX")
   425     blocks["SBOTMIX"].add_entry((1, 1,  math.cos(THETAB)))
   426     blocks["SBOTMIX"].add_entry((1, 2, -math.sin(THETAB)))
   427     blocks["SBOTMIX"].add_entry((2, 1,  math.sin(THETAB)))
   428     blocks["SBOTMIX"].add_entry((2, 2,  math.cos(THETAB)))
   429     blocks["STAUMIX"] = Block("STAUMIX")
   430     blocks["STAUMIX"].add_entry((1, 1,  math.cos(THETAL)))
   431     blocks["STAUMIX"].add_entry((1, 2, -math.sin(THETAL)))
   432     blocks["STAUMIX"].add_entry((2, 1,  math.sin(THETAL)))
   433     blocks["STAUMIX"].add_entry((2, 2,  math.cos(THETAL)))
   434     #
   435     ATSS, ABSS, ALSS = getnextvalidlineitems()
   436     blocks["AU"] = Block("AU")
   437     blocks["AU"].add_entry((3, 3, ATSS))
   438     blocks["AD"] = Block("AD")
   439     blocks["AD"].add_entry((3, 3, ABSS))
   440     blocks["AE"] = Block("AE")
   441     blocks["AE"].add_entry((3, 3, ALSS))
   442     #
   443     MUSS = getnextvalidlineitems()[0]
   444     blocks["MINPAR"].add_entry((4, MUSS))
   445     #
   446     return blocks, decays
   449 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
   450     """
   451     Write an SLHA file from the supplied blocks and decays dicts.
   453     Other keyword parameters are passed to writeSLHA.
   454     """
   455     f = open(spcfilename, "w")
   456     f.write(writeSLHA(blocks, decays, kwargs))
   457     f.close()
   460 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   463 def writeSLHA(blocks, decays, ignorenobr=False):
   464     """
   465     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   466     """
   467     sep = "   "
   468     out = ""
   469     def dict_hier_strs(d, s=""):
   470         if type(d) is dict:
   471             for k, v in sorted(d.iteritems()):
   472                 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
   473                     yield s2
   474         else:
   475             yield s + sep + _autostr(d)
   476     ## Blocks
   477     for bname, b in sorted(blocks.iteritems()):
   478         namestr = b.name
   479         if b.q is not None:
   480             namestr += " Q= %e" % b.q
   481         out += "BLOCK %s\n" % namestr
   482         for s in dict_hier_strs(b.entries):
   483             out += sep + s + "\n"
   484         out += "\n"
   485     ## Decays
   486     for pid, particle in sorted(decays.iteritems()):
   487         out += "DECAY %d %e\n" % (particle.pid, particle.totalwidth or -1)
   488         for d in sorted(particle.decays):
   489             if d.br > 0.0 or not ignorenobr:
   490                 products_str = "   ".join(map(str, d.ids))
   491                 out += sep + "%e" % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
   492         out += "\n"
   493     return out
   497 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
   498     """
   499     Write an ISAWIG file from the supplied blocks and decays dicts.
   501     Other keyword parameters are passed to writeISAWIG.
   503     TODO: Handle RPV SUSY
   504     """
   505     f = open(isafilename, "w")
   506     f.write(writeISAWIG(blocks, decays, kwargs))
   507     f.close()
   510 def writeISAWIG(blocks, decays, ignorenobr=False):
   511     """
   512     Return an ISAWIG definition as a string, from the supplied blocks and decays dicts.
   514     ISAWIG parsing based on the HERWIG SUSY specification format, from
   515     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   517     If the ignorenobr parameter is True, do not write decay entries with a
   518     branching ratio of zero.
   519     """
   520     ## PDG MC ID codes mapped to HERWIG SUSY ID codes, based on
   521     ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   522     PDGID2HERWIGID = {}
   523     PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (ADDED)
   524     PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   525     PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   526     PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   527     PDGID2HERWIGID[      37] = 206 ## HIGGS+
   528     PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   529     PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   530     PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   531     PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   532     PDGID2HERWIGID[-1000002] = 408 ## SSUL
   533     PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   534     PDGID2HERWIGID[-1000003] = 409 ## SSSL
   535     PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   536     PDGID2HERWIGID[-1000004] = 410 ## SSCL
   537     PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   538     PDGID2HERWIGID[-1000005] = 411 ## SSB1
   539     PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   540     PDGID2HERWIGID[-1000006] = 412 ## SST1
   541     PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   542     PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   543     PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   544     PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   545     PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   546     PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   547     PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   548     PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   549     PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   550     PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   551     PDGID2HERWIGID[ 2000006] = 418 ## SST2
   552     PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   553     PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   554     PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   555     PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   556     PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   557     PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   558     PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   559     PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   560     PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   561     PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   562     PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   563     PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   564     PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   565     PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   566     PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   567     PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   568     PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   569     PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   570     PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   571     PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   572     PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   573     PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   574     PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   575     PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   576     PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   577     PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   578     PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   579     PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   580     PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   581     PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   582     PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   583     PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   584     PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   585     PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   586     PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   588     masses = blocks["MASS"].entries
   590     ## Init output string
   591     out = ""
   593     ## First write out masses section:
   594     ##   Number of SUSY + top particles
   595     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   596     ##   repeated for each particle
   597     ## IDHW is the HERWIG identity code.
   598     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   599     massout = ""
   600     for pid in masses.keys():
   601         lifetime = -1
   602         try:
   603             width = decays[pid].totalwidth
   604             if width and width > 0:
   605                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   606         except:
   607             pass
   608         massout += "%d %e %e\n" % (PDGID2HERWIGID.get(pid, pid), masses[pid], lifetime)
   609     out += "%d\n" % massout.count("\n")
   610     out += massout
   612     assert(len(masses) == len(decays))
   614     ## Next each particles decay modes together with their branching ratios and matrix element codes
   615     ##   Number of decay modes for a given particle (IDK)
   616     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   617     ##     repeated for each mode.
   618     ##   Repeated for each particle.
   619     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   620     ## the decay mode. NME is a code for the matrix element to be used, either from the
   621     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   622     for i, pid in enumerate(decays.keys()):
   623         # if not decays.has_key(pid):
   624         #     continue
   625         hwid = PDGID2HERWIGID.get(pid, pid)
   626         decayout = ""
   627         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   628         for i_d, d in enumerate(decays[pid].decays):
   629             ## Skip decay if it has no branching ratio
   630             if ignorenobr and d.br == 0:
   631                 continue
   633             ## Identify decay matrix element to use
   634             ## From std HW docs, or from this pair:
   635             ## Two new matrix element codes have been added for these new decays:
   636             ##    NME =	200 	3 body top quark via charged Higgs
   637             ##    	300 	3 body R-parity violating gaugino and gluino decays
   638             nme = 0
   639             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   640             if abs(pid) in (6, 12):
   641                 nme = 100
   642             ## Extra SUSY MEs
   643             if len(d.ids) == 3:
   644                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   645                 pass
   646             decayout += "%d %e %d " % (hwid, d.br, nme)
   648             def is_quark(pid):
   649                 return (abs(pid) in range(1, 7))
   651             def is_lepton(pid):
   652                 return (abs(pid) in range(11, 17))
   654             def is_squark(pid):
   655                 if abs(pid) in range(1000001, 1000007):
   656                     return True
   657                 if abs(pid) in range(2000001, 2000007):
   658                     return True
   659                 return False
   661             def is_slepton(pid):
   662                 if abs(pid) in range(1000011, 1000017):
   663                     return True
   664                 if abs(pid) in range(2000011, 2000016, 2):
   665                     return True
   666                 return False
   668             def is_gaugino(pid):
   669                 if abs(pid) in range(1000022, 1000026):
   670                     return True
   671                 if abs(pid) in (1000035, 1000037):
   672                     return True
   673                 return False
   675             def is_susy(pid):
   676                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   678             absids = map(abs, d.ids)
   680             ## Order decay products as required by HERWIG
   681             ## Top
   682             if abs(pid) == 6:
   683                 def cmp_bottomlast(a, b):
   684                     """Comparison function which always puts b/bbar last"""
   685                     if abs(a) == 5:
   686                         return True
   687                     if abs(b) == 5:
   688                         return False
   689                     return cmp(a, b)
   690                 if len(absids) == 2:
   691                     ## 2 body mode, to Higgs: Higgs; Bottom
   692                     if (25 in absids or 26 in absids) and 5 in absids:
   693                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   694                 elif len(absids) == 3:
   695                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   696                     if 37 in absids or 23 in absids:
   697                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   698             ## Gluino
   699             elif abs(pid) == 1000021:
   700                 if len(absids) == 2:
   701                     ## 2 body mode
   702                     ## without gluon: any order
   703                     ## with gluon: gluon; colour neutral
   704                     if 21 in absids:
   705                         def cmp_gluonfirst(a, b):
   706                             """Comparison function which always puts gluon first"""
   707                             if a == 21:
   708                                 return False
   709                             if b == 21:
   710                                 return True
   711                             return cmp(a, b)
   712                         d.ids = sorted(d.ids, key=cmp_gluonfirst)
   713                 elif len(absids) == 3:
   714                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   715                     def cmp_quarkslast(a, b):
   716                         """Comparison function which always puts quarks last"""
   717                         if is_quark(a):
   718                             return True
   719                         if is_quark(b):
   720                             return False
   721                         return cmp(a, b)
   722                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   723             ## Squark/Slepton
   724             elif is_squark(pid) or is_slepton(pid):
   725                 def cmp_susy_quark_lepton(a, b):
   726                     if is_susy(a):
   727                         return False
   728                     if is_susy(b):
   729                         return True
   730                     if is_quark(a):
   731                         return False
   732                     if is_quark(b):
   733                         return True
   734                     return cmp(a, b)
   735                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   736                 ##                                                      Gluino       lepton
   737                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   738                 ## Squark
   739                 ##   2 body modes: Lepton Number Violated               quark     lepton
   740                 ##                 Baryon number violated               quark     quark
   741                 ## Slepton
   742                 ##   2 body modes: Lepton Number Violated               q or qbar
   743                 d.ids = sorted(d.ids, key=cmp_bottomlast)
   744             ## Higgs
   745             elif pid in (25, 26):
   746                 # TODO: Includes SUSY Higgses?
   747                 ## Higgs
   748                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   749                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   750                 ##   3 body modes:                                      colour neutral       q or qbar
   751                 if len(absids) == 3:
   752                     def cmp_quarkslast(a, b):
   753                         """Comparison function which always puts quarks last"""
   754                         if is_quark(a):
   755                             return True
   756                         if is_quark(b):
   757                             return False
   758                         return cmp(a, b)
   759                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   760             elif is_gaugino(pid):
   761                 # TODO: Is there actually anything to do here?
   762                 ## Gaugino
   763                 ##   2 body modes: Squark-quark                         q or sq
   764                 ##                 Slepton-lepton                       l or sl
   765                 ##
   766                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
   767                 ##                                                                           l or lbar
   768                 if len(absids) == 3:
   769                     def cmp_quarkslast(a, b):
   770                         """Comparison function which always puts quarks last"""
   771                         if is_quark(a):
   772                             return True
   773                         if is_quark(b):
   774                             return False
   775                         return cmp(a, b)
   776                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   778             # TODO: Gaugino/Gluino
   779             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
   781             ## Pad out IDs list with zeros
   782             ids = [0,0,0,0,0]
   783             for i, pid in enumerate(d.ids):
   784                 ids[i] = pid
   785             ids = map(str, ids)
   786             decayout += " ".join(ids) + "\n"
   787         decayout = "%d\n" % decayout.count("\n") + decayout
   788         out += decayout
   790     ## Now the SUSY parameters
   791     ## TANB, ALPHAH:
   792     out += "%e %e\n" % (blocks["MINPAR"].entries[3], blocks["ALPHA"].entries)
   793     ## Neutralino mixing matrix
   794     nmix = blocks["NMIX"].entries
   795     for i in xrange(1, 5):
   796         out += "%e %e %e %e\n" % (nmix[i][1], nmix[i][2], nmix[i][3], nmix[i][4])
   797     ## Chargino mixing matrices V and U
   798     vmix = blocks["VMIX"].entries
   799     out += "%e %e %e %e\n" % (vmix[1][1], vmix[1][2], vmix[2][1], vmix[2][2])
   800     umix = blocks["UMIX"].entries
   801     out += "%e %e %e %e\n" % (umix[1][1], umix[1][2], umix[2][1], umix[2][2])
   802     # THETAT,THETAB,THETAL
   803     import math
   804     out += "%e %e %e\n" % (math.acos(blocks["STOPMIX"].entries[1][1]),
   805                            math.acos(blocks["SBOTMIX"].entries[1][1]),
   806                            math.acos(blocks["STAUMIX"].entries[1][1]))
   807     # ATSS,ABSS,ALSS
   808     out += "%e %e %e\n" % (blocks["AU"].entries[3][3],
   809                            blocks["AD"].entries[3][3],
   810                            blocks["AE"].entries[3][3])
   811     # MUSS == sign(mu)
   812     out += "%f\n" % blocks["MINPAR"].entries[4]
   814     ## TODO: Handle RPV SUSY
   816     return out
   820 if __name__ == "__main__":
   821     import sys
   822     for a in sys.argv[1:]:
   823         if a.endswith(".isa"):
   824             blocks, decays = readISAWIGFile(a)
   825         else:
   826             blocks, decays = readSLHAFile(a)
   828         for bname, b in sorted(blocks.iteritems()):
   829             print b
   830             print
   832         print blocks.keys()
   834         print blocks["MASS"].entries[25]
   835         print
   837         for p in sorted(decays.values()):
   838             print p
   839             print
   841         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial