pyslha.py

Sun, 28 Apr 2013 22:15:24 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 22:15:24 +0200
changeset 204
bef82eaef56e
parent 203
ce90a0dace07
child 205
536967068fe2
permissions
-rw-r--r--

Adding a dict arg (compatible with dict.update) to the Block constructor.

     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 keys(self):
   178         """Access the block item keys."""
   179         return self.entries.keys()
   181     def items(self, key=None):
   182         """Access the block items as (key,value) tuples.
   184         Note: The Python 3 dict attribute 'items()' is used rather than the
   185         'old' Python 2 'iteritems()' name for forward-looking compatibility.\
   186         """
   187         return self.entries.iteritems()
   189     def __len__(self):
   190         return len(self.entries)
   192     def __iter(self):
   193         return self.entries.__iter__()
   195     def __getitem__(self, key):
   196         return self.entries[key]
   198     def __setitem__(self, key, value):
   199         self.entries[key] = value
   201     def __cmp__(self, other):
   202         return cmp(self.name, other.name) and cmp(self.entries, other.entries)
   204     def __repr__(self):
   205         s = self.name
   206         if self.q is not None:
   207             s += " (Q=%s)" % self.q
   208         def _format_kv(k, v):
   209             if type(k) is not tuple:
   210                 s = "%r" % k
   211             else:
   212                 s = ",".join("%r" % subindex for subindex in k)
   213             s += " : %r" % v
   214             return s
   215         s += " { " + "; ".join(_format_kv(k, v) for k, v in self.items()) + " }"
   216         return s
   219 class Decay(object):
   220     """
   221     Object representing a decay entry on a particle decribed by the SLHA file.
   222     'Decay' objects are not a direct representation of a DECAY block in an SLHA
   223     file... that role, somewhat confusingly, is taken by the Particle class.
   225     Decay objects have three properties: a branching ratio, br, an nda number
   226     (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
   227     decay occurs. The PDG ID of the particle whose decay this represents may
   228     also be stored, but this is normally known via the Particle in which the
   229     decay is stored.
   230     """
   231     def __init__(self, br, nda, ids, parentid=None):
   232         self.parentid = parentid
   233         self.br = br
   234         self.nda = nda
   235         self.ids = ids
   236         assert(self.nda == len(self.ids))
   238     def __cmp__(self, other):
   239         return cmp(other.br, self.br)
   241     def __str__(self):
   242         return "%.8e %s" % (self.br, self.ids)
   244     def __repr__(self):
   245         return self.__str__()
   248 class Particle(object):
   249     """
   250     Representation of a single, specific particle, decay block from an SLHA
   251     file.  These objects are not themselves called 'Decay', since that concept
   252     applies more naturally to the various decays found inside this
   253     object. Particle classes store the PDG ID (pid) of the particle being
   254     represented, and optionally the mass (mass) and total decay width
   255     (totalwidth) of that particle in the SLHA scenario. Masses may also be found
   256     via the MASS block, from which the Particle.mass property is filled, if at
   257     all. They also store a list of Decay objects (decays) which are probably the
   258     item of most interest.
   259     """
   260     def __init__(self, pid, totalwidth=None, mass=None):
   261         self.pid = pid
   262         self.totalwidth = totalwidth
   263         self.mass = mass
   264         self.decays = []
   266     def add_decay(self, br, nda, ids):
   267         self.decays.append(Decay(br, nda, ids))
   268         self.decays.sort()
   270     def __cmp__(self, other):
   271         if abs(self.pid) == abs(other.pid):
   272             return cmp(self.pid, other.pid)
   273         return cmp(abs(self.pid), abs(other.pid))
   275     def __str__(self):
   276         s = str(self.pid)
   277         if self.mass is not None:
   278             s += " : mass = %.8e GeV" % self.mass
   279         if self.totalwidth is not None:
   280             s += " : total width = %.8e GeV" % self.totalwidth
   281         for d in self.decays:
   282             if d.br > 0.0:
   283                 s += "\n  %s" % d
   284         return s
   286     def __repr__(self):
   287         return self.__str__()
   293 def readSLHA(spcstr, ignorenobr=False):
   294     """
   295     Read an SLHA definition from a string, returning dictionaries of blocks and
   296     decays.
   298     If the ignorenobr parameter is True, do not store decay entries with a
   299     branching ratio of zero.
   300     """
   301     blocks = _mkdict()
   302     decays = _mkdict()
   303     #
   304     import re
   305     currentblock = None
   306     currentdecay = None
   307     for line in spcstr.splitlines():
   308         ## Handle (ignore) comment lines
   309         # TODO: Store block/entry comments
   310         if line.startswith("#"):
   311             continue
   312         if "#" in line:
   313             line = line[:line.index("#")]
   315         ## Handle BLOCK/DECAY start lines
   316         if line.upper().startswith("BLOCK"):
   317             #print line
   318             match = re.match(r"BLOCK\s+(\w+)(\s+Q\s*=\s*.+)?", line.upper())
   319             if not match:
   320                 continue
   321             blockname = match.group(1)
   322             qstr = match.group(2)
   323             if qstr is not None:
   324                 qstr = qstr[qstr.find("=")+1:].strip()
   325             currentblock = blockname
   326             currentdecay = None
   327             blocks[blockname] = Block(blockname, q=qstr)
   328         elif line.upper().startswith("DECAY"):
   329             match = re.match(r"DECAY\s+(-?\d+)\s+([\d\.E+-]+|NAN).*", line.upper())
   330             if not match:
   331                 continue
   332             pdgid = int(match.group(1))
   333             width = float(match.group(2)) if match.group(2) != "NAN" else None
   334             currentblock = "DECAY"
   335             currentdecay = pdgid
   336             decays[pdgid] = Particle(pdgid, width)
   337         else:
   338             ## In-block line
   339             if currentblock is not None:
   340                 items = line.split()
   341                 if len(items) < 1:
   342                     continue
   343                 if currentblock != "DECAY":
   344                     blocks[currentblock].add_entry(items)
   345                 else:
   346                     br = float(items[0]) if items[0].upper() != "NAN" else None
   347                     nda = int(items[1])
   348                     ids = map(int, items[2:])
   349                     if br > 0.0 or not ignorenobr: # br == None is < 0
   350                         decays[currentdecay].add_decay(br, nda, ids)
   352     ## Try to populate Particle masses from the MASS block
   353     # print blocks.keys()
   354     try:
   355         for pid in blocks["MASS"].keys():
   356             if decays.has_key(pid):
   357                 decays[pid].mass = blocks["MASS"][pid]
   358     except:
   359         raise ParseError("No MASS block found: cannot populate particle masses")
   361     return blocks, decays
   366 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   369 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
   370     """
   371     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   372     """
   373     # TODO: Pay attention to space-padding and minus signs for column alignment
   374     fmte = "%." + str(precision) + "e"
   375     sep = "   "
   376     blockstrs = []
   377     ## Blocks
   378     for bname, b in blocks.iteritems():
   379         namestr = b.name
   380         if b.q is not None:
   381             namestr += (" Q= " + fmte) % float(b.q)
   382         blockstr = "BLOCK %s\n" % namestr
   383         entrystrs = []
   384         for k, v in b.items():
   385             entrystr = ""
   386             if type(k) == tuple:
   387                 entrystr += sep.join(_autostr(i) for i in k)
   388             elif k is not None:
   389                 entrystr += _autostr(k)
   390             entrystr += sep + _autostr(v) # TODO: apply precision formatting for floats
   391             entrystrs.append(entrystr)
   392         blockstr += "\n".join(entrystrs)
   393         blockstrs.append(blockstr)
   394         ##
   395     ## Decays
   396     for pid, particle in decays.iteritems():
   397         blockstr = ("DECAY %d " + fmte + "\n") % (particle.pid, particle.totalwidth or -1)
   398         decaystrs = []
   399         for d in particle.decays:
   400             if d.br > 0.0 or not ignorenobr:
   401                 products_str = sep.join(map(str, d.ids))
   402                 decaystr = sep + (fmte % d.br) + sep + ("%d" % len(d.ids)) + sep + products_str
   403                 decaystrs.append(decaystr)
   404         blockstr += "\n".join(decaystrs)
   405         blockstrs.append(blockstr)
   406     ## Total result
   407     return "\n\n".join(blockstrs)
   411 ###############################################################################
   412 ## PDG <-> HERWIG particle ID code translations for ISAWIG handling
   414 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
   415 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   416 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   417 _HERWIGID2PDGID = {}
   418 _HERWIGID2PDGID[7]   = -1
   419 _HERWIGID2PDGID[8]   = -2
   420 _HERWIGID2PDGID[9]   = -3
   421 _HERWIGID2PDGID[10]  = -4
   422 _HERWIGID2PDGID[11]  = -5
   423 _HERWIGID2PDGID[12]  = -6
   424 _HERWIGID2PDGID[13]  =  21
   425 _HERWIGID2PDGID[59]  =  22
   426 _HERWIGID2PDGID[121] =  11
   427 _HERWIGID2PDGID[122] =  12
   428 _HERWIGID2PDGID[123] =  13
   429 _HERWIGID2PDGID[124] =  14
   430 _HERWIGID2PDGID[125] =  15
   431 _HERWIGID2PDGID[126] =  16
   432 _HERWIGID2PDGID[127] = -11
   433 _HERWIGID2PDGID[128] = -12
   434 _HERWIGID2PDGID[129] = -13
   435 _HERWIGID2PDGID[130] = -14
   436 _HERWIGID2PDGID[131] = -15
   437 _HERWIGID2PDGID[132] = -16
   438 _HERWIGID2PDGID[198] =  24 # W+
   439 _HERWIGID2PDGID[199] = -24 # W-
   440 _HERWIGID2PDGID[200] =  23 # Z0
   441 _HERWIGID2PDGID[201] =  25 ## SM HIGGS
   442 _HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   443 _HERWIGID2PDGID[204] =  35 ## HIGGSH0
   444 _HERWIGID2PDGID[205] =  36 ## HIGGSA0
   445 _HERWIGID2PDGID[206] =  37 ## HIGGS+
   446 _HERWIGID2PDGID[207] = -37 ## HIGGS-
   447 _HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   448 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   449 _HERWIGID2PDGID[402] =  1000002 ## SSULBR
   450 _HERWIGID2PDGID[408] = -1000002 ## SSUL
   451 _HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   452 _HERWIGID2PDGID[409] = -1000003 ## SSSL
   453 _HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   454 _HERWIGID2PDGID[410] = -1000004 ## SSCL
   455 _HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   456 _HERWIGID2PDGID[411] = -1000005 ## SSB1
   457 _HERWIGID2PDGID[406] =  1000006 ## SST1BR
   458 _HERWIGID2PDGID[412] = -1000006 ## SST1
   459 _HERWIGID2PDGID[413] =  2000001 ## SSDR
   460 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   461 _HERWIGID2PDGID[414] =  2000002 ## SSUR
   462 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
   463 _HERWIGID2PDGID[415] =  2000003 ## SSSR
   464 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   465 _HERWIGID2PDGID[416] =  2000004 ## SSCR
   466 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   467 _HERWIGID2PDGID[417] =  2000005 ## SSB2
   468 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   469 _HERWIGID2PDGID[418] =  2000006 ## SST2
   470 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
   471 _HERWIGID2PDGID[425] =  1000011 ## SSEL-
   472 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
   473 _HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   474 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   475 _HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   476 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   477 _HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   478 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   479 _HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   480 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   481 _HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   482 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   483 _HERWIGID2PDGID[437] =  2000011 ## SSEL-
   484 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
   485 _HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   486 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   487 _HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   488 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   489 _HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   490 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   491 _HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   492 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   493 _HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   494 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   495 _HERWIGID2PDGID[449] =  1000021 ## GLUINO
   496 _HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   497 _HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   498 _HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   499 _HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   500 _HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   501 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   502 _HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   503 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   504 _HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   506 def herwigid2pdgid(hwid):
   507     """
   508     Convert a particle ID code in the HERWIG internal IDHW format (as used by
   509     ISAWIG) into its equivalent in the standard PDG ID code definition.
   510     """
   511     return _HERWIGID2PDGID.get(hwid, hwid)
   514 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   515 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   516 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   517 _PDGID2HERWIGID = {}
   518 _PDGID2HERWIGID[      -1] = 7
   519 _PDGID2HERWIGID[      -2] = 8
   520 _PDGID2HERWIGID[      -3] = 9
   521 _PDGID2HERWIGID[      -4] = 10
   522 _PDGID2HERWIGID[      -5] = 11
   523 _PDGID2HERWIGID[      -6] = 12
   524 _PDGID2HERWIGID[      21] = 13
   525 _PDGID2HERWIGID[      22] = 59
   526 _PDGID2HERWIGID[      11] = 121
   527 _PDGID2HERWIGID[      12] = 122
   528 _PDGID2HERWIGID[      13] = 123
   529 _PDGID2HERWIGID[      14] = 124
   530 _PDGID2HERWIGID[      15] = 125
   531 _PDGID2HERWIGID[      16] = 126
   532 _PDGID2HERWIGID[     -11] = 127
   533 _PDGID2HERWIGID[     -12] = 128
   534 _PDGID2HERWIGID[     -13] = 129
   535 _PDGID2HERWIGID[     -14] = 130
   536 _PDGID2HERWIGID[     -15] = 131
   537 _PDGID2HERWIGID[     -16] = 132
   538 _PDGID2HERWIGID[      24] = 198 ## W+
   539 _PDGID2HERWIGID[     -24] = 199 ## W-
   540 _PDGID2HERWIGID[      23] = 200 ## Z
   541 _PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
   542 _PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   543 _PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   544 _PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   545 _PDGID2HERWIGID[      37] = 206 ## HIGGS+
   546 _PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   547 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   548 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   549 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   550 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
   551 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   552 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
   553 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   554 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
   555 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   556 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
   557 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   558 _PDGID2HERWIGID[-1000006] = 412 ## SST1
   559 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   560 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   561 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   562 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   563 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   564 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   565 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   566 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   567 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   568 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   569 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
   570 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   571 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   572 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   573 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   574 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   575 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   576 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   577 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   578 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   579 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   580 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   581 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   582 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   583 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   584 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   585 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   586 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   587 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   588 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   589 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   590 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   591 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   592 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   593 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   594 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   595 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   596 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   597 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   598 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   599 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   600 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   601 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   602 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   603 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   604 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   606 def pdgid2herwigid(pdgid):
   607     """
   608     Convert a particle ID code in the standard PDG ID code definition into
   609     its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
   610     """
   611     return _PDGID2HERWIGID.get(pdgid, pdgid)
   614 ###############################################################################
   615 ## ISAWIG format reading/writing
   618 def readISAWIG(isastr, ignorenobr=False):
   619     """
   620     Read a spectrum definition from a string in the ISAWIG format, returning
   621     dictionaries of blocks and decays. While this is not an SLHA format, it is
   622     informally supported as a useful mechanism for converting ISAWIG spectra to
   623     SLHA.
   625     ISAWIG parsing based on the HERWIG SUSY specification format, from
   626     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   628     If the ignorenobr parameter is True, do not store decay entries with a
   629     branching ratio of zero.
   630     """
   632     blocks = _mkdict()
   633     decays = _mkdict()
   634     LINES = isastr.splitlines()
   636     def getnextvalidline():
   637         while LINES:
   638             s = LINES.pop(0).strip()
   639             ## Return None if EOF reached
   640             if len(s) == 0:
   641                 continue
   642             ## Strip comments
   643             if "#" in s:
   644                 s = s[:s.index("#")].strip()
   645             ## Return if non-empty
   646             if len(s) > 0:
   647                 return s
   649     def getnextvalidlineitems():
   650         return map(_autotype, getnextvalidline().split())
   652     ## Populate MASS block and create decaying particle objects
   653     masses = Block("MASS")
   654     numentries = int(getnextvalidline())
   655     for i in xrange(numentries):
   656         hwid, mass, lifetime = getnextvalidlineitems()
   657         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   658         pdgid = herwigid2pdgid(hwid)
   659         masses[pdgid] = mass
   660         decays[pdgid] = Particle(pdgid, width, mass)
   661         #print pdgid, mass, width
   662     blocks["MASS"] = masses
   664     ## Populate decays
   665     for n in xrange(numentries):
   666         numdecays = int(getnextvalidline())
   667         for d in xrange(numdecays):
   668             #print n, numentries-1, d, numdecays-1
   669             decayitems = getnextvalidlineitems()
   670             hwid = decayitems[0]
   671             pdgid = herwigid2pdgid(hwid)
   672             br = decayitems[1]
   673             nme = decayitems[2]
   674             daughter_hwids = decayitems[3:]
   675             daughter_pdgids = []
   676             for hw in daughter_hwids:
   677                 if hw != 0:
   678                     daughter_pdgids.append(herwigid2pdgid(hw))
   679             if not decays.has_key(pdgid):
   680                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   681                 decays[pdgid] = Particle(pdgid)
   682             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   685     ## Now the SUSY parameters
   686     TANB, ALPHAH = getnextvalidlineitems()
   687     blocks["MINPAR"] = Block("MINPAR")
   688     blocks["MINPAR"][3] = TANB
   689     blocks["ALPHA"] = Block("ALPHA")
   690     blocks["ALPHA"][None] = ALPHAH
   691     #
   692     ## Neutralino mixing matrix
   693     blocks["NMIX"] = Block("NMIX")
   694     for i in xrange(1, 5):
   695         nmix_i = getnextvalidlineitems()
   696         for j, v in enumerate(nmix_i):
   697             blocks["NMIX"][i, j+1] = v
   698     #
   699     ## Chargino mixing matrices V and U
   700     blocks["VMIX"] = Block("VMIX")
   701     vmix = getnextvalidlineitems()
   702     blocks["VMIX"][1, 1] = vmix[0]
   703     blocks["VMIX"][1, 2] = vmix[1]
   704     blocks["VMIX"][2, 1] = vmix[2]
   705     blocks["VMIX"][2, 2] = vmix[3]
   706     blocks["UMIX"] = Block("UMIX")
   707     umix = getnextvalidlineitems()
   708     blocks["UMIX"][1, 1] = umix[0]
   709     blocks["UMIX"][1, 2] = umix[1]
   710     blocks["UMIX"][2, 1] = umix[2]
   711     blocks["UMIX"][2, 2] = umix[3]
   712     #
   713     THETAT, THETAB, THETAL = getnextvalidlineitems()
   714     import math
   715     blocks["STOPMIX"] = Block("STOPMIX")
   716     blocks["STOPMIX"][1, 1] =  math.cos(THETAT)
   717     blocks["STOPMIX"][1, 2] = -math.sin(THETAT)
   718     blocks["STOPMIX"][2, 1] =  math.sin(THETAT)
   719     blocks["STOPMIX"][2, 2] =  math.cos(THETAT)
   720     blocks["SBOTMIX"] = Block("SBOTMIX")
   721     blocks["SBOTMIX"][1, 1] =  math.cos(THETAB)
   722     blocks["SBOTMIX"][1, 2] = -math.sin(THETAB)
   723     blocks["SBOTMIX"][2, 1] =  math.sin(THETAB)
   724     blocks["SBOTMIX"][2, 2] =  math.cos(THETAB)
   725     blocks["STAUMIX"] = Block("STAUMIX")
   726     blocks["STAUMIX"][1, 1] =  math.cos(THETAL)
   727     blocks["STAUMIX"][1, 2] = -math.sin(THETAL)
   728     blocks["STAUMIX"][2, 1] =  math.sin(THETAL)
   729     blocks["STAUMIX"][2, 2] =  math.cos(THETAL)
   730     #
   731     ATSS, ABSS, ALSS = getnextvalidlineitems()
   732     blocks["AU"] = Block("AU")
   733     blocks["AU"][3, 3] = ATSS))
   734     blocks["AD"] = Block("AD")
   735     blocks["AD"][3, 3] = ABSS))
   736     blocks["AE"] = Block("AE")
   737     blocks["AE"][3, 3] = ALSS))
   738     #
   739     MUSS = getnextvalidlineitems()[0]
   740     blocks["MINPAR"][4] = MUSS
   741     #
   743     # TODO: Parse RPV boolean and couplings into SLHA2 blocks
   745     return blocks, decays
   748 def writeISAJET(blocks, decays, outname, ignorenobr=False, precision=8):
   749     """
   750     Return a SUSY spectrum definition in the format required for input by ISAJET,
   751     as a string, from the supplied blocks and decays dicts.
   753     The outname parameter specifies the desired output filename from ISAJET: this
   754     will appear in the first line of the return value.
   756     If the ignorenobr parameter is True, do not write decay entries with a
   757     branching ratio of zero.
   758     """
   759     fmte = "%." + str(precision) + "e"
   761     masses = blocks["MASS"]
   763     ## Init output string
   764     out = ""
   766     ## First line is the output name
   767     out += "'%s'" % outname + "\n"
   769     ## Next the top mass
   770     out += fmte % masses[6] + "\n"
   772     ## Next the top mass
   773     out += fmte % masses[6] + "\n"
   775     ## mSUGRA parameters (one line)
   776     # e.g. 1273.78,713.286,804.721,4.82337
   778     ## Masses and trilinear couplings (3 lines)
   779     # e.g. 1163.14,1114.15,1118.99,374.664,209.593
   780     # e.g. 1069.54,1112.7,919.908,374.556,209.381,-972.817,-326.745,-406.494
   781     # e.g. 1163.14,1114.15,1118.99,374.712,210.328
   783     ## RPV couplings (?? lines)
   784     # e.g. 232.615,445.477
   786     ## Etc ???!!!
   787     # e.g. /
   788     # e.g. n
   789     # e.g. y
   790     # e.g. y
   791     # e.g. 0.047441 3.80202e-23 0 0 0 2.17356e-22 0 0 5.23773e-09
   792     # e.g. y
   793     # 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
   794     # e.g. n
   795     # e.g. 'susy_RPV_stau_BC1scan_m560_tanb05.txt'
   797     return out
   800 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
   801     """
   802     Return a SUSY spectrum definition in the format produced by ISAWIG for inut to HERWIG
   803     as a string, from the supplied SLHA blocks and decays dicts.
   805     ISAWIG parsing based on the HERWIG SUSY specification format, from
   806     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   808     If the ignorenobr parameter is True, do not write decay entries with a
   809     branching ratio of zero.
   810     """
   811     fmte = "%." + str(precision) + "e"
   813     masses = blocks["MASS"]
   815     ## Init output string
   816     out = ""
   818     ## First write out masses section:
   819     ##   Number of SUSY + top particles
   820     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   821     ##   repeated for each particle
   822     ## IDHW is the HERWIG identity code.
   823     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   824     massout = ""
   825     for pid in masses.keys():
   826         lifetime = -1
   827         try:
   828             width = decays[pid].totalwidth
   829             if width and width > 0:
   830                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   831         except:
   832             pass
   833         massout += ("%d " + fmte + " " + fmte + "\n") % (pdgid2herwigid(pid), masses[pid], lifetime)
   834     out += "%d\n" % massout.count("\n")
   835     out += massout
   837     assert(len(masses) == len(decays))
   839     ## Next each particles decay modes together with their branching ratios and matrix element codes
   840     ##   Number of decay modes for a given particle (IDK)
   841     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   842     ##     repeated for each mode.
   843     ##   Repeated for each particle.
   844     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   845     ## the decay mode. NME is a code for the matrix element to be used, either from the
   846     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   847     for i, pid in enumerate(decays.keys()):
   848         # if not decays.has_key(pid):
   849         #     continue
   850         hwid = pdgid2herwigid(pid)
   851         decayout = ""
   852         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   853         for i_d, d in enumerate(decays[pid].decays):
   854             ## Skip decay if it has no branching ratio
   855             if ignorenobr and d.br == 0:
   856                 continue
   858             ## Identify decay matrix element to use
   859             ## From std HW docs, or from this pair:
   860             ## Two new matrix element codes have been added for these new decays:
   861             ##    NME =	200 	3 body top quark via charged Higgs
   862             ##    	300 	3 body R-parity violating gaugino and gluino decays
   863             nme = 0
   864             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   865             if abs(pid) in (6, 12):
   866                 nme = 100
   867             ## Extra SUSY MEs
   868             if len(d.ids) == 3:
   869                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   870                 pass
   871             decayout += ("%d " + fmte + " %d ") % (hwid, d.br, nme)
   873             def is_quark(pid):
   874                 return (abs(pid) in range(1, 7))
   876             def is_lepton(pid):
   877                 return (abs(pid) in range(11, 17))
   879             def is_squark(pid):
   880                 if abs(pid) in range(1000001, 1000007):
   881                     return True
   882                 if abs(pid) in range(2000001, 2000007):
   883                     return True
   884                 return False
   886             def is_slepton(pid):
   887                 if abs(pid) in range(1000011, 1000017):
   888                     return True
   889                 if abs(pid) in range(2000011, 2000016, 2):
   890                     return True
   891                 return False
   893             def is_gaugino(pid):
   894                 if abs(pid) in range(1000022, 1000026):
   895                     return True
   896                 if abs(pid) in (1000035, 1000037):
   897                     return True
   898                 return False
   900             def is_susy(pid):
   901                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   903             absids = map(abs, d.ids)
   905             ## Order decay products as required by HERWIG
   906             ## Top
   907             if abs(pid) == 6:
   908                 def cmp_bottomlast(a, b):
   909                     """Comparison function which always puts b/bbar last"""
   910                     if abs(a) == 5:
   911                         return True
   912                     if abs(b) == 5:
   913                         return False
   914                     return cmp(a, b)
   915                 if len(absids) == 2:
   916                     ## 2 body mode, to Higgs: Higgs; Bottom
   917                     if (25 in absids or 26 in absids) and 5 in absids:
   918                         d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   919                 elif len(absids) == 3:
   920                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   921                     if 37 in absids or 23 in absids:
   922                         d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   923             ## Gluino
   924             elif abs(pid) == 1000021:
   925                 if len(absids) == 2:
   926                     ## 2 body mode
   927                     ## without gluon: any order
   928                     ## with gluon: gluon; colour neutral
   929                     if 21 in absids:
   930                         def cmp_gluonfirst(a, b):
   931                             """Comparison function which always puts gluon first"""
   932                             if a == 21:
   933                                 return False
   934                             if b == 21:
   935                                 return True
   936                             return cmp(a, b)
   937                         d.ids = sorted(d.ids, cmp=cmp_gluonfirst)
   938                 elif len(absids) == 3:
   939                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   940                     def cmp_quarkslast(a, b):
   941                         """Comparison function which always puts quarks last"""
   942                         if is_quark(a):
   943                             return True
   944                         if is_quark(b):
   945                             return False
   946                         return cmp(a, b)
   947                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
   948             ## Squark/Slepton
   949             elif is_squark(pid) or is_slepton(pid):
   950                 def cmp_susy_quark_lepton(a, b):
   951                     if is_susy(a):
   952                         return False
   953                     if is_susy(b):
   954                         return True
   955                     if is_quark(a):
   956                         return False
   957                     if is_quark(b):
   958                         return True
   959                     return cmp(a, b)
   960                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   961                 ##                                                      Gluino       lepton
   962                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   963                 ## Squark
   964                 ##   2 body modes: Lepton Number Violated               quark     lepton
   965                 ##                 Baryon number violated               quark     quark
   966                 ## Slepton
   967                 ##   2 body modes: Lepton Number Violated               q or qbar
   968                 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
   969             ## Higgs
   970             elif pid in (25, 26):
   971                 # TODO: Includes SUSY Higgses?
   972                 ## Higgs
   973                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   974                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   975                 ##   3 body modes:                                      colour neutral       q or qbar
   976                 if len(absids) == 3:
   977                     def cmp_quarkslast(a, b):
   978                         """Comparison function which always puts quarks last"""
   979                         if is_quark(a):
   980                             return True
   981                         if is_quark(b):
   982                             return False
   983                         return cmp(a, b)
   984                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
   985             elif is_gaugino(pid):
   986                 # TODO: Is there actually anything to do here?
   987                 ## Gaugino
   988                 ##   2 body modes: Squark-quark                         q or sq
   989                 ##                 Slepton-lepton                       l or sl
   990                 ##
   991                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
   992                 ##                                                                           l or lbar
   993                 if len(absids) == 3:
   994                     def cmp_quarkslast(a, b):
   995                         """Comparison function which always puts quarks last"""
   996                         if is_quark(a):
   997                             return True
   998                         if is_quark(b):
   999                             return False
  1000                         return cmp(a, b)
  1001                     d.ids = sorted(d.ids, cmp=cmp_quarkslast)
  1003             # TODO: Gaugino/Gluino
  1004             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
  1006             ## Pad out IDs list with zeros
  1007             ids = [0,0,0,0,0]
  1008             for i, pid in enumerate(d.ids):
  1009                 ids[i] = pid
  1010             ids = map(str, ids)
  1011             decayout += " ".join(ids) + "\n"
  1012         decayout = "%d\n" % decayout.count("\n") + decayout
  1013         out += decayout
  1015     ## Now the SUSY parameters
  1016     ## TANB, ALPHAH:
  1017     out += (fmte + " " + fmte + "\n") % (blocks["MINPAR"].value(3), blocks["ALPHA"].value())
  1018     ## Neutralino mixing matrix
  1019     nmix = blocks["NMIX"]
  1020     for i in xrange(1, 5):
  1021         out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (nmix[i,1], nmix[i,2], nmix[i,3], nmix[i,4])
  1022     ## Chargino mixing matrices V and U
  1023     vmix = blocks["VMIX"]
  1024     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (vmix[1,1], vmix[1,2], vmix[2,1], vmix[2,2])
  1025     umix = blocks["UMIX"]
  1026     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (umix[1,1], umix[1,2], umix[2,1], umix[2,2])
  1027     ## THETAT,THETAB,THETAL
  1028     import math
  1029     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (math.acos(blocks["STOPMIX"][1,1]),
  1030                                                             math.acos(blocks["SBOTMIX"][1,1]),
  1031                                                             math.acos(blocks["STAUMIX"][1,1]))
  1032     ## ATSS,ABSS,ALSS
  1033     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (blocks["AU"][3,3],
  1034                                                             blocks["AD"][3,3],
  1035                                                             blocks["AE"][3,3])
  1036     ## MUSS == sign(mu)
  1037     out += "%f\n" % blocks["MINPAR"][4]
  1039     ## RPV SUSY
  1040     isRPV = False
  1041     out += "%d\n" % isRPV
  1042     # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
  1043     # TODO: Get the index orderings right
  1044     # if isRPV: ...
  1046     return out
  1049 ###############################################################################
  1050 ## File-level functions
  1053 def readSLHAFile(spcfilename, **kwargs):
  1054     """
  1055     Read an SLHA file, returning dictionaries of blocks and decays.
  1057     Other keyword parameters are passed to readSLHA.
  1058     """
  1059     f = open(spcfilename, "r")
  1060     rtn = readSLHA(f.read(), kwargs)
  1061     f.close()
  1062     return rtn
  1065 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
  1066     """
  1067     Write an SLHA file from the supplied blocks and decays dicts.
  1069     Other keyword parameters are passed to writeSLHA.
  1070     """
  1071     f = open(spcfilename, "w")
  1072     f.write(writeSLHA(blocks, decays, kwargs))
  1073     f.close()
  1076 def readISAWIGFile(isafilename, **kwargs):
  1077     """
  1078     Read a spectrum definition from a file in the ISAWIG format, returning
  1079     dictionaries of blocks and decays. While this is not an SLHA format, it is
  1080     informally supported as a useful mechanism for converting ISAWIG spectra to
  1081     SLHA.
  1083     Other keyword parameters are passed to readSLHA.
  1084     """
  1085     f = open(isafilename, "r")
  1086     rtn = readISAWIG(f.read(), kwargs)
  1087     f.close()
  1088     return rtn
  1091 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
  1092     """
  1093     Write an ISAWIG file from the supplied blocks and decays dicts.
  1095     Other keyword parameters are passed to writeISAWIG.
  1096     """
  1097     f = open(isafilename, "w")
  1098     f.write(writeISAWIG(blocks, decays, kwargs))
  1099     f.close()
  1102 def writeISAJETFile(isafilename, blocks, decays, **kwargs):
  1103     """
  1104     Write an ISAJET file from the supplied blocks and decays dicts (see writeISAJET).
  1106     Other keyword parameters are passed to writeISAJET.
  1107     """
  1108     f = open(isafilename, "w")
  1109     f.write(writeISAWIG(blocks, decays, kwargs))
  1110     f.close()
  1114 ###############################################################################
  1115 ## Main function for module testing
  1118 if __name__ == "__main__":
  1119     import sys
  1120     for a in sys.argv[1:]:
  1121         if a.endswith(".isa"):
  1122             blocks, decays = readISAWIGFile(a)
  1123         else:
  1124             blocks, decays = readSLHAFile(a)
  1126         for bname, b in sorted(blocks.iteritems()):
  1127             print b
  1128             print
  1130         print blocks.keys()
  1132         print blocks["MASS"][25]
  1133         print
  1135         for p in sorted(decays.values()):
  1136             print p
  1137             print
  1139         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial