pyslha.py

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

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

Tidy

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

mercurial