pyslha.py

Mon, 29 Apr 2013 14:35:12 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 29 Apr 2013 14:35:12 +0200
changeset 210
0f4f5472b7d8
parent 208
67e368e5f414
child 211
91f559c01cf7
permissions
-rw-r--r--

Split writeSLHA implementation into distinct writeSLHA{Blocks,Decays}.

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

mercurial