pyslha.py

Mon, 29 Apr 2013 14:54:16 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 29 Apr 2013 14:54:16 +0200
changeset 211
91f559c01cf7
parent 210
0f4f5472b7d8
child 212
3a6db3deedef
permissions
-rw-r--r--

Clean-up of precision formatting in output.

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

mercurial