pyslha.py

Sun, 28 Apr 2013 22:13:19 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 22:13:19 +0200
changeset 203
ce90a0dace07
parent 202
d3d069f4549a
child 204
bef82eaef56e
permissions
-rw-r--r--

Adding Block.set_value(*args) and Block documentation.

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

mercurial