pyslha.py

Thu, 25 Apr 2013 11:43:16 +0200

author
Andy Buckley <andy@insectnation.org>
date
Thu, 25 Apr 2013 11:43:16 +0200
changeset 188
e706e3b2a647
parent 185
6cae67cd3cc1
child 189
60c73b489420
permissions
-rw-r--r--

Noted some suggested TODOs

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

mercurial