pyslha.py

Sat, 28 May 2011 23:20:53 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sat, 28 May 2011 23:20:53 +0100
changeset 145
718d8a3355f0
parent 144
c549f909067c
child 147
55fc50e4ac53
permissions
-rw-r--r--

Clarify that SLHA2 actually seems to be fine

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

mercurial