pyslha.py

Sat, 26 Feb 2011 23:21:12 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sat, 26 Feb 2011 23:21:12 +0100
changeset 120
bffefe12df80
parent 118
5ab517422296
child 123
c85e29bc13c4
permissions
-rw-r--r--

Remove --show-gravitino: we'll always show it from now

     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.8"
    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
   256 def readISAWIG(isastr, ignorenobr=False):
   257     """
   258     Read a spectrum definition from a string in the ISAWIG format, returning
   259     dictionaries of blocks and decays. While this is not an SLHA format, it is
   260     informally supported as a useful mechanism for converting ISAWIG spectra to
   261     SLHA.
   263     ISAWIG parsing based on the HERWIG SUSY specification format, from
   264     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   266     If the ignorenobr parameter is True, do not store decay entries with a
   267     branching ratio of zero.
   268     """
   270     ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   271     ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   272     ## + the IDPDG array and section 4.13 of the HERWIG manual.
   273     HERWIGID2PDGID = {}
   274     HERWIGID2PDGID[7]   = -1
   275     HERWIGID2PDGID[8]   = -2
   276     HERWIGID2PDGID[9]   = -3
   277     HERWIGID2PDGID[10]  = -4
   278     HERWIGID2PDGID[11]  = -5
   279     HERWIGID2PDGID[12]  = -6
   280     HERWIGID2PDGID[13]  =  21
   281     HERWIGID2PDGID[59]  =  22
   282     HERWIGID2PDGID[121] =  11
   283     HERWIGID2PDGID[122] =  12
   284     HERWIGID2PDGID[123] =  13
   285     HERWIGID2PDGID[124] =  14
   286     HERWIGID2PDGID[125] =  15
   287     HERWIGID2PDGID[126] =  16
   288     HERWIGID2PDGID[127] = -11
   289     HERWIGID2PDGID[128] = -12
   290     HERWIGID2PDGID[129] = -13
   291     HERWIGID2PDGID[130] = -14
   292     HERWIGID2PDGID[131] = -15
   293     HERWIGID2PDGID[132] = -16
   294     HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   295     HERWIGID2PDGID[204] =  35 ## HIGGSH0
   296     HERWIGID2PDGID[205] =  36 ## HIGGSA0
   297     HERWIGID2PDGID[206] =  37 ## HIGGS+
   298     HERWIGID2PDGID[207] = -37 ## HIGGS-
   299     HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   300     HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   301     HERWIGID2PDGID[402] =  1000002 ## SSULBR
   302     HERWIGID2PDGID[408] = -1000002 ## SSUL
   303     HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   304     HERWIGID2PDGID[409] = -1000003 ## SSSL
   305     HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   306     HERWIGID2PDGID[410] = -1000004 ## SSCL
   307     HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   308     HERWIGID2PDGID[411] = -1000005 ## SSB1
   309     HERWIGID2PDGID[406] =  1000006 ## SST1BR
   310     HERWIGID2PDGID[412] = -1000006 ## SST1
   311     HERWIGID2PDGID[413] =  2000001 ## SSDR
   312     HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   313     HERWIGID2PDGID[414] =  2000002 ## SSUR
   314     HERWIGID2PDGID[420] = -2000002 ## SSURBR
   315     HERWIGID2PDGID[415] =  2000003 ## SSSR
   316     HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   317     HERWIGID2PDGID[416] =  2000004 ## SSCR
   318     HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   319     HERWIGID2PDGID[417] =  2000005 ## SSB2
   320     HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   321     HERWIGID2PDGID[418] =  2000006 ## SST2
   322     HERWIGID2PDGID[424] = -2000006 ## SST2BR
   323     HERWIGID2PDGID[425] =  1000011 ## SSEL-
   324     HERWIGID2PDGID[431] = -1000011 ## SSEL+
   325     HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   326     HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   327     HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   328     HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   329     HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   330     HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   331     HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   332     HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   333     HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   334     HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   335     HERWIGID2PDGID[437] =  2000011 ## SSEL-
   336     HERWIGID2PDGID[443] = -2000011 ## SSEL+
   337     HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   338     HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   339     HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   340     HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   341     HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   342     HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   343     HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   344     HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   345     HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   346     HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   347     HERWIGID2PDGID[449] =  1000021 ## GLUINO
   348     HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   349     HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   350     HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   351     HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   352     HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   353     HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   354     HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   355     HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   356     HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   358     blocks = {}
   359     decays = {}
   360     LINES = isastr.splitlines()
   362     def getnextvalidline():
   363         while LINES:
   364             s = LINES.pop(0).strip()
   365             ## Return None if EOF reached
   366             if len(s) == 0:
   367                 continue
   368             ## Strip comments
   369             if "#" in s:
   370                 s = s[:s.index("#")].strip()
   371             ## Return if non-empty
   372             if len(s) > 0:
   373                 return s
   375     def getnextvalidlineitems():
   376         return map(_autotype, getnextvalidline().split())
   378     ## Populate MASS block and create decaying particle objects
   379     masses = Block("MASS")
   380     numentries = int(getnextvalidline())
   381     for i in xrange(numentries):
   382         hwid, mass, lifetime = getnextvalidlineitems()
   383         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   384         pdgid = HERWIGID2PDGID.get(hwid, hwid)
   385         masses.add_entry((pdgid, mass))
   386         decays[pdgid] = Particle(pdgid, width, mass)
   387         #print pdgid, mass, width
   388     blocks["MASS"] = masses
   390     ## Populate decays
   391     for n in xrange(numentries):
   392         numdecays = int(getnextvalidline())
   393         for d in xrange(numdecays):
   394             #print n, numentries-1, d, numdecays-1
   395             decayitems = getnextvalidlineitems()
   396             hwid = decayitems[0]
   397             pdgid = HERWIGID2PDGID.get(hwid, hwid)
   398             br = decayitems[1]
   399             nme = decayitems[2]
   400             daughter_hwids = decayitems[3:]
   401             daughter_pdgids = []
   402             for hw in daughter_hwids:
   403                 if hw != 0:
   404                     daughter_pdgids.append(HERWIGID2PDGID.get(hw, hw))
   405             if not decays.has_key(pdgid):
   406                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   407                 decays[pdgid] = Particle(pdgid)
   408             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   411     ## Now the SUSY parameters
   412     TANB, ALPHAH = getnextvalidlineitems()
   413     blocks["MINPAR"] = Block("MINPAR")
   414     blocks["MINPAR"].add_entry((3, TANB))
   415     blocks["ALPHA"] = Block("ALPHA")
   416     blocks["ALPHA"].entries = ALPHAH
   417     #
   418     ## Neutralino mixing matrix
   419     blocks["NMIX"] = Block("NMIX")
   420     for i in xrange(1, 5):
   421         nmix_i = getnextvalidlineitems()
   422         for j, v in enumerate(nmix_i):
   423             blocks["NMIX"].add_entry((i, j+1, v))
   424     #
   425     ## Chargino mixing matrices V and U
   426     blocks["VMIX"] = Block("VMIX")
   427     vmix = getnextvalidlineitems()
   428     blocks["VMIX"].add_entry((1, 1, vmix[0]))
   429     blocks["VMIX"].add_entry((1, 2, vmix[1]))
   430     blocks["VMIX"].add_entry((2, 1, vmix[2]))
   431     blocks["VMIX"].add_entry((2, 2, vmix[3]))
   432     blocks["UMIX"] = Block("UMIX")
   433     umix = getnextvalidlineitems()
   434     blocks["UMIX"].add_entry((1, 1, umix[0]))
   435     blocks["UMIX"].add_entry((1, 2, umix[1]))
   436     blocks["UMIX"].add_entry((2, 1, umix[2]))
   437     blocks["UMIX"].add_entry((2, 2, umix[3]))
   438     #
   439     THETAT, THETAB, THETAL = getnextvalidlineitems()
   440     import math
   441     blocks["STOPMIX"] = Block("STOPMIX")
   442     blocks["STOPMIX"].add_entry((1, 1,  math.cos(THETAT)))
   443     blocks["STOPMIX"].add_entry((1, 2, -math.sin(THETAT)))
   444     blocks["STOPMIX"].add_entry((2, 1,  math.sin(THETAT)))
   445     blocks["STOPMIX"].add_entry((2, 2,  math.cos(THETAT)))
   446     blocks["SBOTMIX"] = Block("SBOTMIX")
   447     blocks["SBOTMIX"].add_entry((1, 1,  math.cos(THETAB)))
   448     blocks["SBOTMIX"].add_entry((1, 2, -math.sin(THETAB)))
   449     blocks["SBOTMIX"].add_entry((2, 1,  math.sin(THETAB)))
   450     blocks["SBOTMIX"].add_entry((2, 2,  math.cos(THETAB)))
   451     blocks["STAUMIX"] = Block("STAUMIX")
   452     blocks["STAUMIX"].add_entry((1, 1,  math.cos(THETAL)))
   453     blocks["STAUMIX"].add_entry((1, 2, -math.sin(THETAL)))
   454     blocks["STAUMIX"].add_entry((2, 1,  math.sin(THETAL)))
   455     blocks["STAUMIX"].add_entry((2, 2,  math.cos(THETAL)))
   456     #
   457     ATSS, ABSS, ALSS = getnextvalidlineitems()
   458     blocks["AU"] = Block("AU")
   459     blocks["AU"].add_entry((3, 3, ATSS))
   460     blocks["AD"] = Block("AD")
   461     blocks["AD"].add_entry((3, 3, ABSS))
   462     blocks["AE"] = Block("AE")
   463     blocks["AE"].add_entry((3, 3, ALSS))
   464     #
   465     MUSS = getnextvalidlineitems()[0]
   466     blocks["MINPAR"].add_entry((4, MUSS))
   467     #
   468     return blocks, decays
   471 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
   472     """
   473     Write an SLHA file from the supplied blocks and decays dicts.
   475     Other keyword parameters are passed to writeSLHA.
   476     """
   477     f = open(spcfilename, "w")
   478     f.write(writeSLHA(blocks, decays, kwargs))
   479     f.close()
   482 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   485 def writeSLHA(blocks, decays, ignorenobr=False):
   486     """
   487     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   488     """
   489     sep = "   "
   490     out = ""
   491     def dict_hier_strs(d, s=""):
   492         if type(d) is dict:
   493             for k, v in sorted(d.iteritems()):
   494                 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
   495                     yield s2
   496         else:
   497             yield s + sep + _autostr(d)
   498     ## Blocks
   499     for bname, b in sorted(blocks.iteritems()):
   500         namestr = b.name
   501         if b.q is not None:
   502             namestr += " Q= %e" % b.q
   503         out += "BLOCK %s\n" % namestr
   504         for s in dict_hier_strs(b.entries):
   505             out += sep + s + "\n"
   506         out += "\n"
   507     ## Decays
   508     for pid, particle in sorted(decays.iteritems()):
   509         out += "DECAY %d %e\n" % (particle.pid, particle.totalwidth or -1)
   510         for d in sorted(particle.decays):
   511             if d.br > 0.0 or not ignorenobr:
   512                 products_str = "   ".join(map(str, d.ids))
   513                 out += sep + "%e" % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
   514         out += "\n"
   515     return out
   519 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
   520     """
   521     Write an ISAWIG file from the supplied blocks and decays dicts.
   523     Other keyword parameters are passed to writeISAWIG.
   525     TODO: Handle RPV SUSY
   526     """
   527     f = open(isafilename, "w")
   528     f.write(writeISAWIG(blocks, decays, kwargs))
   529     f.close()
   532 def writeISAWIG(blocks, decays, ignorenobr=False):
   533     """
   534     Return an ISAWIG definition as a string, from the supplied blocks and decays dicts.
   536     ISAWIG parsing based on the HERWIG SUSY specification format, from
   537     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   539     If the ignorenobr parameter is True, do not write decay entries with a
   540     branching ratio of zero.
   541     """
   542     ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   543     ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   544     ## + the IDPDG array and section 4.13 of the HERWIG manual.
   545     PDGID2HERWIGID = {}
   546     PDGID2HERWIGID[      -1] = 7
   547     PDGID2HERWIGID[      -2] = 8
   548     PDGID2HERWIGID[      -3] = 9
   549     PDGID2HERWIGID[      -4] = 10
   550     PDGID2HERWIGID[      -5] = 11
   551     PDGID2HERWIGID[      -6] = 12
   552     PDGID2HERWIGID[      21] = 13
   553     PDGID2HERWIGID[      22] = 59
   554     PDGID2HERWIGID[      11] = 121
   555     PDGID2HERWIGID[      12] = 122
   556     PDGID2HERWIGID[      13] = 123
   557     PDGID2HERWIGID[      14] = 124
   558     PDGID2HERWIGID[      15] = 125
   559     PDGID2HERWIGID[      16] = 126
   560     PDGID2HERWIGID[     -11] = 127
   561     PDGID2HERWIGID[     -12] = 128
   562     PDGID2HERWIGID[     -13] = 129
   563     PDGID2HERWIGID[     -14] = 130
   564     PDGID2HERWIGID[     -15] = 131
   565     PDGID2HERWIGID[     -16] = 132
   566     PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW)
   567     PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   568     PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   569     PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   570     PDGID2HERWIGID[      37] = 206 ## HIGGS+
   571     PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   572     PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   573     PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   574     PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   575     PDGID2HERWIGID[-1000002] = 408 ## SSUL
   576     PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   577     PDGID2HERWIGID[-1000003] = 409 ## SSSL
   578     PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   579     PDGID2HERWIGID[-1000004] = 410 ## SSCL
   580     PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   581     PDGID2HERWIGID[-1000005] = 411 ## SSB1
   582     PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   583     PDGID2HERWIGID[-1000006] = 412 ## SST1
   584     PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   585     PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   586     PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   587     PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   588     PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   589     PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   590     PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   591     PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   592     PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   593     PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   594     PDGID2HERWIGID[ 2000006] = 418 ## SST2
   595     PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   596     PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   597     PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   598     PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   599     PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   600     PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   601     PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   602     PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   603     PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   604     PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   605     PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   606     PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   607     PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   608     PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   609     PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   610     PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   611     PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   612     PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   613     PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   614     PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   615     PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   616     PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   617     PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   618     PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   619     PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   620     PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   621     PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   622     PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   623     PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   624     PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   625     PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   626     PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   627     PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   628     PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   629     PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   631     masses = blocks["MASS"].entries
   633     ## Init output string
   634     out = ""
   636     ## First write out masses section:
   637     ##   Number of SUSY + top particles
   638     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   639     ##   repeated for each particle
   640     ## IDHW is the HERWIG identity code.
   641     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   642     massout = ""
   643     for pid in masses.keys():
   644         lifetime = -1
   645         try:
   646             width = decays[pid].totalwidth
   647             if width and width > 0:
   648                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   649         except:
   650             pass
   651         massout += "%d %e %e\n" % (PDGID2HERWIGID.get(pid, pid), masses[pid], lifetime)
   652     out += "%d\n" % massout.count("\n")
   653     out += massout
   655     assert(len(masses) == len(decays))
   657     ## Next each particles decay modes together with their branching ratios and matrix element codes
   658     ##   Number of decay modes for a given particle (IDK)
   659     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   660     ##     repeated for each mode.
   661     ##   Repeated for each particle.
   662     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   663     ## the decay mode. NME is a code for the matrix element to be used, either from the
   664     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   665     for i, pid in enumerate(decays.keys()):
   666         # if not decays.has_key(pid):
   667         #     continue
   668         hwid = PDGID2HERWIGID.get(pid, pid)
   669         decayout = ""
   670         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   671         for i_d, d in enumerate(decays[pid].decays):
   672             ## Skip decay if it has no branching ratio
   673             if ignorenobr and d.br == 0:
   674                 continue
   676             ## Identify decay matrix element to use
   677             ## From std HW docs, or from this pair:
   678             ## Two new matrix element codes have been added for these new decays:
   679             ##    NME =	200 	3 body top quark via charged Higgs
   680             ##    	300 	3 body R-parity violating gaugino and gluino decays
   681             nme = 0
   682             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   683             if abs(pid) in (6, 12):
   684                 nme = 100
   685             ## Extra SUSY MEs
   686             if len(d.ids) == 3:
   687                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   688                 pass
   689             decayout += "%d %e %d " % (hwid, d.br, nme)
   691             def is_quark(pid):
   692                 return (abs(pid) in range(1, 7))
   694             def is_lepton(pid):
   695                 return (abs(pid) in range(11, 17))
   697             def is_squark(pid):
   698                 if abs(pid) in range(1000001, 1000007):
   699                     return True
   700                 if abs(pid) in range(2000001, 2000007):
   701                     return True
   702                 return False
   704             def is_slepton(pid):
   705                 if abs(pid) in range(1000011, 1000017):
   706                     return True
   707                 if abs(pid) in range(2000011, 2000016, 2):
   708                     return True
   709                 return False
   711             def is_gaugino(pid):
   712                 if abs(pid) in range(1000022, 1000026):
   713                     return True
   714                 if abs(pid) in (1000035, 1000037):
   715                     return True
   716                 return False
   718             def is_susy(pid):
   719                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   721             absids = map(abs, d.ids)
   723             ## Order decay products as required by HERWIG
   724             ## Top
   725             if abs(pid) == 6:
   726                 def cmp_bottomlast(a, b):
   727                     """Comparison function which always puts b/bbar last"""
   728                     if abs(a) == 5:
   729                         return True
   730                     if abs(b) == 5:
   731                         return False
   732                     return cmp(a, b)
   733                 if len(absids) == 2:
   734                     ## 2 body mode, to Higgs: Higgs; Bottom
   735                     if (25 in absids or 26 in absids) and 5 in absids:
   736                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   737                 elif len(absids) == 3:
   738                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   739                     if 37 in absids or 23 in absids:
   740                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   741             ## Gluino
   742             elif abs(pid) == 1000021:
   743                 if len(absids) == 2:
   744                     ## 2 body mode
   745                     ## without gluon: any order
   746                     ## with gluon: gluon; colour neutral
   747                     if 21 in absids:
   748                         def cmp_gluonfirst(a, b):
   749                             """Comparison function which always puts gluon first"""
   750                             if a == 21:
   751                                 return False
   752                             if b == 21:
   753                                 return True
   754                             return cmp(a, b)
   755                         d.ids = sorted(d.ids, key=cmp_gluonfirst)
   756                 elif len(absids) == 3:
   757                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   758                     def cmp_quarkslast(a, b):
   759                         """Comparison function which always puts quarks last"""
   760                         if is_quark(a):
   761                             return True
   762                         if is_quark(b):
   763                             return False
   764                         return cmp(a, b)
   765                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   766             ## Squark/Slepton
   767             elif is_squark(pid) or is_slepton(pid):
   768                 def cmp_susy_quark_lepton(a, b):
   769                     if is_susy(a):
   770                         return False
   771                     if is_susy(b):
   772                         return True
   773                     if is_quark(a):
   774                         return False
   775                     if is_quark(b):
   776                         return True
   777                     return cmp(a, b)
   778                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   779                 ##                                                      Gluino       lepton
   780                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   781                 ## Squark
   782                 ##   2 body modes: Lepton Number Violated               quark     lepton
   783                 ##                 Baryon number violated               quark     quark
   784                 ## Slepton
   785                 ##   2 body modes: Lepton Number Violated               q or qbar
   786                 d.ids = sorted(d.ids, key=cmp_bottomlast)
   787             ## Higgs
   788             elif pid in (25, 26):
   789                 # TODO: Includes SUSY Higgses?
   790                 ## Higgs
   791                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   792                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   793                 ##   3 body modes:                                      colour neutral       q or qbar
   794                 if len(absids) == 3:
   795                     def cmp_quarkslast(a, b):
   796                         """Comparison function which always puts quarks last"""
   797                         if is_quark(a):
   798                             return True
   799                         if is_quark(b):
   800                             return False
   801                         return cmp(a, b)
   802                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   803             elif is_gaugino(pid):
   804                 # TODO: Is there actually anything to do here?
   805                 ## Gaugino
   806                 ##   2 body modes: Squark-quark                         q or sq
   807                 ##                 Slepton-lepton                       l or sl
   808                 ##
   809                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
   810                 ##                                                                           l or lbar
   811                 if len(absids) == 3:
   812                     def cmp_quarkslast(a, b):
   813                         """Comparison function which always puts quarks last"""
   814                         if is_quark(a):
   815                             return True
   816                         if is_quark(b):
   817                             return False
   818                         return cmp(a, b)
   819                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   821             # TODO: Gaugino/Gluino
   822             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
   824             ## Pad out IDs list with zeros
   825             ids = [0,0,0,0,0]
   826             for i, pid in enumerate(d.ids):
   827                 ids[i] = pid
   828             ids = map(str, ids)
   829             decayout += " ".join(ids) + "\n"
   830         decayout = "%d\n" % decayout.count("\n") + decayout
   831         out += decayout
   833     ## Now the SUSY parameters
   834     ## TANB, ALPHAH:
   835     out += "%e %e\n" % (blocks["MINPAR"].entries[3], blocks["ALPHA"].entries)
   836     ## Neutralino mixing matrix
   837     nmix = blocks["NMIX"].entries
   838     for i in xrange(1, 5):
   839         out += "%e %e %e %e\n" % (nmix[i][1], nmix[i][2], nmix[i][3], nmix[i][4])
   840     ## Chargino mixing matrices V and U
   841     vmix = blocks["VMIX"].entries
   842     out += "%e %e %e %e\n" % (vmix[1][1], vmix[1][2], vmix[2][1], vmix[2][2])
   843     umix = blocks["UMIX"].entries
   844     out += "%e %e %e %e\n" % (umix[1][1], umix[1][2], umix[2][1], umix[2][2])
   845     # THETAT,THETAB,THETAL
   846     import math
   847     out += "%e %e %e\n" % (math.acos(blocks["STOPMIX"].entries[1][1]),
   848                            math.acos(blocks["SBOTMIX"].entries[1][1]),
   849                            math.acos(blocks["STAUMIX"].entries[1][1]))
   850     # ATSS,ABSS,ALSS
   851     out += "%e %e %e\n" % (blocks["AU"].entries[3][3],
   852                            blocks["AD"].entries[3][3],
   853                            blocks["AE"].entries[3][3])
   854     # MUSS == sign(mu)
   855     out += "%f\n" % blocks["MINPAR"].entries[4]
   857     ## TODO: Handle RPV SUSY
   859     return out
   863 if __name__ == "__main__":
   864     import sys
   865     for a in sys.argv[1:]:
   866         if a.endswith(".isa"):
   867             blocks, decays = readISAWIGFile(a)
   868         else:
   869             blocks, decays = readSLHAFile(a)
   871         for bname, b in sorted(blocks.iteritems()):
   872             print b
   873             print
   875         print blocks.keys()
   877         print blocks["MASS"].entries[25]
   878         print
   880         for p in sorted(decays.values()):
   881             print p
   882             print
   884         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial