pyslha.py

Mon, 15 Jul 2013 16:50:38 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 15 Jul 2013 16:50:38 +0200
changeset 244
269bc4611cdd
parent 242
41ec8fb352b3
child 245
ca52fe2c0d71
permissions
-rw-r--r--

Add credit to ChangeLog

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

mercurial