pyslha.py

Sun, 28 Apr 2013 19:21:33 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 19:21:33 +0200
changeset 196
7caa9b877b18
parent 195
880af4ab57ba
child 197
fdcc9cd1e00d
permissions
-rw-r--r--

More tidying and tweaking

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

mercurial