pyslha.py

Sun, 28 Apr 2013 18:38:46 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 18:38:46 +0200
changeset 194
b6709078072b
parent 193
40d024dac179
child 195
880af4ab57ba
permissions
-rw-r--r--

New single/multi-value design *seems* to be working

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

mercurial