pyslha.py

Mon, 07 Mar 2011 10:48:37 +0000

author
Andy Buckley <andy@insectnation.org>
date
Mon, 07 Mar 2011 10:48:37 +0000
changeset 134
323754f1d261
parent 131
12a1338fa626
child 137
0aecd3e7f444
permissions
-rw-r--r--

Use more semantic handling of the format string specifier, and try to avoid using constructs not available in SLC5's native Python

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

mercurial