pyslha.py

Sun, 28 Apr 2013 22:32:50 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 22:32:50 +0200
changeset 208
67e368e5f414
parent 207
2990422f1738
child 210
0f4f5472b7d8
permissions
-rw-r--r--

More docs

     1 #! /usr/bin/env python
     3 """\
     4 A simple but flexible handler of the SUSY Les Houches Accord (SLHA) data format.
     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 I'm aware is also
    12 fully compatible with SLHA2: the block structures are read and accessed
    13 generically. If you have any problems, please provide an example input file and
    14 I'll happily 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 Users of version 1.x should note that the interface has changed a little in
    21 version 2.0.0 and onward, in particular in the interface of the Block objects,
    22 which are now more dict-like: entries can be added and accessed via the usual
    23 square-brackets indexing operators, including for multiple indices as is common
    24 for mixing matrices e.g. NMIX[1,2] as opposed to the old NMIX.entries[1][2]
    25 way. This does break backward compatibility but is a big improvement both for
    26 internal code sanity and usability. The Block interface also now supplies
    27 dict-like has_key(), keys(), and items() methods, as well as more specialist
    28 value(), set_value() and is_single_valued() methods for improved access to ALPHA
    29 and any other unindexed blocks.
    31 Please cite PySLHA via the URL http://insectnation.org/projects/pyslha if it is
    32 useful to you in the production of an academic paper. A PySLHA publication will
    33 be provided at some point... watch this space!
    36 TODOs:
    38   For 2.0.0:
    39    * Output column alignment cosmetics
    40    * Precision setting obedience in SLHA output of values
    42   For 2.1.0:
    43    * Preserve comments from read -> write (needs full-line/inline comment separation?)
    44    * Split writeSLHA into writeSLHA{Blocks,Decays}
    46   Later/maybe:
    47    * Identify HERWIG decay matrix element to use in ISAWIG
    48    * Handle RPV SUSY in ISAWIG
    49 """
    51 __author__ = "Andy Buckley <andy.buckley@cern.ch"
    52 __version__ = "2.0.0a0"
    55 def _mkdict():
    56     """Try to return an empty ordered dict, but fall back to normal dict if necessary"""
    57     try:
    58         from collections import OrderedDict
    59         return OrderedDict()
    60     except:
    61         try:
    62             from ordereddict import OrderedDict
    63             return OrderedDict()
    64         except:
    65             return dict()
    67 def _autotype(var):
    68     """Automatically convert strings to numerical types if possible."""
    69     if type(var) is not str:
    70         return var
    71     if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
    72         return int(var)
    73     try:
    74         f = float(var)
    75         return f
    76     except ValueError:
    77         return var
    79 def _autostr(var, precision=8):
    80     """Automatically format numerical types as the right sort of string."""
    81     if type(var) is float:
    82         return ("%." + str(precision) + "e") % var
    83     return str(var)
    86 class AccessError(Exception):
    87     "Exception object to be raised when a SLHA block is accessed in an invalid way"
    88     def __init__(self, errmsg):
    89         self.msg = errmsg
    90     def __str__(self):
    91         return self.msg
    93 class ParseError(Exception):
    94     "Exception object to be raised when a spectrum file/string is malformed"
    95     def __init__(self, errmsg):
    96         self.msg = errmsg
    97     def __str__(self):
    98         return self.msg
   102 class Block(object):
   103     """
   104     Object representation of any BLOCK elements read from an SLHA file.
   106     Blocks have a name, may have an associated Q value, and contain a collection
   107     of data entries, each indexed by one or more keys. Entries in the dictionary
   108     are stored as numeric types (int or float) when a cast from the string in
   109     the file has been possible.
   111     Block is closely related to a Python dict (and, in fact, is implemented via
   112     an OrderedDict when possible). The preferred methods of entry access use the
   113     dict-like [] operator for getting and setting, and the keys() and items()
   114     methods for iteration. Purely iterating over the object behaves like keys(),
   115     as for an ordinary dict.
   117     Multiple indices are possible, especially for entries in mixing matrix
   118     blocks. These are now implemented in the natural way, e.g. for access to the
   119     (1,2) element of a mxing matrix block, use bmix[1,2] = 0.123 and print
   120     bmix[1,2]. The value() and set_value() functions behave similarly.
   122     It is possible, although not usual, to store unindexed values in a
   123     block. This is only supported when that entry is the only one in the block,
   124     and it is stored in the normal fashion but with None as the lookup key. The
   125     value() method may be used without a key argument to retrieve this value, if
   126     the is_single_valued() method returns True, and similarly the set_value()
   127     method may be used to set it if only one argument is supplied and the object
   128     is compatible.
   129     """
   130     def __init__(self, name, q=None, entries=None):
   131         self.name = name
   132         self.entries = _mkdict()
   133         if entries is not None:
   134             self.entries.update(entries)
   135         self.q = _autotype(q)
   137     def add_entry(self, entry):
   138         """Add an entry to the block from an iterable (i.e. list or tuple) or
   139         whitespace-separated string.
   141         Indexing will be determined automatically such that there is always a
   142         single-element value: multi-value or None indices may be constructed
   143         implicitly.
   144         """
   145         ## If the argument is a single string, split it and proceed
   146         if type(entry) is str:
   147             entry = entry.split()
   148         ## Check that the arg is an iterable
   149         if not hasattr(entry, "__iter__"):
   150             raise AccessError("Block entries must be iterable")
   151         ## Auto-convert the types in the list
   152         entry = map(_autotype, entry)
   153         ## Add the entry to the map, with appropriate indices
   154         if len(entry) == 1:
   155             self.entries[None] = entry[0]
   156         elif len(entry) == 2:
   157             self.entries[entry[0]] = entry[1]
   158         else:
   159             self.entries[tuple(entry[:-1])] = entry[-1]
   161     def is_single_valued(self):
   162         """Return true if there is only one entry, and it has no index: the
   163         'value()' attribute may be used in that case without an argument."""
   164         return len(self.entries) == 1 and self.entries.keys()[0] is None
   166     def value(self, key=None):
   167         """Get a value from the block with the supplied key.
   169         If no key is given, then the block must contain only one non-indexed
   170         value otherwise an AccessError exception will be raised.\
   171         """
   172         if key == None and not self.is_single_valued():
   173             raise AccessError("Tried to access unique value of multi-value block")
   174         return self.entries[key]
   176     def set_value(self, *args):
   177         """Set a value in the block with supplied key, val args (or just val for
   178         an unindexed block).
   180         If no key is given, then the block must contain only one non-indexed
   181         value otherwise an AccessError exception will be raised.\
   182         """
   183         if len(args) == 0:
   184             raise AccessError("set_value() called without arguments")
   185         elif len(args) == 1:
   186             if len(self.entries) > 0 and not self.is_single_valued():
   187                 raise AccessError("Tried to set a unique value on a multi-value block")
   188             self.entries[None] = args[0]
   189         else:
   190             self.entries[tuple(args[:-1])] = args[-1]
   192     def has_key(self, key):
   193         """Does the block have the given key?"""
   194         return self.entries.has_key(key)
   196     def keys(self):
   197         """Access the block item keys."""
   198         return self.entries.keys()
   200     def items(self, key=None):
   201         """Access the block items as (key,value) tuples.
   203         Note: The Python 3 dict attribute 'items()' is used rather than the
   204         'old' Python 2 'iteritems()' name for forward-looking compatibility.\
   205         """
   206         return self.entries.iteritems()
   208     def __len__(self):
   209         return len(self.entries)
   211     def __iter(self):
   212         return self.entries.__iter__()
   214     def __getitem__(self, key):
   215         return self.entries[key]
   217     def __setitem__(self, key, value):
   218         self.entries[key] = value
   220     def __cmp__(self, other):
   221         return cmp(self.name, other.name) and cmp(self.entries, other.entries)
   223     def __repr__(self):
   224         s = self.name
   225         if self.q is not None:
   226             s += " (Q=%s)" % self.q
   227         def _format_kv(k, v):
   228             if type(k) is not tuple:
   229                 s = "%r" % k
   230             else:
   231                 s = ",".join("%r" % subindex for subindex in k)
   232             s += " : %r" % v
   233             return s
   234         s += " { " + "; ".join(_format_kv(k, v) for k, v in self.items()) + " }"
   235         return s
   238 class Decay(object):
   239     """
   240     Object representing a decay entry on a particle decribed by the SLHA file.
   241     'Decay' objects are not a direct representation of a DECAY block in an SLHA
   242     file... that role, somewhat confusingly, is taken by the Particle class.
   244     Decay objects have three properties: a branching ratio, br, an nda number
   245     (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
   246     decay occurs. The PDG ID of the particle whose decay this represents may
   247     also be stored, but this is normally known via the Particle in which the
   248     decay is stored.
   249     """
   250     def __init__(self, br, nda, ids, parentid=None):
   251         self.parentid = parentid
   252         self.br = br
   253         self.nda = nda
   254         self.ids = ids
   255         assert(self.nda == len(self.ids))
   257     def __cmp__(self, other):
   258         return cmp(other.br, self.br)
   260     def __str__(self):
   261         return "%.8e %s" % (self.br, self.ids)
   263     def __repr__(self):
   264         return self.__str__()
   267 class Particle(object):
   268     """
   269     Representation of a single, specific particle, decay block from an SLHA
   270     file.  These objects are not themselves called 'Decay', since that concept
   271     applies more naturally to the various decays found inside this
   272     object. Particle classes store the PDG ID (pid) of the particle being
   273     represented, and optionally the mass (mass) and total decay width
   274     (totalwidth) of that particle in the SLHA scenario. Masses may also be found
   275     via the MASS block, from which the Particle.mass property is filled, if at
   276     all. They also store a list of Decay objects (decays) which are probably the
   277     item of most interest.
   278     """
   279     def __init__(self, pid, totalwidth=None, mass=None):
   280         self.pid = pid
   281         self.totalwidth = totalwidth
   282         self.mass = mass
   283         self.decays = []
   285     def add_decay(self, br, nda, ids):
   286         self.decays.append(Decay(br, nda, ids))
   287         self.decays.sort()
   289     def __cmp__(self, other):
   290         if abs(self.pid) == abs(other.pid):
   291             return cmp(self.pid, other.pid)
   292         return cmp(abs(self.pid), abs(other.pid))
   294     def __str__(self):
   295         s = str(self.pid)
   296         if self.mass is not None:
   297             s += " : mass = %.8e GeV" % self.mass
   298         if self.totalwidth is not None:
   299             s += " : total width = %.8e GeV" % self.totalwidth
   300         for d in self.decays:
   301             if d.br > 0.0:
   302                 s += "\n  %s" % d
   303         return s
   305     def __repr__(self):
   306         return self.__str__()
   312 def readSLHA(spcstr, ignorenobr=False):
   313     """
   314     Read an SLHA definition from a string, returning dictionaries of blocks and
   315     decays.
   317     If the ignorenobr parameter is True, do not store decay entries with a
   318     branching ratio of zero.
   319     """
   320     blocks = _mkdict()
   321     decays = _mkdict()
   322     #
   323     import re
   324     currentblock = None
   325     currentdecay = None
   326     for line in spcstr.splitlines():
   327         ## Handle (ignore) comment lines
   328         # TODO: Store block/entry comments
   329         if line.startswith("#"):
   330             continue
   331         if "#" in line:
   332             line = line[:line.index("#")]
   334         ## Handle BLOCK/DECAY start lines
   335         if line.upper().startswith("BLOCK"):
   336             #print line
   337             match = re.match(r"BLOCK\s+(\w+)(\s+Q\s*=\s*.+)?", line.upper())
   338             if not match:
   339                 continue
   340             blockname = match.group(1)
   341             qstr = match.group(2)
   342             if qstr is not None:
   343                 qstr = qstr[qstr.find("=")+1:].strip()
   344             currentblock = blockname
   345             currentdecay = None
   346             blocks[blockname] = Block(blockname, q=qstr)
   347         elif line.upper().startswith("DECAY"):
   348             match = re.match(r"DECAY\s+(-?\d+)\s+([\d\.E+-]+|NAN).*", line.upper())
   349             if not match:
   350                 continue
   351             pdgid = int(match.group(1))
   352             width = float(match.group(2)) if match.group(2) != "NAN" else None
   353             currentblock = "DECAY"
   354             currentdecay = pdgid
   355             decays[pdgid] = Particle(pdgid, width)
   356         else:
   357             ## In-block line
   358             if currentblock is not None:
   359                 items = line.split()
   360                 if len(items) < 1:
   361                     continue
   362                 if currentblock != "DECAY":
   363                     blocks[currentblock].add_entry(items)
   364                 else:
   365                     br = float(items[0]) if items[0].upper() != "NAN" else None
   366                     nda = int(items[1])
   367                     ids = map(int, items[2:])
   368                     if br > 0.0 or not ignorenobr: # br == None is < 0
   369                         decays[currentdecay].add_decay(br, nda, ids)
   371     ## Try to populate Particle masses from the MASS block
   372     # print blocks.keys()
   373     try:
   374         for pid in blocks["MASS"].keys():
   375             if decays.has_key(pid):
   376                 decays[pid].mass = blocks["MASS"][pid]
   377     except:
   378         raise ParseError("No MASS block found: cannot populate particle masses")
   380     return blocks, decays
   385 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   388 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
   389     """
   390     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   391     """
   392     # TODO: Pay attention to space-padding and minus signs for column alignment
   393     fmte = "%." + str(precision) + "e"
   394     sep = "   "
   395     blockstrs = []
   396     ## Blocks
   397     for bname, b in blocks.iteritems():
   398         namestr = b.name
   399         if b.q is not None:
   400             namestr += (" Q= " + fmte) % float(b.q)
   401         blockstr = "BLOCK %s\n" % namestr
   402         entrystrs = []
   403         for k, v in b.items():
   404             entrystr = ""
   405             if type(k) == tuple:
   406                 entrystr += sep.join(_autostr(i) for i in k)
   407             elif k is not None:
   408                 entrystr += _autostr(k)
   409             entrystr += sep + _autostr(v) # TODO: apply precision formatting for floats
   410             entrystrs.append(entrystr)
   411         blockstr += "\n".join(entrystrs)
   412         blockstrs.append(blockstr)
   413         ##
   414     ## Decays
   415     for pid, particle in decays.iteritems():
   416         blockstr = ("DECAY %d " + fmte + "\n") % (particle.pid, particle.totalwidth or -1)
   417         decaystrs = []
   418         for d in particle.decays:
   419             if d.br > 0.0 or not ignorenobr:
   420                 products_str = sep.join(map(str, d.ids))
   421                 decaystr = sep + (fmte % d.br) + sep + ("%d" % len(d.ids)) + sep + products_str
   422                 decaystrs.append(decaystr)
   423         blockstr += "\n".join(decaystrs)
   424         blockstrs.append(blockstr)
   425     ## Total result
   426     return "\n\n".join(blockstrs)
   430 ###############################################################################
   431 ## PDG <-> HERWIG particle ID code translations for ISAWIG handling
   433 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
   434 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   435 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   436 _HERWIGID2PDGID = {}
   437 _HERWIGID2PDGID[7]   = -1
   438 _HERWIGID2PDGID[8]   = -2
   439 _HERWIGID2PDGID[9]   = -3
   440 _HERWIGID2PDGID[10]  = -4
   441 _HERWIGID2PDGID[11]  = -5
   442 _HERWIGID2PDGID[12]  = -6
   443 _HERWIGID2PDGID[13]  =  21
   444 _HERWIGID2PDGID[59]  =  22
   445 _HERWIGID2PDGID[121] =  11
   446 _HERWIGID2PDGID[122] =  12
   447 _HERWIGID2PDGID[123] =  13
   448 _HERWIGID2PDGID[124] =  14
   449 _HERWIGID2PDGID[125] =  15
   450 _HERWIGID2PDGID[126] =  16
   451 _HERWIGID2PDGID[127] = -11
   452 _HERWIGID2PDGID[128] = -12
   453 _HERWIGID2PDGID[129] = -13
   454 _HERWIGID2PDGID[130] = -14
   455 _HERWIGID2PDGID[131] = -15
   456 _HERWIGID2PDGID[132] = -16
   457 _HERWIGID2PDGID[198] =  24 # W+
   458 _HERWIGID2PDGID[199] = -24 # W-
   459 _HERWIGID2PDGID[200] =  23 # Z0
   460 _HERWIGID2PDGID[201] =  25 ## SM HIGGS
   461 _HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   462 _HERWIGID2PDGID[204] =  35 ## HIGGSH0
   463 _HERWIGID2PDGID[205] =  36 ## HIGGSA0
   464 _HERWIGID2PDGID[206] =  37 ## HIGGS+
   465 _HERWIGID2PDGID[207] = -37 ## HIGGS-
   466 _HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   467 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   468 _HERWIGID2PDGID[402] =  1000002 ## SSULBR
   469 _HERWIGID2PDGID[408] = -1000002 ## SSUL
   470 _HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   471 _HERWIGID2PDGID[409] = -1000003 ## SSSL
   472 _HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   473 _HERWIGID2PDGID[410] = -1000004 ## SSCL
   474 _HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   475 _HERWIGID2PDGID[411] = -1000005 ## SSB1
   476 _HERWIGID2PDGID[406] =  1000006 ## SST1BR
   477 _HERWIGID2PDGID[412] = -1000006 ## SST1
   478 _HERWIGID2PDGID[413] =  2000001 ## SSDR
   479 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   480 _HERWIGID2PDGID[414] =  2000002 ## SSUR
   481 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
   482 _HERWIGID2PDGID[415] =  2000003 ## SSSR
   483 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   484 _HERWIGID2PDGID[416] =  2000004 ## SSCR
   485 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   486 _HERWIGID2PDGID[417] =  2000005 ## SSB2
   487 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   488 _HERWIGID2PDGID[418] =  2000006 ## SST2
   489 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
   490 _HERWIGID2PDGID[425] =  1000011 ## SSEL-
   491 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
   492 _HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   493 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   494 _HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   495 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   496 _HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   497 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   498 _HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   499 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   500 _HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   501 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   502 _HERWIGID2PDGID[437] =  2000011 ## SSEL-
   503 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
   504 _HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   505 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   506 _HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   507 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   508 _HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   509 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   510 _HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   511 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   512 _HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   513 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   514 _HERWIGID2PDGID[449] =  1000021 ## GLUINO
   515 _HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   516 _HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   517 _HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   518 _HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   519 _HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   520 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   521 _HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   522 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   523 _HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   525 def herwigid2pdgid(hwid):
   526     """
   527     Convert a particle ID code in the HERWIG internal IDHW format (as used by
   528     ISAWIG) into its equivalent in the standard PDG ID code definition.
   529     """
   530     return _HERWIGID2PDGID.get(hwid, hwid)
   533 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   534 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   535 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   536 _PDGID2HERWIGID = {}
   537 _PDGID2HERWIGID[      -1] = 7
   538 _PDGID2HERWIGID[      -2] = 8
   539 _PDGID2HERWIGID[      -3] = 9
   540 _PDGID2HERWIGID[      -4] = 10
   541 _PDGID2HERWIGID[      -5] = 11
   542 _PDGID2HERWIGID[      -6] = 12
   543 _PDGID2HERWIGID[      21] = 13
   544 _PDGID2HERWIGID[      22] = 59
   545 _PDGID2HERWIGID[      11] = 121
   546 _PDGID2HERWIGID[      12] = 122
   547 _PDGID2HERWIGID[      13] = 123
   548 _PDGID2HERWIGID[      14] = 124
   549 _PDGID2HERWIGID[      15] = 125
   550 _PDGID2HERWIGID[      16] = 126
   551 _PDGID2HERWIGID[     -11] = 127
   552 _PDGID2HERWIGID[     -12] = 128
   553 _PDGID2HERWIGID[     -13] = 129
   554 _PDGID2HERWIGID[     -14] = 130
   555 _PDGID2HERWIGID[     -15] = 131
   556 _PDGID2HERWIGID[     -16] = 132
   557 _PDGID2HERWIGID[      24] = 198 ## W+
   558 _PDGID2HERWIGID[     -24] = 199 ## W-
   559 _PDGID2HERWIGID[      23] = 200 ## Z
   560 _PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
   561 _PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   562 _PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   563 _PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   564 _PDGID2HERWIGID[      37] = 206 ## HIGGS+
   565 _PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   566 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   567 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   568 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   569 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
   570 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   571 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
   572 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   573 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
   574 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   575 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
   576 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   577 _PDGID2HERWIGID[-1000006] = 412 ## SST1
   578 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   579 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   580 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   581 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   582 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   583 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   584 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   585 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   586 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   587 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   588 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
   589 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   590 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   591 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   592 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   593 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   594 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   595 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   596 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   597 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   598 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   599 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   600 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   601 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   602 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   603 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   604 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   605 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   606 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   607 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   608 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   609 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   610 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   611 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   612 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   613 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   614 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   615 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   616 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   617 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   618 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   619 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   620 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   621 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   622 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   623 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   625 def pdgid2herwigid(pdgid):
   626     """
   627     Convert a particle ID code in the standard PDG ID code definition into
   628     its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
   629     """
   630     return _PDGID2HERWIGID.get(pdgid, pdgid)
   633 ###############################################################################
   634 ## ISAWIG format reading/writing
   637 def readISAWIG(isastr, ignorenobr=False):
   638     """
   639     Read a spectrum definition from a string in the ISAWIG format, returning
   640     dictionaries of blocks and decays. While this is not an SLHA format, it is
   641     informally supported as a useful mechanism for converting ISAWIG spectra to
   642     SLHA.
   644     ISAWIG parsing based on the HERWIG SUSY specification format, from
   645     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   647     If the ignorenobr parameter is True, do not store decay entries with a
   648     branching ratio of zero.
   649     """
   651     blocks = _mkdict()
   652     decays = _mkdict()
   653     LINES = isastr.splitlines()
   655     def getnextvalidline():
   656         while LINES:
   657             s = LINES.pop(0).strip()
   658             # print "*", s, "*"
   659             ## Return None if EOF reached
   660             if len(s) == 0:
   661                 continue
   662             ## Strip comments
   663             if "#" in s:
   664                 s = s[:s.index("#")].strip()
   665             ## Return if non-empty
   666             if len(s) > 0:
   667                 return s
   669     def getnextvalidlineitems():
   670         return map(_autotype, getnextvalidline().split())
   672     ## Populate MASS block and create decaying particle objects
   673     masses = Block("MASS")
   674     numentries = int(getnextvalidline())
   675     for i in xrange(numentries):
   676         hwid, mass, lifetime = getnextvalidlineitems()
   677         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   678         pdgid = herwigid2pdgid(hwid)
   679         masses[pdgid] = mass
   680         decays[pdgid] = Particle(pdgid, width, mass)
   681         #print pdgid, mass, width
   682     blocks["MASS"] = masses
   684     ## Populate decays
   685     for n in xrange(numentries):
   686         numdecays = int(getnextvalidline())
   687         for d in xrange(numdecays):
   688             #print n, numentries-1, d, numdecays-1
   689             decayitems = getnextvalidlineitems()
   690             hwid = decayitems[0]
   691             pdgid = herwigid2pdgid(hwid)
   692             br = decayitems[1]
   693             nme = decayitems[2]
   694             daughter_hwids = decayitems[3:]
   695             daughter_pdgids = []
   696             for hw in daughter_hwids:
   697                 if hw != 0:
   698                     daughter_pdgids.append(herwigid2pdgid(hw))
   699             if not decays.has_key(pdgid):
   700                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   701                 decays[pdgid] = Particle(pdgid)
   702             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   705     ## Now the SUSY parameters
   706     TANB, ALPHAH = getnextvalidlineitems()
   707     blocks["MINPAR"] = Block("MINPAR")
   708     blocks["MINPAR"][3] = TANB
   709     blocks["ALPHA"] = Block("ALPHA")
   710     blocks["ALPHA"][None] = ALPHAH
   711     #
   712     ## Neutralino mixing matrix
   713     blocks["NMIX"] = Block("NMIX")
   714     for i in xrange(1, 5):
   715         nmix_i = getnextvalidlineitems()
   716         for j, v in enumerate(nmix_i):
   717             blocks["NMIX"][i, j+1] = v
   718     #
   719     ## Chargino mixing matrices V and U
   720     blocks["VMIX"] = Block("VMIX")
   721     vmix = getnextvalidlineitems()
   722     blocks["VMIX"][1, 1] = vmix[0]
   723     blocks["VMIX"][1, 2] = vmix[1]
   724     blocks["VMIX"][2, 1] = vmix[2]
   725     blocks["VMIX"][2, 2] = vmix[3]
   726     blocks["UMIX"] = Block("UMIX")
   727     umix = getnextvalidlineitems()
   728     blocks["UMIX"][1, 1] = umix[0]
   729     blocks["UMIX"][1, 2] = umix[1]
   730     blocks["UMIX"][2, 1] = umix[2]
   731     blocks["UMIX"][2, 2] = umix[3]
   732     #
   733     THETAT, THETAB, THETAL = getnextvalidlineitems()
   734     import math
   735     blocks["STOPMIX"] = Block("STOPMIX")
   736     blocks["STOPMIX"][1, 1] =  math.cos(THETAT)
   737     blocks["STOPMIX"][1, 2] = -math.sin(THETAT)
   738     blocks["STOPMIX"][2, 1] =  math.sin(THETAT)
   739     blocks["STOPMIX"][2, 2] =  math.cos(THETAT)
   740     blocks["SBOTMIX"] = Block("SBOTMIX")
   741     blocks["SBOTMIX"][1, 1] =  math.cos(THETAB)
   742     blocks["SBOTMIX"][1, 2] = -math.sin(THETAB)
   743     blocks["SBOTMIX"][2, 1] =  math.sin(THETAB)
   744     blocks["SBOTMIX"][2, 2] =  math.cos(THETAB)
   745     blocks["STAUMIX"] = Block("STAUMIX")
   746     blocks["STAUMIX"][1, 1] =  math.cos(THETAL)
   747     blocks["STAUMIX"][1, 2] = -math.sin(THETAL)
   748     blocks["STAUMIX"][2, 1] =  math.sin(THETAL)
   749     blocks["STAUMIX"][2, 2] =  math.cos(THETAL)
   750     #
   751     ATSS, ABSS, ALSS = getnextvalidlineitems()
   752     blocks["AU"] = Block("AU")
   753     blocks["AU"][3, 3] = ATSS
   754     blocks["AD"] = Block("AD")
   755     blocks["AD"][3, 3] = ABSS
   756     blocks["AE"] = Block("AE")
   757     blocks["AE"][3, 3] = ALSS
   758     #
   759     MUSS = getnextvalidlineitems()[0]
   760     blocks["MINPAR"][4] = MUSS
   761     #
   763     # TODO: Parse RPV boolean and couplings into SLHA2 blocks
   765     return blocks, decays
   768 def writeISAJET(blocks, decays, outname, ignorenobr=False, precision=8):
   769     """
   770     Return a SUSY spectrum definition in the format required for input by ISAJET,
   771     as a string, from the supplied blocks and decays dicts.
   773     The outname parameter specifies the desired output filename from ISAJET: this
   774     will appear in the first line of the return value.
   776     If the ignorenobr parameter is True, do not write decay entries with a
   777     branching ratio of zero.
   778     """
   779     fmte = "%." + str(precision) + "e"
   781     masses = blocks["MASS"]
   783     ## Init output string
   784     out = ""
   786     ## First line is the output name
   787     out += "'%s'" % outname + "\n"
   789     ## Next the top mass
   790     out += fmte % masses[6] + "\n"
   792     ## Next the top mass
   793     out += fmte % masses[6] + "\n"
   795     ## mSUGRA parameters (one line)
   796     # e.g. 1273.78,713.286,804.721,4.82337
   798     ## Masses and trilinear couplings (3 lines)
   799     # e.g. 1163.14,1114.15,1118.99,374.664,209.593
   800     # e.g. 1069.54,1112.7,919.908,374.556,209.381,-972.817,-326.745,-406.494
   801     # e.g. 1163.14,1114.15,1118.99,374.712,210.328
   803     ## RPV couplings (?? lines)
   804     # e.g. 232.615,445.477
   806     ## Etc ???!!!
   807     # e.g. /
   808     # e.g. n
   809     # e.g. y
   810     # e.g. y
   811     # e.g. 0.047441 3.80202e-23 0 0 0 2.17356e-22 0 0 5.23773e-09
   812     # e.g. y
   813     # e.g. 3.35297e-25 0 0 0 7.34125e-24 0 0 0 3.17951e-22 8.07984e-12 0 0 0 1.76906e-10 0 0 0 7.66184e-09 0 0 0 0 0 0 0 0 0
   814     # e.g. n
   815     # e.g. 'susy_RPV_stau_BC1scan_m560_tanb05.txt'
   817     return out
   820 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
   821     """
   822     Return a SUSY spectrum definition in the format produced by ISAWIG for inut to HERWIG
   823     as a string, from the supplied SLHA blocks and decays dicts.
   825     ISAWIG parsing based on the HERWIG SUSY specification format, from
   826     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   828     If the ignorenobr parameter is True, do not write decay entries with a
   829     branching ratio of zero.
   830     """
   831     fmte = "%." + str(precision) + "e"
   833     masses = blocks["MASS"]
   835     ## Init output string
   836     out = ""
   838     ## First write out masses section:
   839     ##   Number of SUSY + top particles
   840     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   841     ##   repeated for each particle
   842     ## IDHW is the HERWIG identity code.
   843     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   844     massout = ""
   845     for pid in masses.keys():
   846         lifetime = -1
   847         try:
   848             width = decays[pid].totalwidth
   849             if width and width > 0:
   850                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   851         except:
   852             pass
   853         massout += ("%d " + fmte + " " + fmte + "\n") % (pdgid2herwigid(pid), masses[pid], lifetime)
   854     out += "%d\n" % massout.count("\n")
   855     out += massout
   857     assert(len(masses) == len(decays))
   859     ## Next each particles decay modes together with their branching ratios and matrix element codes
   860     ##   Number of decay modes for a given particle (IDK)
   861     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   862     ##     repeated for each mode.
   863     ##   Repeated for each particle.
   864     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   865     ## the decay mode. NME is a code for the matrix element to be used, either from the
   866     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   867     for i, pid in enumerate(decays.keys()):
   868         # if not decays.has_key(pid):
   869         #     continue
   870         hwid = pdgid2herwigid(pid)
   871         decayout = ""
   872         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   873         for i_d, d in enumerate(decays[pid].decays):
   874             ## Skip decay if it has no branching ratio
   875             if ignorenobr and d.br == 0:
   876                 continue
   878             ## Identify decay matrix element to use
   879             ## From std HW docs, or from this pair:
   880             ## Two new matrix element codes have been added for these new decays:
   881             ##    NME =	200 	3 body top quark via charged Higgs
   882             ##    	300 	3 body R-parity violating gaugino and gluino decays
   883             nme = 0
   884             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   885             if abs(pid) in (6, 12):
   886                 nme = 100
   887             ## Extra SUSY MEs
   888             if len(d.ids) == 3:
   889                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   890                 pass
   891             decayout += ("%d " + fmte + " %d ") % (hwid, d.br, nme)
   893             def is_quark(pid):
   894                 return (abs(pid) in range(1, 7))
   896             def is_lepton(pid):
   897                 return (abs(pid) in range(11, 17))
   899             def is_squark(pid):
   900                 if abs(pid) in range(1000001, 1000007):
   901                     return True
   902                 if abs(pid) in range(2000001, 2000007):
   903                     return True
   904                 return False
   906             def is_slepton(pid):
   907                 if abs(pid) in range(1000011, 1000017):
   908                     return True
   909                 if abs(pid) in range(2000011, 2000016, 2):
   910                     return True
   911                 return False
   913             def is_gaugino(pid):
   914                 if abs(pid) in range(1000022, 1000026):
   915                     return True
   916                 if abs(pid) in (1000035, 1000037):
   917                     return True
   918                 return False
   920             def is_susy(pid):
   921                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   923             absids = map(abs, d.ids)
   925             ## Order decay products as required by HERWIG
   926             ## Top
   927             if abs(pid) == 6:
   928                 def cmp_bottomlast(a, b):
   929                     """Comparison function which always puts b/bbar last"""
   930                     if abs(a) == 5:
   931                         return True
   932                     if abs(b) == 5:
   933                         return False
   934                     return cmp(a, b)
   935                 if len(absids) == 2:
   936                     ## 2 body mode, to Higgs: Higgs; Bottom
   937                     if (25 in absids or 26 in absids) and 5 in absids:
   938                         d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   939                 elif len(absids) == 3:
   940                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   941                     if 37 in absids or 23 in absids:
   942                         d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   943             ## Gluino
   944             elif abs(pid) == 1000021:
   945                 if len(absids) == 2:
   946                     ## 2 body mode
   947                     ## without gluon: any order
   948                     ## with gluon: gluon; colour neutral
   949                     if 21 in absids:
   950                         def cmp_gluonfirst(a, b):
   951                             """Comparison function which always puts gluon first"""
   952                             if a == 21:
   953                                 return False
   954                             if b == 21:
   955                                 return True
   956                             return cmp(a, b)
   957                         d.ids = sorted(d.ids, cmp=cmp_gluonfirst)
   958                 elif len(absids) == 3:
   959                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   960                     def cmp_quarkslast(a, b):
   961                         """Comparison function which always puts quarks last"""
   962                         if is_quark(a):
   963                             return True
   964                         if is_quark(b):
   965                             return False
   966                         return cmp(a, b)
   967                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
   968             ## Squark/Slepton
   969             elif is_squark(pid) or is_slepton(pid):
   970                 def cmp_susy_quark_lepton(a, b):
   971                     if is_susy(a):
   972                         return False
   973                     if is_susy(b):
   974                         return True
   975                     if is_quark(a):
   976                         return False
   977                     if is_quark(b):
   978                         return True
   979                     return cmp(a, b)
   980                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   981                 ##                                                      Gluino       lepton
   982                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   983                 ## Squark
   984                 ##   2 body modes: Lepton Number Violated               quark     lepton
   985                 ##                 Baryon number violated               quark     quark
   986                 ## Slepton
   987                 ##   2 body modes: Lepton Number Violated               q or qbar
   988                 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   989             ## Higgs
   990             elif pid in (25, 26):
   991                 # TODO: Includes SUSY Higgses?
   992                 ## Higgs
   993                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   994                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   995                 ##   3 body modes:                                      colour neutral       q or qbar
   996                 if len(absids) == 3:
   997                     def cmp_quarkslast(a, b):
   998                         """Comparison function which always puts quarks last"""
   999                         if is_quark(a):
  1000                             return True
  1001                         if is_quark(b):
  1002                             return False
  1003                         return cmp(a, b)
  1004                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
  1005             elif is_gaugino(pid):
  1006                 # TODO: Is there actually anything to do here?
  1007                 ## Gaugino
  1008                 ##   2 body modes: Squark-quark                         q or sq
  1009                 ##                 Slepton-lepton                       l or sl
  1010                 ##
  1011                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
  1012                 ##                                                                           l or lbar
  1013                 if len(absids) == 3:
  1014                     def cmp_quarkslast(a, b):
  1015                         """Comparison function which always puts quarks last"""
  1016                         if is_quark(a):
  1017                             return True
  1018                         if is_quark(b):
  1019                             return False
  1020                         return cmp(a, b)
  1021                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
  1023             # TODO: Gaugino/Gluino
  1024             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
  1026             ## Pad out IDs list with zeros
  1027             ids = [0,0,0,0,0]
  1028             for i, pid in enumerate(d.ids):
  1029                 ids[i] = pid
  1030             ids = map(str, ids)
  1031             decayout += " ".join(ids) + "\n"
  1032         decayout = "%d\n" % decayout.count("\n") + decayout
  1033         out += decayout
  1035     ## Now the SUSY parameters
  1036     ## TANB, ALPHAH:
  1037     out += (fmte + " " + fmte + "\n") % (blocks["MINPAR"].value(3), blocks["ALPHA"].value())
  1038     ## Neutralino mixing matrix
  1039     nmix = blocks["NMIX"]
  1040     for i in xrange(1, 5):
  1041         out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (nmix[i,1], nmix[i,2], nmix[i,3], nmix[i,4])
  1042     ## Chargino mixing matrices V and U
  1043     vmix = blocks["VMIX"]
  1044     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (vmix[1,1], vmix[1,2], vmix[2,1], vmix[2,2])
  1045     umix = blocks["UMIX"]
  1046     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (umix[1,1], umix[1,2], umix[2,1], umix[2,2])
  1047     ## THETAT,THETAB,THETAL
  1048     import math
  1049     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (math.acos(blocks["STOPMIX"][1,1]),
  1050                                                             math.acos(blocks["SBOTMIX"][1,1]),
  1051                                                             math.acos(blocks["STAUMIX"][1,1]))
  1052     ## ATSS,ABSS,ALSS
  1053     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (blocks["AU"][3,3],
  1054                                                             blocks["AD"][3,3],
  1055                                                             blocks["AE"][3,3])
  1056     ## MUSS == sign(mu)
  1057     out += "%f\n" % blocks["MINPAR"][4]
  1059     ## RPV SUSY
  1060     isRPV = False
  1061     out += "%d\n" % isRPV
  1062     # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
  1063     # TODO: Get the index orderings right
  1064     # if isRPV: ...
  1066     return out
  1069 ###############################################################################
  1070 ## File-level functions
  1073 def readSLHAFile(spcfilename, **kwargs):
  1074     """
  1075     Read an SLHA file, returning dictionaries of blocks and decays.
  1077     Other keyword parameters are passed to readSLHA.
  1078     """
  1079     f = open(spcfilename, "r")
  1080     rtn = readSLHA(f.read(), kwargs)
  1081     f.close()
  1082     return rtn
  1085 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
  1086     """
  1087     Write an SLHA file from the supplied blocks and decays dicts.
  1089     Other keyword parameters are passed to writeSLHA.
  1090     """
  1091     f = open(spcfilename, "w")
  1092     f.write(writeSLHA(blocks, decays, kwargs))
  1093     f.close()
  1096 def readISAWIGFile(isafilename, **kwargs):
  1097     """
  1098     Read a spectrum definition from a file in the ISAWIG format, returning
  1099     dictionaries of blocks and decays. While this is not an SLHA format, it is
  1100     informally supported as a useful mechanism for converting ISAWIG spectra to
  1101     SLHA.
  1103     Other keyword parameters are passed to readSLHA.
  1104     """
  1105     f = open(isafilename, "r")
  1106     rtn = readISAWIG(f.read(), kwargs)
  1107     f.close()
  1108     return rtn
  1111 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
  1112     """
  1113     Write an ISAWIG file from the supplied blocks and decays dicts.
  1115     Other keyword parameters are passed to writeISAWIG.
  1116     """
  1117     f = open(isafilename, "w")
  1118     f.write(writeISAWIG(blocks, decays, kwargs))
  1119     f.close()
  1122 def writeISAJETFile(isafilename, blocks, decays, **kwargs):
  1123     """
  1124     Write an ISAJET file from the supplied blocks and decays dicts (see writeISAJET).
  1126     Other keyword parameters are passed to writeISAJET.
  1127     """
  1128     f = open(isafilename, "w")
  1129     f.write(writeISAWIG(blocks, decays, kwargs))
  1130     f.close()
  1134 ###############################################################################
  1135 ## Main function for module testing
  1138 if __name__ == "__main__":
  1139     import sys
  1140     for a in sys.argv[1:]:
  1141         if a.endswith(".isa"):
  1142             blocks, decays = readISAWIGFile(a)
  1143         else:
  1144             blocks, decays = readSLHAFile(a)
  1146         for bname, b in sorted(blocks.iteritems()):
  1147             print b
  1148             print
  1150         print blocks.keys()
  1152         print blocks["MASS"][25]
  1153         print
  1155         for p in sorted(decays.values()):
  1156             print p
  1157             print
  1159         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial