pyslha.py

Sun, 28 Apr 2013 22:17:25 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 22:17:25 +0200
changeset 205
536967068fe2
parent 204
bef82eaef56e
child 207
2990422f1738
permissions
-rw-r--r--

Typo fix

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

mercurial