pyslha.py

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

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 22:24:20 +0200
changeset 207
2990422f1738
parent 205
536967068fe2
child 208
67e368e5f414
permissions
-rw-r--r--

Adding missing spectrum files

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

mercurial