pyslha.py

Fri, 26 Apr 2013 20:57:42 +0200

author
Andy Buckley <andy@insectnation.org>
date
Fri, 26 Apr 2013 20:57:42 +0200
changeset 192
453a523cba25
parent 189
60c73b489420
child 193
40d024dac179
permissions
-rw-r--r--

Converting entries storage to avoid recursion, by using tuples as dict keys if necessary. Adding direct entry access via the [] operator to the Block class's entries. Improved SLHA output formatting.

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

mercurial