pyslha.py

Sun, 10 Apr 2011 12:18:08 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sun, 10 Apr 2011 12:18:08 +0100
changeset 140
81b5c01084f4
parent 139
c50a904f9cb4
child 141
869c16f9093b
permissions
-rw-r--r--

Adding little test script and bumping version for release (one little typo bug to fix first)

     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, EPS and PNG via LaTeX and the TikZ
    15 graphics package, and as LaTeX/TikZ source for direct embedding into documents or
    16 user-tweaking of the generated output.
    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.2.3"
    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
    87     def __repr__(self):
    88         return self.__str__()
    91 class Decay(object):
    92     """
    93     Object representing a decay entry on a particle decribed by the SLHA file.
    94     'Decay' objects are not a direct representation of a DECAY block in an SLHA
    95     file... that role, somewhat confusingly, is taken by the Particle class.
    97     Decay objects have three properties: a branching ratio, br, an nda number
    98     (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
    99     decay occurs. The PDG ID of the particle whose decay this represents may
   100     also be stored, but this is normally known via the Particle in which the
   101     decay is stored.
   102     """
   103     def __init__(self, br, nda, ids, parentid=None):
   104         self.parentid = parentid
   105         self.br = br
   106         self.nda = nda
   107         self.ids = ids
   108         assert(self.nda == len(self.ids))
   110     def __cmp__(self, other):
   111         return cmp(other.br, self.br)
   113     def __str__(self):
   114         return "%e %s" % (self.br, self.ids)
   116     def __repr__(self):
   117         return self.__str__()
   120 class Particle(object):
   121     """
   122     Representation of a single, specific particle, decay block from an SLHA
   123     file.  These objects are not themselves called 'Decay', since that concept
   124     applies more naturally to the various decays found inside this
   125     object. Particle classes store the PDG ID (pid) of the particle being
   126     represented, and optionally the mass (mass) and total decay width
   127     (totalwidth) of that particle in the SLHA scenario. Masses may also be found
   128     via the MASS block, from which the Particle.mass property is filled, if at
   129     all. They also store a list of Decay objects (decays) which are probably the
   130     item of most interest.
   131     """
   132     def __init__(self, pid, totalwidth=None, mass=None):
   133         self.pid = pid
   134         self.totalwidth = totalwidth
   135         self.mass = mass
   136         self.decays = []
   138     def add_decay(self, br, nda, ids):
   139         self.decays.append(Decay(br, nda, ids))
   140         self.decays.sort()
   142     def __cmp__(self, other):
   143         if abs(self.pid) == abs(other.pid):
   144             return cmp(self.pid, other.pid)
   145         return cmp(abs(self.pid), abs(other.pid))
   147     def __str__(self):
   148         s = str(self.pid)
   149         if self.mass is not None:
   150             s += " : mass = %e GeV" % self.mass
   151         if self.totalwidth is not None:
   152             s += " : total width = %e GeV" % self.totalwidth
   153         for d in self.decays:
   154             if d.br > 0.0:
   155                 s += "\n  %s" % d
   156         return s
   158     def __repr__(self):
   159         return self.__str__()
   163 def readSLHAFile(spcfilename, **kwargs):
   164     """
   165     Read an SLHA file, returning dictionaries of blocks and decays.
   167     Other keyword parameters are passed to readSLHA.
   168     """
   169     f = open(spcfilename, "r")
   170     rtn = readSLHA(f.read(), kwargs)
   171     f.close()
   172     return rtn
   175 def readSLHA(spcstr, ignorenobr=False):
   176     """
   177     Read an SLHA definition from a string, returning dictionaries of blocks and
   178     decays.
   180     If the ignorenobr parameter is True, do not store decay entries with a
   181     branching ratio of zero.
   182     """
   183     blocks = {}
   184     decays = {}
   185     #
   186     import re
   187     currentblock = None
   188     currentdecay = None
   189     for line in spcstr.splitlines():
   190         ## Handle (ignore) comment lines
   191         if line.startswith("#"):
   192             continue
   193         if "#" in line:
   194             line = line[:line.index("#")]
   196         ## Handle BLOCK/DECAY start lines
   197         if line.upper().startswith("BLOCK"):
   198             #print line
   199             match = re.match(r"BLOCK\s+(\w+)(\s+Q=\s*.+)?", line.upper())
   200             if not match:
   201                 continue
   202             blockname = match.group(1)
   203             qstr = match.group(2)
   204             if qstr is not None:
   205                 qstr = qstr[2:].strip()
   206             currentblock = blockname
   207             currentdecay = None
   208             blocks[blockname] = Block(blockname, q=qstr)
   209         elif line.upper().startswith("DECAY"):
   210             match = re.match(r"DECAY\s+(\d+)\s+([\d\.E+-]+).*", line.upper())
   211             if not match:
   212                 continue
   213             pdgid = int(match.group(1))
   214             width = float(match.group(2))
   215             currentblock = "DECAY"
   216             currentdecay = pdgid
   217             decays[pdgid] = Particle(pdgid, width)
   218         else:
   219             ## In-block line
   220             if currentblock is not None:
   221                 items = line.split()
   222                 if len(items) < 1:
   223                     continue
   224                 if currentblock != "DECAY":
   225                     if len(items) < 2:
   226                         ## Treat the ALPHA block differently
   227                         blocks[currentblock].value = _autotype(items[0])
   228                         blocks[currentblock].entries = _autotype(items[0])
   229                     else:
   230                         blocks[currentblock].add_entry(items)
   231                 else:
   232                     br = float(items[0])
   233                     nda = int(items[1])
   234                     ids = map(int, items[2:])
   235                     if br > 0.0 or not ignorenobr:
   236                         decays[currentdecay].add_decay(br, nda, ids)
   238     ## Try to populate Particle masses from the MASS block
   239     # print blocks.keys()
   240     try:
   241         for pid in blocks["MASS"].entries.keys():
   242             if decays.has_key(pid):
   243                 decays[pid].mass = blocks["MASS"].entries[pid]
   244     except:
   245         raise Exception("No MASS block found, from which to populate particle masses")
   247     return blocks, decays
   250 def readISAWIGFile(isafilename, **kwargs):
   251     """
   252     Read a spectrum definition from a file in the ISAWIG format, returning
   253     dictionaries of blocks and decays. While this is not an SLHA format, it is
   254     informally supported as a useful mechanism for converting ISAWIG spectra to
   255     SLHA.
   257     Other keyword parameters are passed to readSLHA.
   258     """
   259     f = open(isafilename, "r")
   260     rtn = readISAWIG(f.read(), kwargs)
   261     f.close()
   262     return rtn
   265 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
   266     """
   267     Write an SLHA file from the supplied blocks and decays dicts.
   269     Other keyword parameters are passed to writeSLHA.
   270     """
   271     f = open(spcfilename, "w")
   272     f.write(writeSLHA(blocks, decays, kwargs))
   273     f.close()
   276 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   279 def writeSLHA(blocks, decays, ignorenobr=False):
   280     """
   281     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   282     """
   283     sep = "   "
   284     out = ""
   285     def dict_hier_strs(d, s=""):
   286         if type(d) is dict:
   287             for k, v in sorted(d.iteritems()):
   288                 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
   289                     yield s2
   290         else:
   291             yield s + sep + _autostr(d)
   292     ## Blocks
   293     for bname, b in sorted(blocks.iteritems()):
   294         namestr = b.name
   295         if b.q is not None:
   296             namestr += " Q= %e" % b.q
   297         out += "BLOCK %s\n" % namestr
   298         for s in dict_hier_strs(b.entries):
   299             out += sep + s + "\n"
   300         out += "\n"
   301     ## Decays
   302     for pid, particle in sorted(decays.iteritems()):
   303         out += "DECAY %d %e\n" % (particle.pid, particle.totalwidth or -1)
   304         for d in sorted(particle.decays):
   305             if d.br > 0.0 or not ignorenobr:
   306                 products_str = "   ".join(map(str, d.ids))
   307                 out += sep + "%e" % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
   308         out += "\n"
   309     return out
   313 ###############################################################################
   314 ## ISAWIG handling
   316 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
   317 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   318 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   319 _HERWIGID2PDGID = {}
   320 _HERWIGID2PDGID[7]   = -1
   321 _HERWIGID2PDGID[8]   = -2
   322 _HERWIGID2PDGID[9]   = -3
   323 _HERWIGID2PDGID[10]  = -4
   324 _HERWIGID2PDGID[11]  = -5
   325 _HERWIGID2PDGID[12]  = -6
   326 _HERWIGID2PDGID[13]  =  21
   327 _HERWIGID2PDGID[59]  =  22
   328 _HERWIGID2PDGID[121] =  11
   329 _HERWIGID2PDGID[122] =  12
   330 _HERWIGID2PDGID[123] =  13
   331 _HERWIGID2PDGID[124] =  14
   332 _HERWIGID2PDGID[125] =  15
   333 _HERWIGID2PDGID[126] =  16
   334 _HERWIGID2PDGID[127] = -11
   335 _HERWIGID2PDGID[128] = -12
   336 _HERWIGID2PDGID[129] = -13
   337 _HERWIGID2PDGID[130] = -14
   338 _HERWIGID2PDGID[131] = -15
   339 _HERWIGID2PDGID[132] = -16
   340 _HERWIGID2PDGID[198] =  24 # W+
   341 _HERWIGID2PDGID[199] = -24 # W-
   342 _HERWIGID2PDGID[200] =  23 # Z0
   343 _HERWIGID2PDGID[201] =  25 ## SM HIGGS
   344 _HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   345 _HERWIGID2PDGID[204] =  35 ## HIGGSH0
   346 _HERWIGID2PDGID[205] =  36 ## HIGGSA0
   347 _HERWIGID2PDGID[206] =  37 ## HIGGS+
   348 _HERWIGID2PDGID[207] = -37 ## HIGGS-
   349 _HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   350 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   351 _HERWIGID2PDGID[402] =  1000002 ## SSULBR
   352 _HERWIGID2PDGID[408] = -1000002 ## SSUL
   353 _HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   354 _HERWIGID2PDGID[409] = -1000003 ## SSSL
   355 _HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   356 _HERWIGID2PDGID[410] = -1000004 ## SSCL
   357 _HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   358 _HERWIGID2PDGID[411] = -1000005 ## SSB1
   359 _HERWIGID2PDGID[406] =  1000006 ## SST1BR
   360 _HERWIGID2PDGID[412] = -1000006 ## SST1
   361 _HERWIGID2PDGID[413] =  2000001 ## SSDR
   362 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   363 _HERWIGID2PDGID[414] =  2000002 ## SSUR
   364 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
   365 _HERWIGID2PDGID[415] =  2000003 ## SSSR
   366 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   367 _HERWIGID2PDGID[416] =  2000004 ## SSCR
   368 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   369 _HERWIGID2PDGID[417] =  2000005 ## SSB2
   370 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   371 _HERWIGID2PDGID[418] =  2000006 ## SST2
   372 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
   373 _HERWIGID2PDGID[425] =  1000011 ## SSEL-
   374 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
   375 _HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   376 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   377 _HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   378 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   379 _HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   380 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   381 _HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   382 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   383 _HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   384 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   385 _HERWIGID2PDGID[437] =  2000011 ## SSEL-
   386 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
   387 _HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   388 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   389 _HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   390 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   391 _HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   392 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   393 _HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   394 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   395 _HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   396 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   397 _HERWIGID2PDGID[449] =  1000021 ## GLUINO
   398 _HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   399 _HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   400 _HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   401 _HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   402 _HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   403 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   404 _HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   405 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   406 _HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   408 def herwigid2pdgid(hwid):
   409     """
   410     Convert a particle ID code in the HERWIG internal IDHW format (as used by
   411     ISAWIG) into its equivalent in the standard PDG ID code definition.
   412     """
   413     return _HERWIGID2PDGID.get(hwid, hwid)
   416 def readISAWIG(isastr, ignorenobr=False):
   417     """
   418     Read a spectrum definition from a string in the ISAWIG format, returning
   419     dictionaries of blocks and decays. While this is not an SLHA format, it is
   420     informally supported as a useful mechanism for converting ISAWIG spectra to
   421     SLHA.
   423     ISAWIG parsing based on the HERWIG SUSY specification format, from
   424     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   426     If the ignorenobr parameter is True, do not store decay entries with a
   427     branching ratio of zero.
   428     """
   431     blocks = {}
   432     decays = {}
   433     LINES = isastr.splitlines()
   435     def getnextvalidline():
   436         while LINES:
   437             s = LINES.pop(0).strip()
   438             ## Return None if EOF reached
   439             if len(s) == 0:
   440                 continue
   441             ## Strip comments
   442             if "#" in s:
   443                 s = s[:s.index("#")].strip()
   444             ## Return if non-empty
   445             if len(s) > 0:
   446                 return s
   448     def getnextvalidlineitems():
   449         return map(_autotype, getnextvalidline().split())
   451     ## Populate MASS block and create decaying particle objects
   452     masses = Block("MASS")
   453     numentries = int(getnextvalidline())
   454     for i in xrange(numentries):
   455         hwid, mass, lifetime = getnextvalidlineitems()
   456         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   457         pdgid = herwigid2pdgid(hwid)
   458         masses.add_entry((pdgid, mass))
   459         decays[pdgid] = Particle(pdgid, width, mass)
   460         #print pdgid, mass, width
   461     blocks["MASS"] = masses
   463     ## Populate decays
   464     for n in xrange(numentries):
   465         numdecays = int(getnextvalidline())
   466         for d in xrange(numdecays):
   467             #print n, numentries-1, d, numdecays-1
   468             decayitems = getnextvalidlineitems()
   469             hwid = decayitems[0]
   470             pdgid = herwigid2pdgid(hwid)
   471             br = decayitems[1]
   472             nme = decayitems[2]
   473             daughter_hwids = decayitems[3:]
   474             daughter_pdgids = []
   475             for hw in daughter_hwids:
   476                 if hw != 0:
   477                     daughter_pdgids.append(herwigid2pdgid(hw))
   478             if not decays.has_key(pdgid):
   479                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   480                 decays[pdgid] = Particle(pdgid)
   481             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   484     ## Now the SUSY parameters
   485     TANB, ALPHAH = getnextvalidlineitems()
   486     blocks["MINPAR"] = Block("MINPAR")
   487     blocks["MINPAR"].add_entry((3, TANB))
   488     blocks["ALPHA"] = Block("ALPHA")
   489     blocks["ALPHA"].entries = ALPHAH
   490     #
   491     ## Neutralino mixing matrix
   492     blocks["NMIX"] = Block("NMIX")
   493     for i in xrange(1, 5):
   494         nmix_i = getnextvalidlineitems()
   495         for j, v in enumerate(nmix_i):
   496             blocks["NMIX"].add_entry((i, j+1, v))
   497     #
   498     ## Chargino mixing matrices V and U
   499     blocks["VMIX"] = Block("VMIX")
   500     vmix = getnextvalidlineitems()
   501     blocks["VMIX"].add_entry((1, 1, vmix[0]))
   502     blocks["VMIX"].add_entry((1, 2, vmix[1]))
   503     blocks["VMIX"].add_entry((2, 1, vmix[2]))
   504     blocks["VMIX"].add_entry((2, 2, vmix[3]))
   505     blocks["UMIX"] = Block("UMIX")
   506     umix = getnextvalidlineitems()
   507     blocks["UMIX"].add_entry((1, 1, umix[0]))
   508     blocks["UMIX"].add_entry((1, 2, umix[1]))
   509     blocks["UMIX"].add_entry((2, 1, umix[2]))
   510     blocks["UMIX"].add_entry((2, 2, umix[3]))
   511     #
   512     THETAT, THETAB, THETAL = getnextvalidlineitems()
   513     import math
   514     blocks["STOPMIX"] = Block("STOPMIX")
   515     blocks["STOPMIX"].add_entry((1, 1,  math.cos(THETAT)))
   516     blocks["STOPMIX"].add_entry((1, 2, -math.sin(THETAT)))
   517     blocks["STOPMIX"].add_entry((2, 1,  math.sin(THETAT)))
   518     blocks["STOPMIX"].add_entry((2, 2,  math.cos(THETAT)))
   519     blocks["SBOTMIX"] = Block("SBOTMIX")
   520     blocks["SBOTMIX"].add_entry((1, 1,  math.cos(THETAB)))
   521     blocks["SBOTMIX"].add_entry((1, 2, -math.sin(THETAB)))
   522     blocks["SBOTMIX"].add_entry((2, 1,  math.sin(THETAB)))
   523     blocks["SBOTMIX"].add_entry((2, 2,  math.cos(THETAB)))
   524     blocks["STAUMIX"] = Block("STAUMIX")
   525     blocks["STAUMIX"].add_entry((1, 1,  math.cos(THETAL)))
   526     blocks["STAUMIX"].add_entry((1, 2, -math.sin(THETAL)))
   527     blocks["STAUMIX"].add_entry((2, 1,  math.sin(THETAL)))
   528     blocks["STAUMIX"].add_entry((2, 2,  math.cos(THETAL)))
   529     #
   530     ATSS, ABSS, ALSS = getnextvalidlineitems()
   531     blocks["AU"] = Block("AU")
   532     blocks["AU"].add_entry((3, 3, ATSS))
   533     blocks["AD"] = Block("AD")
   534     blocks["AD"].add_entry((3, 3, ABSS))
   535     blocks["AE"] = Block("AE")
   536     blocks["AE"].add_entry((3, 3, ALSS))
   537     #
   538     MUSS = getnextvalidlineitems()[0]
   539     blocks["MINPAR"].add_entry((4, MUSS))
   540     #
   541     return blocks, decays
   544 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   545 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   546 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   547 _PDGID2HERWIGID = {}
   548 _PDGID2HERWIGID[      -1] = 7
   549 _PDGID2HERWIGID[      -2] = 8
   550 _PDGID2HERWIGID[      -3] = 9
   551 _PDGID2HERWIGID[      -4] = 10
   552 _PDGID2HERWIGID[      -5] = 11
   553 _PDGID2HERWIGID[      -6] = 12
   554 _PDGID2HERWIGID[      21] = 13
   555 _PDGID2HERWIGID[      22] = 59
   556 _PDGID2HERWIGID[      11] = 121
   557 _PDGID2HERWIGID[      12] = 122
   558 _PDGID2HERWIGID[      13] = 123
   559 _PDGID2HERWIGID[      14] = 124
   560 _PDGID2HERWIGID[      15] = 125
   561 _PDGID2HERWIGID[      16] = 126
   562 _PDGID2HERWIGID[     -11] = 127
   563 _PDGID2HERWIGID[     -12] = 128
   564 _PDGID2HERWIGID[     -13] = 129
   565 _PDGID2HERWIGID[     -14] = 130
   566 _PDGID2HERWIGID[     -15] = 131
   567 _PDGID2HERWIGID[     -16] = 132
   568 _PDGID2HERWIGID[      24] = 198 ## W+
   569 _PDGID2HERWIGID[     -24] = 199 ## W-
   570 _PDGID2HERWIGID[      23] = 200 ## Z
   571 _PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
   572 _PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   573 _PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   574 _PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   575 _PDGID2HERWIGID[      37] = 206 ## HIGGS+
   576 _PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   577 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   578 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   579 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   580 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
   581 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   582 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
   583 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   584 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
   585 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   586 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
   587 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   588 _PDGID2HERWIGID[-1000006] = 412 ## SST1
   589 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   590 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   591 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   592 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   593 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   594 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   595 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   596 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   597 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   598 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   599 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
   600 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   601 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   602 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   603 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   604 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   605 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   606 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   607 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   608 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   609 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   610 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   611 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   612 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   613 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   614 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   615 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   616 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   617 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   618 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   619 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   620 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   621 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   622 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   623 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   624 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   625 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   626 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   627 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   628 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   629 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   630 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   631 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   632 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   633 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   634 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   636 def pdgid2herwigid(pdgid):
   637     """
   638     Convert a particle ID code in the standard PDG ID code definition into
   639     its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
   640     """
   641     return _PDGID2HERWIGID.get(pdgid, pdgid)
   646 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
   647     """
   648     Write an ISAWIG file from the supplied blocks and decays dicts.
   650     Other keyword parameters are passed to writeISAWIG.
   652     TODO: Handle RPV SUSY
   653     """
   654     f = open(isafilename, "w")
   655     f.write(writeISAWIG(blocks, decays, kwargs))
   656     f.close()
   659 def writeISAWIG(blocks, decays, ignorenobr=False):
   660     """
   661     Return an ISAWIG definition as a string, from the supplied blocks and decays dicts.
   663     ISAWIG parsing based on the HERWIG SUSY specification format, from
   664     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   666     If the ignorenobr parameter is True, do not write decay entries with a
   667     branching ratio of zero.
   668     """
   670     masses = blocks["MASS"].entries
   672     ## Init output string
   673     out = ""
   675     ## First write out masses section:
   676     ##   Number of SUSY + top particles
   677     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   678     ##   repeated for each particle
   679     ## IDHW is the HERWIG identity code.
   680     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   681     massout = ""
   682     for pid in masses.keys():
   683         lifetime = -1
   684         try:
   685             width = decays[pid].totalwidth
   686             if width and width > 0:
   687                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   688         except:
   689             pass
   690         massout += "%d %e %e\n" % (pdgid2herwigid(pid), masses[pid], lifetime)
   691     out += "%d\n" % massout.count("\n")
   692     out += massout
   694     assert(len(masses) == len(decays))
   696     ## Next each particles decay modes together with their branching ratios and matrix element codes
   697     ##   Number of decay modes for a given particle (IDK)
   698     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   699     ##     repeated for each mode.
   700     ##   Repeated for each particle.
   701     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   702     ## the decay mode. NME is a code for the matrix element to be used, either from the
   703     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   704     for i, pid in enumerate(decays.keys()):
   705         # if not decays.has_key(pid):
   706         #     continue
   707         hwid = pdgid2herwigid(pid)
   708         decayout = ""
   709         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   710         for i_d, d in enumerate(decays[pid].decays):
   711             ## Skip decay if it has no branching ratio
   712             if ignorenobr and d.br == 0:
   713                 continue
   715             ## Identify decay matrix element to use
   716             ## From std HW docs, or from this pair:
   717             ## Two new matrix element codes have been added for these new decays:
   718             ##    NME =	200 	3 body top quark via charged Higgs
   719             ##    	300 	3 body R-parity violating gaugino and gluino decays
   720             nme = 0
   721             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   722             if abs(pid) in (6, 12):
   723                 nme = 100
   724             ## Extra SUSY MEs
   725             if len(d.ids) == 3:
   726                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   727                 pass
   728             decayout += "%d %e %d " % (hwid, d.br, nme)
   730             def is_quark(pid):
   731                 return (abs(pid) in range(1, 7))
   733             def is_lepton(pid):
   734                 return (abs(pid) in range(11, 17))
   736             def is_squark(pid):
   737                 if abs(pid) in range(1000001, 1000007):
   738                     return True
   739                 if abs(pid) in range(2000001, 2000007):
   740                     return True
   741                 return False
   743             def is_slepton(pid):
   744                 if abs(pid) in range(1000011, 1000017):
   745                     return True
   746                 if abs(pid) in range(2000011, 2000016, 2):
   747                     return True
   748                 return False
   750             def is_gaugino(pid):
   751                 if abs(pid) in range(1000022, 1000026):
   752                     return True
   753                 if abs(pid) in (1000035, 1000037):
   754                     return True
   755                 return False
   757             def is_susy(pid):
   758                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   760             absids = map(abs, d.ids)
   762             ## Order decay products as required by HERWIG
   763             ## Top
   764             if abs(pid) == 6:
   765                 def cmp_bottomlast(a, b):
   766                     """Comparison function which always puts b/bbar last"""
   767                     if abs(a) == 5:
   768                         return True
   769                     if abs(b) == 5:
   770                         return False
   771                     return cmp(a, b)
   772                 if len(absids) == 2:
   773                     ## 2 body mode, to Higgs: Higgs; Bottom
   774                     if (25 in absids or 26 in absids) and 5 in absids:
   775                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   776                 elif len(absids) == 3:
   777                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   778                     if 37 in absids or 23 in absids:
   779                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   780             ## Gluino
   781             elif abs(pid) == 1000021:
   782                 if len(absids) == 2:
   783                     ## 2 body mode
   784                     ## without gluon: any order
   785                     ## with gluon: gluon; colour neutral
   786                     if 21 in absids:
   787                         def cmp_gluonfirst(a, b):
   788                             """Comparison function which always puts gluon first"""
   789                             if a == 21:
   790                                 return False
   791                             if b == 21:
   792                                 return True
   793                             return cmp(a, b)
   794                         d.ids = sorted(d.ids, key=cmp_gluonfirst)
   795                 elif len(absids) == 3:
   796                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   797                     def cmp_quarkslast(a, b):
   798                         """Comparison function which always puts quarks last"""
   799                         if is_quark(a):
   800                             return True
   801                         if is_quark(b):
   802                             return False
   803                         return cmp(a, b)
   804                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   805             ## Squark/Slepton
   806             elif is_squark(pid) or is_slepton(pid):
   807                 def cmp_susy_quark_lepton(a, b):
   808                     if is_susy(a):
   809                         return False
   810                     if is_susy(b):
   811                         return True
   812                     if is_quark(a):
   813                         return False
   814                     if is_quark(b):
   815                         return True
   816                     return cmp(a, b)
   817                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   818                 ##                                                      Gluino       lepton
   819                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   820                 ## Squark
   821                 ##   2 body modes: Lepton Number Violated               quark     lepton
   822                 ##                 Baryon number violated               quark     quark
   823                 ## Slepton
   824                 ##   2 body modes: Lepton Number Violated               q or qbar
   825                 d.ids = sorted(d.ids, key=cmp_bottomlast)
   826             ## Higgs
   827             elif pid in (25, 26):
   828                 # TODO: Includes SUSY Higgses?
   829                 ## Higgs
   830                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   831                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   832                 ##   3 body modes:                                      colour neutral       q or qbar
   833                 if len(absids) == 3:
   834                     def cmp_quarkslast(a, b):
   835                         """Comparison function which always puts quarks last"""
   836                         if is_quark(a):
   837                             return True
   838                         if is_quark(b):
   839                             return False
   840                         return cmp(a, b)
   841                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   842             elif is_gaugino(pid):
   843                 # TODO: Is there actually anything to do here?
   844                 ## Gaugino
   845                 ##   2 body modes: Squark-quark                         q or sq
   846                 ##                 Slepton-lepton                       l or sl
   847                 ##
   848                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
   849                 ##                                                                           l or lbar
   850                 if len(absids) == 3:
   851                     def cmp_quarkslast(a, b):
   852                         """Comparison function which always puts quarks last"""
   853                         if is_quark(a):
   854                             return True
   855                         if is_quark(b):
   856                             return False
   857                         return cmp(a, b)
   858                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   860             # TODO: Gaugino/Gluino
   861             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
   863             ## Pad out IDs list with zeros
   864             ids = [0,0,0,0,0]
   865             for i, pid in enumerate(d.ids):
   866                 ids[i] = pid
   867             ids = map(str, ids)
   868             decayout += " ".join(ids) + "\n"
   869         decayout = "%d\n" % decayout.count("\n") + decayout
   870         out += decayout
   872     ## Now the SUSY parameters
   873     ## TANB, ALPHAH:
   874     out += "%e %e\n" % (blocks["MINPAR"].entries[3], blocks["ALPHA"].entries)
   875     ## Neutralino mixing matrix
   876     nmix = blocks["NMIX"].entries
   877     for i in xrange(1, 5):
   878         out += "%e %e %e %e\n" % (nmix[i][1], nmix[i][2], nmix[i][3], nmix[i][4])
   879     ## Chargino mixing matrices V and U
   880     vmix = blocks["VMIX"].entries
   881     out += "%e %e %e %e\n" % (vmix[1][1], vmix[1][2], vmix[2][1], vmix[2][2])
   882     umix = blocks["UMIX"].entries
   883     out += "%e %e %e %e\n" % (umix[1][1], umix[1][2], umix[2][1], umix[2][2])
   884     # THETAT,THETAB,THETAL
   885     import math
   886     out += "%e %e %e\n" % (math.acos(blocks["STOPMIX"].entries[1][1]),
   887                            math.acos(blocks["SBOTMIX"].entries[1][1]),
   888                            math.acos(blocks["STAUMIX"].entries[1][1]))
   889     # ATSS,ABSS,ALSS
   890     out += "%e %e %e\n" % (blocks["AU"].entries[3][3],
   891                            blocks["AD"].entries[3][3],
   892                            blocks["AE"].entries[3][3])
   893     # MUSS == sign(mu)
   894     out += "%f\n" % blocks["MINPAR"].entries[4]
   896     ## TODO: Handle RPV SUSY
   898     return out
   902 if __name__ == "__main__":
   903     import sys
   904     for a in sys.argv[1:]:
   905         if a.endswith(".isa"):
   906             blocks, decays = readISAWIGFile(a)
   907         else:
   908             blocks, decays = readSLHAFile(a)
   910         for bname, b in sorted(blocks.iteritems()):
   911             print b
   912             print
   914         print blocks.keys()
   916         print blocks["MASS"].entries[25]
   917         print
   919         for p in sorted(decays.values()):
   920             print p
   921             print
   923         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial