pyslha.py

Mon, 15 Jul 2013 14:26:47 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 15 Jul 2013 14:26:47 +0200
changeset 240
84b1c74e0020
parent 239
42f74a1e96cf
child 241
401d5e64948f
permissions
-rw-r--r--

Fix kwargs passing from read/write*File to read/write*

     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 If you use PySLHA, for either model data handling or spectrum visualisation,
    32 please cite the paper: http://arxiv.org/abs/1305.4194
    34 TODOs:
    36   For 2.1.3:
    37    * Allow ignoring of missing MASS block in SLHA reader
    39   For 2.1.4:
    40    * In set_value, if first item is non-int, treat as None-indexed
    41    * Refine value string heuristic for strings with ints in them?
    43   For 2.2.0:
    44    * Add Sphinx docs.
    45    * Preserve comments from read -> write (needs full-line/inline comment
    46      separation). Can use separate comment dicts in Block and Decay, and
    47      attach a multiline .comment attr to the returned/written dicts.
    49   For 3.0.0:
    50    * Add a new Doc type which encapsulates the blocks, decays, xsections,
    51    etc. data sections and will act as a more robust single return value from
    52    the read commands / argument to the write commands. Can also handle
    53    document-level header comments this way and improve the __repr__ resulting
    54    from an interactive read() call without a target variable (i.e. hide the
    55    uninteresting OrderedDict boilerplate).
    57   Later/maybe:
    58    * Identify HERWIG decay matrix element to use in ISAWIG
    59    * Handle RPV SUSY in ISAWIG
    60 """
    62 __author__ = "Andy Buckley <andy.buckley@cern.ch>"
    63 __version__ = "2.1.2"
    67 ###############################################################################
    68 ## Private utility functions
    70 def _mkdict():
    71     """Try to return an empty ordered dict, but fall back to normal dict if necessary"""
    72     try:
    73         from collections import OrderedDict
    74         return OrderedDict()
    75     except:
    76         try:
    77             from ordereddict import OrderedDict
    78             return OrderedDict()
    79         except:
    80             return dict()
    82 def _autotype(var):
    83     """Automatically convert strings to numerical types if possible."""
    84     if type(var) is not str:
    85         return var
    86     if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
    87         return int(var)
    88     try:
    89         f = float(var)
    90         return f
    91     except ValueError:
    92         return var
    94 def _autostr(var, precision=8):
    95     """Automatically format numerical types as the right sort of string."""
    96     if type(var) is float:
    97         return ("% ." + str(precision) + "e") % var
    98     elif not hasattr(var, "__iter__"):
    99         return str(var)
   100     else:
   101         return ",".join(_autostr(subval) for subval in var)
   103 def _autotuple(a):
   104     """Automatically convert the supplied iterable to a scalar or tuple as appropriate."""
   105     if len(a) == 1:
   106         return a[0]
   107     return tuple(a)
   111 ###############################################################################
   112 ## Exceptions
   114 class AccessError(Exception):
   115     "Exception object to be raised when a SLHA block is accessed in an invalid way"
   116     def __init__(self, errmsg):
   117         self.msg = errmsg
   118     def __str__(self):
   119         return self.msg
   121 class ParseError(Exception):
   122     "Exception object to be raised when a spectrum file/string is malformed"
   123     def __init__(self, errmsg):
   124         self.msg = errmsg
   125     def __str__(self):
   126         return self.msg
   130 ###############################################################################
   131 ## The data block, decay and particle classes
   133 class Block(object):
   134     """
   135     Object representation of any BLOCK elements read from an SLHA file.
   137     Blocks have a name, may have an associated Q value, and contain a collection
   138     of data entries, each indexed by one or more keys. Entries in the dictionary
   139     are stored as numeric types (int or float) when a cast from the string in
   140     the file has been possible.
   142     Block is closely related to a Python dict (and, in fact, is implemented via
   143     an OrderedDict when possible). The preferred methods of entry access use the
   144     dict-like [] operator for getting and setting, and the keys() and items()
   145     methods for iteration. Purely iterating over the object behaves like keys(),
   146     as for an ordinary dict.
   148     Multiple (integer) indices are possible, especially for entries in mixing
   149     matrix blocks. These are now implemented in the natural way, e.g. for access
   150     to the (1,2) element of a mixing matrix block, use bmix[1,2] = 0.123 and
   151     print bmix[1,2]. The value() and set_value() functions behave
   152     similarly. Multi-element values are also permitted.
   154     It is possible, although not usual, to store unindexed values in a
   155     block. This is only supported when that entry is the only one in the block,
   156     and it is stored in the normal fashion but with None as the lookup key. The
   157     value() method may be used without a key argument to retrieve this value, if
   158     the is_single_valued() method returns True, and similarly the set_value()
   159     method may be used to set it if only one argument is supplied and the object
   160     is compatible.
   161     """
   162     def __init__(self, name, q=None, entries=None):
   163         self.name = name
   164         self.entries = _mkdict()
   165         if entries is not None:
   166             self.entries.update(entries)
   167         self.q = _autotype(q)
   169     # TODO: Rename? To what?
   170     def add_entry(self, args):
   171         """Add an entry to the block from an iterable (i.e. list or tuple) of
   172         strings, or from a whitespace-separated string.
   174         This method is just for convenience: it splits the single string
   175         argument if necessary and converts the list of strings into numeric
   176         types when possible. For the treatment of the resulting iterable see the
   177         set_value method.
   178         """
   179         ## If the argument is a single string, split it and proceed
   180         if type(args) is str:
   181             args = args.split()
   182         ## Check that the arg is an iterable
   183         if not hasattr(args, "__iter__"):
   184             raise AccessError("Block entries must be iterable")
   185         ## Auto-convert the types in the list
   186         args = map(_autotype, args)
   187         ## Re-join consecutive strings into single entries
   188         i = 0
   189         while i < len(args)-1:
   190             if type(args[i]) is str and type(args[i+1]) is str:
   191                 args[i] += " " + args[i+1]
   192                 del args[i+1]
   193                 continue
   194             i += 1
   195         ## Add the entry to the map, with appropriate indices
   196         self.set_value(*args)
   198     def is_single_valued(self):
   199         """Return true if there is only one entry, and it has no index: the
   200         'value()' attribute may be used in that case without an argument."""
   201         return len(self.entries) == 1 and self.entries.keys()[0] is None
   203     def value(self, key=None):
   204         """Get a value from the block with the supplied key.
   206         If no key is given, then the block must contain only one non-indexed
   207         value otherwise an AccessError exception will be raised.\
   208         """
   209         if key == None and not self.is_single_valued():
   210             raise AccessError("Tried to access unique value of multi-value block")
   211         return self.entries[key]
   213     def set_value(self, *args):
   214         """Set a value in the block via supplied key/value arguments.
   216         Indexing is determined automatically: any leading integers will be
   217         treated as a multi-dimensional index, with the remaining entries being a
   218         (potentially multi-dimensional) value. If all N args are ints, then the
   219         first N-1 are treated as the index and the Nth as the value.
   221         If there is only one arg it will be treated as the value of a
   222         single-valued block. In this case the block must already contain at most
   223         one non-indexed value otherwise an AccessError exception will be
   224         raised.\
   225         """
   226         if len(args) == 0:
   227             raise AccessError("set_value() called without arguments")
   228         elif len(args) == 1:
   229             if len(self.entries) > 0 and not self.is_single_valued():
   230                 raise AccessError("Tried to set a unique value on a multi-value block")
   231             self.entries[None] = args[0]
   232         else:
   233             ## Find the first non-integer -- all previous items are indices
   234             i_first_nonint = -1
   235             for i, x in enumerate(args):
   236                 if type(x) is not int:
   237                     i_first_nonint = i
   238                     break
   239             if i_first_nonint == 0:
   240                 raise AccessError("Attempted to set a block entry with a non-integer(s) index")
   241             else:
   242                 self.entries[_autotuple(args[:i_first_nonint])] = _autotuple(args[i_first_nonint:])
   244     def has_key(self, key):
   245         """Does the block have the given key?"""
   246         return self.entries.has_key(key)
   248     def keys(self):
   249         """Access the block item keys."""
   250         return self.entries.keys()
   252     def values(self):
   253         """Access the block item values."""
   254         return self.entries.values()
   256     def items(self, key=None):
   257         """Access the block items as (key,value) tuples.
   259         Note: The Python 3 dict attribute 'items()' is used rather than the
   260         'old' Python 2 'iteritems()' name for forward-looking compatibility.\
   261         """
   262         return self.entries.iteritems()
   264     def __len__(self):
   265         return len(self.entries)
   267     def __iter(self):
   268         return self.entries.__iter__()
   270     def __getitem__(self, key):
   271         return self.entries[key]
   273     def __setitem__(self, key, value):
   274         if key is not None and type(key) is not int and not all(type(x) is int for x in key):
   275             raise AccessError("Attempted to set a block entry with a non-integer(s) index")
   276         self.entries[key] = value
   278     def __cmp__(self, other):
   279         return cmp(self.name, other.name) and cmp(self.entries, other.entries)
   281     def __repr__(self):
   282         s = self.name
   283         if self.q is not None:
   284             s += " (Q=%s)" % self.q
   285         s += " { " + "; ".join(_autostr(k) + " : " + _autostr(v) for k, v in self.items()) + " }"
   286         return s
   289 class Decay(object):
   290     """
   291     Object representing a decay entry on a particle decribed by the SLHA file.
   292     'Decay' objects are not a direct representation of a DECAY block in an SLHA
   293     file... that role, somewhat confusingly, is taken by the Particle class.
   295     Decay objects have three properties: a branching ratio, br, an nda number
   296     (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
   297     decay occurs. The PDG ID of the particle whose decay this represents may
   298     also be stored, but this is normally known via the Particle in which the
   299     decay is stored.
   300     """
   301     def __init__(self, br, nda, ids, parentid=None):
   302         self.parentid = parentid
   303         self.br = br
   304         self.nda = nda
   305         self.ids = ids
   306         assert(self.nda == len(self.ids))
   308     def __cmp__(self, other):
   309         return cmp(other.br, self.br)
   311     def __repr__(self):
   312         return "% .8e %s" % (self.br, self.ids)
   315 class Particle(object):
   316     """
   317     Representation of a single, specific particle, decay block from an SLHA
   318     file.  These objects are not themselves called 'Decay', since that concept
   319     applies more naturally to the various decays found inside this
   320     object. Particle classes store the PDG ID (pid) of the particle being
   321     represented, and optionally the mass (mass) and total decay width
   322     (totalwidth) of that particle in the SLHA scenario. Masses may also be found
   323     via the MASS block, from which the Particle.mass property is filled, if at
   324     all. They also store a list of Decay objects (decays) which are probably the
   325     item of most interest.
   326     """
   327     def __init__(self, pid, totalwidth=None, mass=None):
   328         self.pid = pid
   329         self.totalwidth = totalwidth
   330         self.mass = mass
   331         self.decays = []
   333     def add_decay(self, br, nda, ids):
   334         self.decays.append(Decay(br, nda, ids))
   335         self.decays.sort()
   337     def __cmp__(self, other):
   338         if abs(self.pid) == abs(other.pid):
   339             return cmp(self.pid, other.pid)
   340         return cmp(abs(self.pid), abs(other.pid))
   342     def __repr__(self):
   343         s = str(self.pid)
   344         if self.mass is not None:
   345             s += " : mass = %.8e GeV" % self.mass
   346         if self.totalwidth is not None:
   347             s += " : total width = %.8e GeV" % self.totalwidth
   348         for d in self.decays:
   349             if d.br > 0.0:
   350                 s += "\n  %s" % d
   351         return s
   355 ###############################################################################
   356 ## SLHA parsing and writing functions
   358 def readSLHA(spcstr, ignorenobr=False):
   359     """
   360     Read an SLHA definition from a string, returning dictionaries of blocks and
   361     decays.
   363     If the ignorenobr parameter is True, do not store decay entries with a
   364     branching ratio of zero.
   365     """
   366     blocks = _mkdict()
   367     decays = _mkdict()
   368     #
   369     import re
   370     currentblock = None
   371     currentdecay = None
   372     for line in spcstr.splitlines():
   373         ## Handle (ignore) comment lines
   374         # TODO: Store block/entry comments
   375         if line.startswith("#"):
   376             continue
   377         if "#" in line:
   378             line = line[:line.index("#")]
   379         ## Ignore empty lines (after comment removal and whitespace trimming)
   380         if not line.strip():
   381             continue
   383         ## Section header lines start with a non-whitespace character, data lines have a whitespace indent
   384         # TODO: Are tabs also allowed for indents? Check the SLHA standard.
   385         if not line.startswith(" "):
   386             # TODO: Should we now strip the line to remove any trailing whitespace?
   388             ## Handle BLOCK start lines
   389             if line.upper().startswith("BLOCK"):
   390                 #print line
   391                 match = re.match(r"BLOCK\s+(\w+)(\s+Q\s*=\s*.+)?", line.upper())
   392                 if not match:
   393                     continue
   394                 blockname = match.group(1)
   395                 qstr = match.group(2)
   396                 if qstr is not None:
   397                     qstr = qstr[qstr.find("=")+1:].strip()
   398                 currentblock = blockname
   399                 currentdecay = None
   400                 blocks[blockname] = Block(blockname, q=qstr)
   401             ## Handle DECAY start lines
   402             elif line.upper().startswith("DECAY"):
   403                 match = re.match(r"DECAY\s+(-?\d+)\s+([\d\.E+-]+|NAN).*", line.upper())
   404                 if not match:
   405                     continue
   406                 pdgid = int(match.group(1))
   407                 width = float(match.group(2)) if match.group(2) != "NAN" else None
   408                 currentblock = "DECAY"
   409                 currentdecay = pdgid
   410                 decays[pdgid] = Particle(pdgid, width)
   411             ## Handle unknown section type start lines (and continue ignoring until a non-header line is found)
   412             elif type(_autotype(line.split()[0])) is str:
   413                 import sys
   414                 sys.stderr.write("Ignoring unknown section type: %s\n" % line.split()[0])
   415                 currentblock = None
   416                 currentdecay = None
   418         ## This non-empty line starts with an indent, hence must be an in-section data line
   419         else:
   420             # TODO: Should we now strip the line to remove the indent (and any trailing whitespace)?
   421             if currentblock is not None:
   422                 items = line.split()
   423                 if len(items) < 1:
   424                     continue
   425                 if currentblock != "DECAY":
   426                     blocks[currentblock].add_entry(items)
   427                 # TODO: Add handling of XSECTION if/when standardised
   428                 else:
   429                     br = float(items[0]) if items[0].upper() != "NAN" else None
   430                     nda = int(items[1])
   431                     ids = map(int, items[2:])
   432                     if br > 0.0 or not ignorenobr: # br == None is < 0
   433                         decays[currentdecay].add_decay(br, nda, ids)
   435     ## Try to populate Particle masses from the MASS block
   436     # print blocks.keys()
   437     try:
   438         for pid in blocks["MASS"].keys():
   439             if decays.has_key(pid):
   440                 decays[pid].mass = blocks["MASS"][pid]
   441     except:
   442         raise ParseError("No MASS block found: cannot populate particle masses")
   444     return blocks, decays
   448 def writeSLHABlocks(blocks, precision=8):
   449     """Return an SLHA definition as a string, from the supplied blocks dict."""
   450     sep = 3 * " "
   451     blockstrs = []
   452     for bname, b in blocks.iteritems():
   453         namestr = b.name
   454         if b.q is not None:
   455             namestr += " Q= " + _autostr(float(b.q), precision)
   456         blockstr = "BLOCK %s\n" % namestr
   457         entrystrs = []
   458         for k, v in b.items():
   459             entrystr = ""
   460             if type(k) == tuple:
   461                 entrystr += sep.join(_autostr(i, precision) for i in k)
   462             elif k is not None:
   463                 entrystr += _autostr(k, precision)
   464             entrystr += sep + _autostr(v, precision)
   465             entrystrs.append(entrystr)
   466         blockstr += "\n".join(entrystrs)
   467         blockstrs.append(blockstr)
   468     return "\n\n".join(blockstrs)
   472 def writeSLHADecays(decays, ignorenobr=False, precision=8):
   473     """Return an SLHA decay definition as a string, from the supplied decays dict."""
   474     sep = 3 * " "
   475     blockstrs = []
   476     for pid, particle in decays.iteritems():
   477         blockstr = ("DECAY %d " % particle.pid) + _autostr(particle.totalwidth or -1, precision) + "\n"
   478         decaystrs = []
   479         for d in particle.decays:
   480             if d.br > 0.0 or not ignorenobr:
   481                 products_str = sep.join("% d" % i for i in d.ids)
   482                 decaystr = sep + _autostr(d.br, precision) + sep + str(len(d.ids)) + sep + products_str
   483                 decaystrs.append(decaystr)
   484         blockstr += "\n".join(decaystrs)
   485         blockstrs.append(blockstr)
   486     return "\n\n".join(blockstrs)
   490 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
   491     """Return an SLHA definition as a string, from the supplied blocks and decays dicts."""
   492     ss = [x for x in (writeSLHABlocks(blocks, precision), writeSLHADecays(decays, ignorenobr, precision)) if x]
   493     return "\n\n".join(ss)
   497 ###############################################################################
   498 ## PDG <-> HERWIG particle ID code translations for ISAWIG handling
   500 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
   501 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   502 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   503 _HERWIGID2PDGID = {}
   504 _HERWIGID2PDGID[7]   = -1
   505 _HERWIGID2PDGID[8]   = -2
   506 _HERWIGID2PDGID[9]   = -3
   507 _HERWIGID2PDGID[10]  = -4
   508 _HERWIGID2PDGID[11]  = -5
   509 _HERWIGID2PDGID[12]  = -6
   510 _HERWIGID2PDGID[13]  =  21
   511 _HERWIGID2PDGID[59]  =  22
   512 _HERWIGID2PDGID[121] =  11
   513 _HERWIGID2PDGID[122] =  12
   514 _HERWIGID2PDGID[123] =  13
   515 _HERWIGID2PDGID[124] =  14
   516 _HERWIGID2PDGID[125] =  15
   517 _HERWIGID2PDGID[126] =  16
   518 _HERWIGID2PDGID[127] = -11
   519 _HERWIGID2PDGID[128] = -12
   520 _HERWIGID2PDGID[129] = -13
   521 _HERWIGID2PDGID[130] = -14
   522 _HERWIGID2PDGID[131] = -15
   523 _HERWIGID2PDGID[132] = -16
   524 _HERWIGID2PDGID[198] =  24 # W+
   525 _HERWIGID2PDGID[199] = -24 # W-
   526 _HERWIGID2PDGID[200] =  23 # Z0
   527 _HERWIGID2PDGID[201] =  25 ## SM HIGGS
   528 _HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   529 _HERWIGID2PDGID[204] =  35 ## HIGGSH0
   530 _HERWIGID2PDGID[205] =  36 ## HIGGSA0
   531 _HERWIGID2PDGID[206] =  37 ## HIGGS+
   532 _HERWIGID2PDGID[207] = -37 ## HIGGS-
   533 _HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   534 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   535 _HERWIGID2PDGID[402] =  1000002 ## SSULBR
   536 _HERWIGID2PDGID[408] = -1000002 ## SSUL
   537 _HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   538 _HERWIGID2PDGID[409] = -1000003 ## SSSL
   539 _HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   540 _HERWIGID2PDGID[410] = -1000004 ## SSCL
   541 _HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   542 _HERWIGID2PDGID[411] = -1000005 ## SSB1
   543 _HERWIGID2PDGID[406] =  1000006 ## SST1BR
   544 _HERWIGID2PDGID[412] = -1000006 ## SST1
   545 _HERWIGID2PDGID[413] =  2000001 ## SSDR
   546 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   547 _HERWIGID2PDGID[414] =  2000002 ## SSUR
   548 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
   549 _HERWIGID2PDGID[415] =  2000003 ## SSSR
   550 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   551 _HERWIGID2PDGID[416] =  2000004 ## SSCR
   552 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   553 _HERWIGID2PDGID[417] =  2000005 ## SSB2
   554 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   555 _HERWIGID2PDGID[418] =  2000006 ## SST2
   556 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
   557 _HERWIGID2PDGID[425] =  1000011 ## SSEL-
   558 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
   559 _HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   560 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   561 _HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   562 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   563 _HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   564 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   565 _HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   566 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   567 _HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   568 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   569 _HERWIGID2PDGID[437] =  2000011 ## SSEL-
   570 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
   571 _HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   572 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   573 _HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   574 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   575 _HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   576 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   577 _HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   578 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   579 _HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   580 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   581 _HERWIGID2PDGID[449] =  1000021 ## GLUINO
   582 _HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   583 _HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   584 _HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   585 _HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   586 _HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   587 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   588 _HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   589 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   590 _HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   592 def herwigid2pdgid(hwid):
   593     """
   594     Convert a particle ID code in the HERWIG internal IDHW format (as used by
   595     ISAWIG) into its equivalent in the standard PDG ID code definition.
   596     """
   597     return _HERWIGID2PDGID.get(hwid, hwid)
   600 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   601 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   602 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   603 _PDGID2HERWIGID = {}
   604 _PDGID2HERWIGID[      -1] = 7
   605 _PDGID2HERWIGID[      -2] = 8
   606 _PDGID2HERWIGID[      -3] = 9
   607 _PDGID2HERWIGID[      -4] = 10
   608 _PDGID2HERWIGID[      -5] = 11
   609 _PDGID2HERWIGID[      -6] = 12
   610 _PDGID2HERWIGID[      21] = 13
   611 _PDGID2HERWIGID[      22] = 59
   612 _PDGID2HERWIGID[      11] = 121
   613 _PDGID2HERWIGID[      12] = 122
   614 _PDGID2HERWIGID[      13] = 123
   615 _PDGID2HERWIGID[      14] = 124
   616 _PDGID2HERWIGID[      15] = 125
   617 _PDGID2HERWIGID[      16] = 126
   618 _PDGID2HERWIGID[     -11] = 127
   619 _PDGID2HERWIGID[     -12] = 128
   620 _PDGID2HERWIGID[     -13] = 129
   621 _PDGID2HERWIGID[     -14] = 130
   622 _PDGID2HERWIGID[     -15] = 131
   623 _PDGID2HERWIGID[     -16] = 132
   624 _PDGID2HERWIGID[      24] = 198 ## W+
   625 _PDGID2HERWIGID[     -24] = 199 ## W-
   626 _PDGID2HERWIGID[      23] = 200 ## Z
   627 _PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
   628 _PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   629 _PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   630 _PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   631 _PDGID2HERWIGID[      37] = 206 ## HIGGS+
   632 _PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   633 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   634 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   635 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   636 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
   637 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   638 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
   639 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   640 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
   641 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   642 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
   643 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   644 _PDGID2HERWIGID[-1000006] = 412 ## SST1
   645 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   646 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   647 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   648 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   649 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   650 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   651 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   652 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   653 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   654 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   655 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
   656 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   657 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   658 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   659 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   660 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   661 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   662 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   663 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   664 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   665 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   666 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   667 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   668 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   669 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   670 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   671 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   672 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   673 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   674 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   675 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   676 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   677 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   678 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   679 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   680 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   681 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   682 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   683 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   684 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   685 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   686 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   687 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   688 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   689 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   690 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   692 def pdgid2herwigid(pdgid):
   693     """
   694     Convert a particle ID code in the standard PDG ID code definition into
   695     its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
   696     """
   697     return _PDGID2HERWIGID.get(pdgid, pdgid)
   700 ###############################################################################
   701 ## ISAWIG format reading/writing
   703 def readISAWIG(isastr, ignorenobr=False):
   704     """
   705     Read a spectrum definition from a string in the ISAWIG format, returning
   706     dictionaries of blocks and decays. While this is not an SLHA format, it is
   707     informally supported as a useful mechanism for converting ISAWIG spectra to
   708     SLHA.
   710     ISAWIG parsing based on the HERWIG SUSY specification format, from
   711     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   713     If the ignorenobr parameter is True, do not store decay entries with a
   714     branching ratio of zero.
   715     """
   717     blocks = _mkdict()
   718     decays = _mkdict()
   719     LINES = isastr.splitlines()
   721     def getnextvalidline():
   722         while LINES:
   723             s = LINES.pop(0).strip()
   724             # print "*", s, "*"
   725             ## Return None if EOF reached
   726             if len(s) == 0:
   727                 continue
   728             ## Strip comments
   729             if "#" in s:
   730                 s = s[:s.index("#")].strip()
   731             ## Return if non-empty
   732             if len(s) > 0:
   733                 return s
   735     def getnextvalidlineitems():
   736         return map(_autotype, getnextvalidline().split())
   738     ## Populate MASS block and create decaying particle objects
   739     masses = Block("MASS")
   740     numentries = int(getnextvalidline())
   741     for i in xrange(numentries):
   742         hwid, mass, lifetime = getnextvalidlineitems()
   743         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   744         pdgid = herwigid2pdgid(hwid)
   745         masses[pdgid] = mass
   746         decays[pdgid] = Particle(pdgid, width, mass)
   747         #print pdgid, mass, width
   748     blocks["MASS"] = masses
   750     ## Populate decays
   751     for n in xrange(numentries):
   752         numdecays = int(getnextvalidline())
   753         for d in xrange(numdecays):
   754             #print n, numentries-1, d, numdecays-1
   755             decayitems = getnextvalidlineitems()
   756             hwid = decayitems[0]
   757             pdgid = herwigid2pdgid(hwid)
   758             br = decayitems[1]
   759             nme = decayitems[2]
   760             daughter_hwids = decayitems[3:]
   761             daughter_pdgids = []
   762             for hw in daughter_hwids:
   763                 if hw != 0:
   764                     daughter_pdgids.append(herwigid2pdgid(hw))
   765             if not decays.has_key(pdgid):
   766                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   767                 decays[pdgid] = Particle(pdgid)
   768             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   771     ## Now the SUSY parameters
   772     TANB, ALPHAH = getnextvalidlineitems()
   773     blocks["MINPAR"] = Block("MINPAR")
   774     blocks["MINPAR"][3] = TANB
   775     blocks["ALPHA"] = Block("ALPHA")
   776     blocks["ALPHA"].set_value(ALPHAH)
   777     #
   778     ## Neutralino mixing matrix
   779     blocks["NMIX"] = Block("NMIX")
   780     for i in xrange(1, 5):
   781         nmix_i = getnextvalidlineitems()
   782         for j, v in enumerate(nmix_i):
   783             blocks["NMIX"][i, j+1] = v
   784     #
   785     ## Chargino mixing matrices V and U
   786     blocks["VMIX"] = Block("VMIX")
   787     vmix = getnextvalidlineitems()
   788     blocks["VMIX"][1, 1] = vmix[0]
   789     blocks["VMIX"][1, 2] = vmix[1]
   790     blocks["VMIX"][2, 1] = vmix[2]
   791     blocks["VMIX"][2, 2] = vmix[3]
   792     blocks["UMIX"] = Block("UMIX")
   793     umix = getnextvalidlineitems()
   794     blocks["UMIX"][1, 1] = umix[0]
   795     blocks["UMIX"][1, 2] = umix[1]
   796     blocks["UMIX"][2, 1] = umix[2]
   797     blocks["UMIX"][2, 2] = umix[3]
   798     #
   799     THETAT, THETAB, THETAL = getnextvalidlineitems()
   800     import math
   801     blocks["STOPMIX"] = Block("STOPMIX")
   802     blocks["STOPMIX"][1, 1] =  math.cos(THETAT)
   803     blocks["STOPMIX"][1, 2] = -math.sin(THETAT)
   804     blocks["STOPMIX"][2, 1] =  math.sin(THETAT)
   805     blocks["STOPMIX"][2, 2] =  math.cos(THETAT)
   806     blocks["SBOTMIX"] = Block("SBOTMIX")
   807     blocks["SBOTMIX"][1, 1] =  math.cos(THETAB)
   808     blocks["SBOTMIX"][1, 2] = -math.sin(THETAB)
   809     blocks["SBOTMIX"][2, 1] =  math.sin(THETAB)
   810     blocks["SBOTMIX"][2, 2] =  math.cos(THETAB)
   811     blocks["STAUMIX"] = Block("STAUMIX")
   812     blocks["STAUMIX"][1, 1] =  math.cos(THETAL)
   813     blocks["STAUMIX"][1, 2] = -math.sin(THETAL)
   814     blocks["STAUMIX"][2, 1] =  math.sin(THETAL)
   815     blocks["STAUMIX"][2, 2] =  math.cos(THETAL)
   816     #
   817     ATSS, ABSS, ALSS = getnextvalidlineitems()
   818     blocks["AU"] = Block("AU")
   819     blocks["AU"][3, 3] = ATSS
   820     blocks["AD"] = Block("AD")
   821     blocks["AD"][3, 3] = ABSS
   822     blocks["AE"] = Block("AE")
   823     blocks["AE"][3, 3] = ALSS
   824     #
   825     MUSS = getnextvalidlineitems()[0]
   826     blocks["MINPAR"][4] = MUSS
   827     #
   829     # TODO: Parse RPV boolean and couplings into SLHA2 blocks
   831     return blocks, decays
   834 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
   835     """
   836     Return a SUSY spectrum definition in the format produced by ISAWIG for inut to HERWIG
   837     as a string, from the supplied SLHA blocks and decays dicts.
   839     ISAWIG parsing based on the HERWIG SUSY specification format, from
   840     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   842     If the ignorenobr parameter is True, do not write decay entries with a
   843     branching ratio of zero.
   844     """
   845     masses = blocks["MASS"]
   847     ## Init output string
   848     out = ""
   850     ## First write out masses section:
   851     ##   Number of SUSY + top particles
   852     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   853     ##   repeated for each particle
   854     ## IDHW is the HERWIG identity code.
   855     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   856     massout = ""
   857     for pid in masses.keys():
   858         lifetime = -1
   859         try:
   860             width = decays[pid].totalwidth
   861             if width and width > 0:
   862                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   863         except:
   864             pass
   865         massout += ("%d " % pdgid2herwigid(pid)) + _autostr(masses[pid], precision) + " " + _autostr(lifetime, precision) + "\n"
   866     out += "%d\n" % massout.count("\n")
   867     out += massout
   869     assert(len(masses) == len(decays))
   871     ## Next each particles decay modes together with their branching ratios and matrix element codes
   872     ##   Number of decay modes for a given particle (IDK)
   873     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   874     ##     repeated for each mode.
   875     ##   Repeated for each particle.
   876     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   877     ## the decay mode. NME is a code for the matrix element to be used, either from the
   878     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   879     for i, pid in enumerate(decays.keys()):
   880         # if not decays.has_key(pid):
   881         #     continue
   882         hwid = pdgid2herwigid(pid)
   883         decayout = ""
   884         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   885         for i_d, d in enumerate(decays[pid].decays):
   886             ## Skip decay if it has no branching ratio
   887             if ignorenobr and d.br == 0:
   888                 continue
   890             ## Identify decay matrix element to use
   891             ## From std HW docs, or from this pair:
   892             ## Two new matrix element codes have been added for these new decays:
   893             ##    NME =	200 	3 body top quark via charged Higgs
   894             ##    	300 	3 body R-parity violating gaugino and gluino decays
   895             nme = 0
   896             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   897             if abs(pid) in (6, 12):
   898                 nme = 100
   899             ## Extra SUSY MEs
   900             if len(d.ids) == 3:
   901                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   902                 pass
   903             decayout += ("%d " % hwid) + _autostr(d.br, precision) + (" %d " % nme)
   905             def is_quark(pid):
   906                 return (abs(pid) in range(1, 7))
   908             def is_lepton(pid):
   909                 return (abs(pid) in range(11, 17))
   911             def is_squark(pid):
   912                 if abs(pid) in range(1000001, 1000007):
   913                     return True
   914                 if abs(pid) in range(2000001, 2000007):
   915                     return True
   916                 return False
   918             def is_slepton(pid):
   919                 if abs(pid) in range(1000011, 1000017):
   920                     return True
   921                 if abs(pid) in range(2000011, 2000016, 2):
   922                     return True
   923                 return False
   925             def is_gaugino(pid):
   926                 if abs(pid) in range(1000022, 1000026):
   927                     return True
   928                 if abs(pid) in (1000035, 1000037):
   929                     return True
   930                 return False
   932             def is_susy(pid):
   933                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   935             absids = map(abs, d.ids)
   937             ## Order decay products as required by HERWIG
   938             ## Top
   939             if abs(pid) == 6:
   940                 def cmp_bottomlast(a, b):
   941                     """Comparison function which always puts b/bbar last"""
   942                     if abs(a) == 5:
   943                         return True
   944                     if abs(b) == 5:
   945                         return False
   946                     return cmp(a, b)
   947                 if len(absids) == 2:
   948                     ## 2 body mode, to Higgs: Higgs; Bottom
   949                     if (25 in absids or 26 in absids) and 5 in absids:
   950                         d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   951                 elif len(absids) == 3:
   952                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   953                     if 37 in absids or 23 in absids:
   954                         d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   955             ## Gluino
   956             elif abs(pid) == 1000021:
   957                 if len(absids) == 2:
   958                     ## 2 body mode
   959                     ## without gluon: any order
   960                     ## with gluon: gluon; colour neutral
   961                     if 21 in absids:
   962                         def cmp_gluonfirst(a, b):
   963                             """Comparison function which always puts gluon first"""
   964                             if a == 21:
   965                                 return False
   966                             if b == 21:
   967                                 return True
   968                             return cmp(a, b)
   969                         d.ids = sorted(d.ids, cmp=cmp_gluonfirst)
   970                 elif len(absids) == 3:
   971                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   972                     def cmp_quarkslast(a, b):
   973                         """Comparison function which always puts quarks last"""
   974                         if is_quark(a):
   975                             return True
   976                         if is_quark(b):
   977                             return False
   978                         return cmp(a, b)
   979                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
   980             ## Squark/Slepton
   981             elif is_squark(pid) or is_slepton(pid):
   982                 def cmp_susy_quark_lepton(a, b):
   983                     if is_susy(a):
   984                         return False
   985                     if is_susy(b):
   986                         return True
   987                     if is_quark(a):
   988                         return False
   989                     if is_quark(b):
   990                         return True
   991                     return cmp(a, b)
   992                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   993                 ##                                                      Gluino       lepton
   994                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   995                 ## Squark
   996                 ##   2 body modes: Lepton Number Violated               quark     lepton
   997                 ##                 Baryon number violated               quark     quark
   998                 ## Slepton
   999                 ##   2 body modes: Lepton Number Violated               q or qbar
  1000                 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
  1001             ## Higgs
  1002             elif pid in (25, 26):
  1003                 # TODO: Includes SUSY Higgses?
  1004                 ## Higgs
  1005                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
  1006                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
  1007                 ##   3 body modes:                                      colour neutral       q or qbar
  1008                 if len(absids) == 3:
  1009                     def cmp_quarkslast(a, b):
  1010                         """Comparison function which always puts quarks last"""
  1011                         if is_quark(a):
  1012                             return True
  1013                         if is_quark(b):
  1014                             return False
  1015                         return cmp(a, b)
  1016                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
  1017             elif is_gaugino(pid):
  1018                 # TODO: Is there actually anything to do here?
  1019                 ## Gaugino
  1020                 ##   2 body modes: Squark-quark                         q or sq
  1021                 ##                 Slepton-lepton                       l or sl
  1022                 ##
  1023                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
  1024                 ##                                                                           l or lbar
  1025                 if len(absids) == 3:
  1026                     def cmp_quarkslast(a, b):
  1027                         """Comparison function which always puts quarks last"""
  1028                         if is_quark(a):
  1029                             return True
  1030                         if is_quark(b):
  1031                             return False
  1032                         return cmp(a, b)
  1033                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
  1035             # TODO: Gaugino/Gluino
  1036             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
  1038             ## Pad out IDs list with zeros
  1039             ids = [0,0,0,0,0]
  1040             for i, pid in enumerate(d.ids):
  1041                 ids[i] = pid
  1042             ids = map(str, ids)
  1043             decayout += " ".join(ids) + "\n"
  1044         decayout = "%d\n" % decayout.count("\n") + decayout
  1045         out += decayout
  1047     ## Now the SUSY parameters
  1048     ## TANB, ALPHAH:
  1049     out += _autostr(blocks["MINPAR"][3], precision) + " " + _autostr(blocks["ALPHA"].value(), precision) + "\n"
  1050     ## Neutralino mixing matrix
  1051     nmix = blocks["NMIX"]
  1052     for i in xrange(1, 5):
  1053         out += _autostr(nmix[i,1], precision) + " " + \
  1054                _autostr(nmix[i,2], precision) + " " + \
  1055                _autostr(nmix[i,3], precision) + " " + \
  1056                _autostr(nmix[i,4], precision) + "\n"
  1057     ## Chargino mixing matrices V and U
  1058     vmix = blocks["VMIX"]
  1059     out += _autostr(vmix[1,1], precision) + " " + \
  1060            _autostr(vmix[1,2], precision) + " " + \
  1061            _autostr(vmix[2,1], precision) + " " + \
  1062            _autostr(vmix[2,2], precision) + "\n"
  1063     umix = blocks["UMIX"]
  1064     out += _autostr(umix[1,1], precision) + " " + \
  1065            _autostr(umix[1,2], precision) + " " + \
  1066            _autostr(umix[2,1], precision) + " " + \
  1067            _autostr(umix[2,2], precision) + "\n"
  1068     ## THETAT,THETAB,THETAL
  1069     import math
  1070     out += _autostr(math.acos(blocks["STOPMIX"][1,1]), precision) + " " + \
  1071            _autostr(math.acos(blocks["SBOTMIX"][1,1]), precision) + " " + \
  1072            _autostr(math.acos(blocks["STAUMIX"][1,1]), precision) + "\n"
  1073     ## ATSS,ABSS,ALSS
  1074     out += _autostr(blocks["AU"][3,3], precision) + " " + \
  1075            _autostr(blocks["AD"][3,3], precision) + " " + \
  1076            _autostr(blocks["AE"][3,3], precision) + "\n"
  1077     ## MUSS == sign(mu)
  1078     out += "%f\n" % blocks["MINPAR"][4]
  1080     ## RPV SUSY
  1081     isRPV = False
  1082     out += "%d\n" % isRPV
  1083     # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
  1084     # TODO: Get the index orderings right
  1085     # if isRPV: ...
  1087     return out
  1091 ###############################################################################
  1092 ## File-level read/write functions
  1094 def readSLHAFile(spcfilename, **kwargs):
  1095     """
  1096     Read an SLHA file, returning dictionaries of blocks and decays.
  1098     Other keyword parameters are passed to readSLHA.
  1099     """
  1100     f = open(spcfilename, "r")
  1101     rtn = readSLHA(f.read(), **kwargs)
  1102     f.close()
  1103     return rtn
  1106 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
  1107     """
  1108     Write an SLHA file from the supplied blocks and decays dicts.
  1110     Other keyword parameters are passed to writeSLHA.
  1111     """
  1112     f = open(spcfilename, "w")
  1113     f.write(writeSLHA(blocks, decays, **kwargs))
  1114     f.close()
  1117 def readISAWIGFile(isafilename, **kwargs):
  1118     """
  1119     Read a spectrum definition from a file in the ISAWIG format, returning
  1120     dictionaries of blocks and decays. While this is not an SLHA format, it is
  1121     informally supported as a useful mechanism for converting ISAWIG spectra to
  1122     SLHA.
  1124     Other keyword parameters are passed to readSLHA.
  1125     """
  1126     f = open(isafilename, "r")
  1127     rtn = readISAWIG(f.read(), **kwargs)
  1128     f.close()
  1129     return rtn
  1132 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
  1133     """
  1134     Write an ISAWIG file from the supplied blocks and decays dicts.
  1136     Other keyword parameters are passed to writeISAWIG.
  1137     """
  1138     f = open(isafilename, "w")
  1139     f.write(writeISAWIG(blocks, decays, **kwargs))
  1140     f.close()
  1144 ###############################################################################
  1145 ## Main function for module testing
  1147 if __name__ == "__main__":
  1148     import sys
  1149     for a in sys.argv[1:]:
  1150         if a.endswith(".isa"):
  1151             blocks, decays = readISAWIGFile(a)
  1152         else:
  1153             blocks, decays = readSLHAFile(a)
  1155         for bname, b in sorted(blocks.iteritems()):
  1156             print b
  1157             print
  1159         print blocks.keys()
  1161         print blocks["MASS"][25]
  1162         print
  1164         for p in sorted(decays.values()):
  1165             print p
  1166             print
  1168         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial