pyslha.py

Sun, 28 Apr 2013 21:20:55 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 21:20:55 +0200
changeset 201
3ac555efdbba
parent 200
651b0fac2163
child 202
d3d069f4549a
permissions
-rw-r--r--

Removing more indirections via .entries

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

mercurial