pyslha.py

Mon, 29 Apr 2013 14:55:44 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 29 Apr 2013 14:55:44 +0200
changeset 212
3a6db3deedef
parent 211
91f559c01cf7
child 213
11aef4bd7802
permissions
-rw-r--r--

Extra thoughts about comment implementation

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

mercurial