pyslha.py

Fri, 26 Apr 2013 01:00:15 +0200

author
Andy Buckley <andy@insectnation.org>
date
Fri, 26 Apr 2013 01:00:15 +0200
changeset 189
60c73b489420
parent 188
e706e3b2a647
child 192
453a523cba25
permissions
-rw-r--r--

Preserving the ordering of blocks, decays, and their their contents if possible, using ordered dicts. Version 1.5.0, since the behaviour has significantly altered.

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

mercurial