pyslha.py

Sun, 28 Apr 2013 18:33:16 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 18:33:16 +0200
changeset 193
40d024dac179
parent 192
453a523cba25
child 194
b6709078072b
permissions
-rw-r--r--

Working out a better way to deal with single/multivalue blocks

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

mercurial