pyslha.py

Sun, 27 Feb 2011 00:14:05 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sun, 27 Feb 2011 00:14:05 +0100
changeset 123
c85e29bc13c4
parent 120
bffefe12df80
child 124
7bd52be093b2
permissions
-rw-r--r--

Enabling simultaneous rendering of multiple input files.

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

mercurial