pyslha.py

Thu, 13 Oct 2011 01:11:47 +0200

author
Andy Buckley <andy@insectnation.org>
date
Thu, 13 Oct 2011 01:11:47 +0200
changeset 160
8acc44158b3f
parent 158
5cb1da0bddac
child 161
9306d1ad087f
permissions
-rw-r--r--

Write out the RPV = 0 term... need to add all the RPV couplings, too

     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  * Split writeSLHA into writeSLHA{Blocks,Decays}
    22  * Identify HERWIG decay matrix element to use in ISAWIG
    23  * Handle RPV SUSY in ISAWIG
    24 """
    26 __author__ = "Andy Buckley <andy.buckley@cern.ch"
    27 __version__ = "1.2.9"
    30 def _autotype(var):
    31     """Automatically convert strings to numerical types if possible."""
    32     if type(var) is not str:
    33         return var
    34     if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
    35         return int(var)
    36     try:
    37         f = float(var)
    38         return f
    39     except ValueError:
    40         return var
    42 def _autostr(var, precision=8):
    43     """Automatically numerical types to the right sort of string."""
    44     if type(var) is float:
    45         return ("%." + str(precision) + "e") % var
    46     return str(var)
    50 class Block(object):
    51     """
    52     Object representation of any BLOCK elements read from the SLHA file.  Blocks
    53     have a name, may have an associated Q value, and then a collection of data
    54     entries, stored as a recursive dictionary. Types in the dictionary are
    55     numeric (int or float) when a cast from the string in the file has been
    56     possible.
    57     """
    58     def __init__(self, name, q=None):
    59         self.name = name
    60         self.entries = {}
    61         self.q = _autotype(q)
    63     def add_entry(self, entry):
    64         #print entry
    65         nextparent = self.entries
    66         if len(entry) < 2:
    67             raise Exception("Block entries must be at least a 2-tuple")
    68         #print "in", entry
    69         entry = map(_autotype, entry)
    70         #print "out", entry
    71         for e in entry[:-2]:
    72             if e is not entry[-1]:
    73                 nextparent = nextparent.setdefault(e, {})
    74         nextparent[entry[-2]] = entry[-1]
    75         #print self.entries
    77     def __cmp__(self, other):
    78         return cmp(self.name, other.name)
    80     def __str__(self):
    81         s = self.name
    82         if self.q is not None:
    83             s += " (Q=%s)" % self.q
    84         s += "\n"
    85         s += str(self.entries)
    86         return s
    88     def __repr__(self):
    89         return self.__str__()
    92 class Decay(object):
    93     """
    94     Object representing a decay entry on a particle decribed by the SLHA file.
    95     'Decay' objects are not a direct representation of a DECAY block in an SLHA
    96     file... that role, somewhat confusingly, is taken by the Particle class.
    98     Decay objects have three properties: a branching ratio, br, an nda number
    99     (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
   100     decay occurs. The PDG ID of the particle whose decay this represents may
   101     also be stored, but this is normally known via the Particle in which the
   102     decay is stored.
   103     """
   104     def __init__(self, br, nda, ids, parentid=None):
   105         self.parentid = parentid
   106         self.br = br
   107         self.nda = nda
   108         self.ids = ids
   109         assert(self.nda == len(self.ids))
   111     def __cmp__(self, other):
   112         return cmp(other.br, self.br)
   114     def __str__(self):
   115         return "%.8e %s" % (self.br, self.ids)
   117     def __repr__(self):
   118         return self.__str__()
   121 class Particle(object):
   122     """
   123     Representation of a single, specific particle, decay block from an SLHA
   124     file.  These objects are not themselves called 'Decay', since that concept
   125     applies more naturally to the various decays found inside this
   126     object. Particle classes store the PDG ID (pid) of the particle being
   127     represented, and optionally the mass (mass) and total decay width
   128     (totalwidth) of that particle in the SLHA scenario. Masses may also be found
   129     via the MASS block, from which the Particle.mass property is filled, if at
   130     all. They also store a list of Decay objects (decays) which are probably the
   131     item of most interest.
   132     """
   133     def __init__(self, pid, totalwidth=None, mass=None):
   134         self.pid = pid
   135         self.totalwidth = totalwidth
   136         self.mass = mass
   137         self.decays = []
   139     def add_decay(self, br, nda, ids):
   140         self.decays.append(Decay(br, nda, ids))
   141         self.decays.sort()
   143     def __cmp__(self, other):
   144         if abs(self.pid) == abs(other.pid):
   145             return cmp(self.pid, other.pid)
   146         return cmp(abs(self.pid), abs(other.pid))
   148     def __str__(self):
   149         s = str(self.pid)
   150         if self.mass is not None:
   151             s += " : mass = %.8e GeV" % self.mass
   152         if self.totalwidth is not None:
   153             s += " : total width = %.8e GeV" % self.totalwidth
   154         for d in self.decays:
   155             if d.br > 0.0:
   156                 s += "\n  %s" % d
   157         return s
   159     def __repr__(self):
   160         return self.__str__()
   164 def readSLHAFile(spcfilename, **kwargs):
   165     """
   166     Read an SLHA file, returning dictionaries of blocks and decays.
   168     Other keyword parameters are passed to readSLHA.
   169     """
   170     f = open(spcfilename, "r")
   171     rtn = readSLHA(f.read(), kwargs)
   172     f.close()
   173     return rtn
   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 = {}
   185     decays = {}
   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+-]+).*", line.upper())
   212             if not match:
   213                 continue
   214             pdgid = int(match.group(1))
   215             width = float(match.group(2))
   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])
   234                     nda = int(items[1])
   235                     ids = map(int, items[2:])
   236                     if br > 0.0 or not ignorenobr:
   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 Exception("No MASS block found, from which to populate particle masses")
   248     return blocks, decays
   251 def readISAWIGFile(isafilename, **kwargs):
   252     """
   253     Read a spectrum definition from a file in the ISAWIG format, returning
   254     dictionaries of blocks and decays. While this is not an SLHA format, it is
   255     informally supported as a useful mechanism for converting ISAWIG spectra to
   256     SLHA.
   258     Other keyword parameters are passed to readSLHA.
   259     """
   260     f = open(isafilename, "r")
   261     rtn = readISAWIG(f.read(), kwargs)
   262     f.close()
   263     return rtn
   266 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
   267     """
   268     Write an SLHA file from the supplied blocks and decays dicts.
   270     Other keyword parameters are passed to writeSLHA.
   271     """
   272     f = open(spcfilename, "w")
   273     f.write(writeSLHA(blocks, decays, kwargs))
   274     f.close()
   277 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   280 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
   281     """
   282     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   283     """
   284     fmte = "%." + str(precision) + "e"
   286     sep = "   "
   287     out = ""
   288     def dict_hier_strs(d, s=""):
   289         if type(d) is dict:
   290             for k, v in sorted(d.iteritems()):
   291                 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
   292                     yield s2
   293         else:
   294             yield s + sep + _autostr(d)
   295     ## Blocks
   296     for bname, b in sorted(blocks.iteritems()):
   297         namestr = b.name
   298         if b.q is not None:
   299             namestr += (" Q= " + fmte) % float(b.q)
   300         out += "BLOCK %s\n" % namestr
   301         for s in dict_hier_strs(b.entries):
   302             out += sep + s + "\n"
   303         out += "\n"
   304     ## Decays
   305     for pid, particle in sorted(decays.iteritems()):
   306         out += ("DECAY %d " + fmte + "\n") % (particle.pid, particle.totalwidth or -1)
   307         for d in sorted(particle.decays):
   308             if d.br > 0.0 or not ignorenobr:
   309                 products_str = "   ".join(map(str, d.ids))
   310                 out += sep + fmte % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
   311         out += "\n"
   312     return out
   316 ###############################################################################
   317 ## ISAWIG handling
   319 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
   320 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   321 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   322 _HERWIGID2PDGID = {}
   323 _HERWIGID2PDGID[7]   = -1
   324 _HERWIGID2PDGID[8]   = -2
   325 _HERWIGID2PDGID[9]   = -3
   326 _HERWIGID2PDGID[10]  = -4
   327 _HERWIGID2PDGID[11]  = -5
   328 _HERWIGID2PDGID[12]  = -6
   329 _HERWIGID2PDGID[13]  =  21
   330 _HERWIGID2PDGID[59]  =  22
   331 _HERWIGID2PDGID[121] =  11
   332 _HERWIGID2PDGID[122] =  12
   333 _HERWIGID2PDGID[123] =  13
   334 _HERWIGID2PDGID[124] =  14
   335 _HERWIGID2PDGID[125] =  15
   336 _HERWIGID2PDGID[126] =  16
   337 _HERWIGID2PDGID[127] = -11
   338 _HERWIGID2PDGID[128] = -12
   339 _HERWIGID2PDGID[129] = -13
   340 _HERWIGID2PDGID[130] = -14
   341 _HERWIGID2PDGID[131] = -15
   342 _HERWIGID2PDGID[132] = -16
   343 _HERWIGID2PDGID[198] =  24 # W+
   344 _HERWIGID2PDGID[199] = -24 # W-
   345 _HERWIGID2PDGID[200] =  23 # Z0
   346 _HERWIGID2PDGID[201] =  25 ## SM HIGGS
   347 _HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   348 _HERWIGID2PDGID[204] =  35 ## HIGGSH0
   349 _HERWIGID2PDGID[205] =  36 ## HIGGSA0
   350 _HERWIGID2PDGID[206] =  37 ## HIGGS+
   351 _HERWIGID2PDGID[207] = -37 ## HIGGS-
   352 _HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   353 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   354 _HERWIGID2PDGID[402] =  1000002 ## SSULBR
   355 _HERWIGID2PDGID[408] = -1000002 ## SSUL
   356 _HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   357 _HERWIGID2PDGID[409] = -1000003 ## SSSL
   358 _HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   359 _HERWIGID2PDGID[410] = -1000004 ## SSCL
   360 _HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   361 _HERWIGID2PDGID[411] = -1000005 ## SSB1
   362 _HERWIGID2PDGID[406] =  1000006 ## SST1BR
   363 _HERWIGID2PDGID[412] = -1000006 ## SST1
   364 _HERWIGID2PDGID[413] =  2000001 ## SSDR
   365 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   366 _HERWIGID2PDGID[414] =  2000002 ## SSUR
   367 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
   368 _HERWIGID2PDGID[415] =  2000003 ## SSSR
   369 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   370 _HERWIGID2PDGID[416] =  2000004 ## SSCR
   371 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   372 _HERWIGID2PDGID[417] =  2000005 ## SSB2
   373 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   374 _HERWIGID2PDGID[418] =  2000006 ## SST2
   375 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
   376 _HERWIGID2PDGID[425] =  1000011 ## SSEL-
   377 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
   378 _HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   379 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   380 _HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   381 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   382 _HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   383 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   384 _HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   385 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   386 _HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   387 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   388 _HERWIGID2PDGID[437] =  2000011 ## SSEL-
   389 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
   390 _HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   391 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   392 _HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   393 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   394 _HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   395 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   396 _HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   397 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   398 _HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   399 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   400 _HERWIGID2PDGID[449] =  1000021 ## GLUINO
   401 _HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   402 _HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   403 _HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   404 _HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   405 _HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   406 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   407 _HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   408 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   409 _HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   411 def herwigid2pdgid(hwid):
   412     """
   413     Convert a particle ID code in the HERWIG internal IDHW format (as used by
   414     ISAWIG) into its equivalent in the standard PDG ID code definition.
   415     """
   416     return _HERWIGID2PDGID.get(hwid, hwid)
   419 def readISAWIG(isastr, ignorenobr=False):
   420     """
   421     Read a spectrum definition from a string in the ISAWIG format, returning
   422     dictionaries of blocks and decays. While this is not an SLHA format, it is
   423     informally supported as a useful mechanism for converting ISAWIG spectra to
   424     SLHA.
   426     ISAWIG parsing based on the HERWIG SUSY specification format, from
   427     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   429     If the ignorenobr parameter is True, do not store decay entries with a
   430     branching ratio of zero.
   431     """
   434     blocks = {}
   435     decays = {}
   436     LINES = isastr.splitlines()
   438     def getnextvalidline():
   439         while LINES:
   440             s = LINES.pop(0).strip()
   441             ## Return None if EOF reached
   442             if len(s) == 0:
   443                 continue
   444             ## Strip comments
   445             if "#" in s:
   446                 s = s[:s.index("#")].strip()
   447             ## Return if non-empty
   448             if len(s) > 0:
   449                 return s
   451     def getnextvalidlineitems():
   452         return map(_autotype, getnextvalidline().split())
   454     ## Populate MASS block and create decaying particle objects
   455     masses = Block("MASS")
   456     numentries = int(getnextvalidline())
   457     for i in xrange(numentries):
   458         hwid, mass, lifetime = getnextvalidlineitems()
   459         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   460         pdgid = herwigid2pdgid(hwid)
   461         masses.add_entry((pdgid, mass))
   462         decays[pdgid] = Particle(pdgid, width, mass)
   463         #print pdgid, mass, width
   464     blocks["MASS"] = masses
   466     ## Populate decays
   467     for n in xrange(numentries):
   468         numdecays = int(getnextvalidline())
   469         for d in xrange(numdecays):
   470             #print n, numentries-1, d, numdecays-1
   471             decayitems = getnextvalidlineitems()
   472             hwid = decayitems[0]
   473             pdgid = herwigid2pdgid(hwid)
   474             br = decayitems[1]
   475             nme = decayitems[2]
   476             daughter_hwids = decayitems[3:]
   477             daughter_pdgids = []
   478             for hw in daughter_hwids:
   479                 if hw != 0:
   480                     daughter_pdgids.append(herwigid2pdgid(hw))
   481             if not decays.has_key(pdgid):
   482                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   483                 decays[pdgid] = Particle(pdgid)
   484             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   487     ## Now the SUSY parameters
   488     TANB, ALPHAH = getnextvalidlineitems()
   489     blocks["MINPAR"] = Block("MINPAR")
   490     blocks["MINPAR"].add_entry((3, TANB))
   491     blocks["ALPHA"] = Block("ALPHA")
   492     blocks["ALPHA"].entries = ALPHAH
   493     #
   494     ## Neutralino mixing matrix
   495     blocks["NMIX"] = Block("NMIX")
   496     for i in xrange(1, 5):
   497         nmix_i = getnextvalidlineitems()
   498         for j, v in enumerate(nmix_i):
   499             blocks["NMIX"].add_entry((i, j+1, v))
   500     #
   501     ## Chargino mixing matrices V and U
   502     blocks["VMIX"] = Block("VMIX")
   503     vmix = getnextvalidlineitems()
   504     blocks["VMIX"].add_entry((1, 1, vmix[0]))
   505     blocks["VMIX"].add_entry((1, 2, vmix[1]))
   506     blocks["VMIX"].add_entry((2, 1, vmix[2]))
   507     blocks["VMIX"].add_entry((2, 2, vmix[3]))
   508     blocks["UMIX"] = Block("UMIX")
   509     umix = getnextvalidlineitems()
   510     blocks["UMIX"].add_entry((1, 1, umix[0]))
   511     blocks["UMIX"].add_entry((1, 2, umix[1]))
   512     blocks["UMIX"].add_entry((2, 1, umix[2]))
   513     blocks["UMIX"].add_entry((2, 2, umix[3]))
   514     #
   515     THETAT, THETAB, THETAL = getnextvalidlineitems()
   516     import math
   517     blocks["STOPMIX"] = Block("STOPMIX")
   518     blocks["STOPMIX"].add_entry((1, 1,  math.cos(THETAT)))
   519     blocks["STOPMIX"].add_entry((1, 2, -math.sin(THETAT)))
   520     blocks["STOPMIX"].add_entry((2, 1,  math.sin(THETAT)))
   521     blocks["STOPMIX"].add_entry((2, 2,  math.cos(THETAT)))
   522     blocks["SBOTMIX"] = Block("SBOTMIX")
   523     blocks["SBOTMIX"].add_entry((1, 1,  math.cos(THETAB)))
   524     blocks["SBOTMIX"].add_entry((1, 2, -math.sin(THETAB)))
   525     blocks["SBOTMIX"].add_entry((2, 1,  math.sin(THETAB)))
   526     blocks["SBOTMIX"].add_entry((2, 2,  math.cos(THETAB)))
   527     blocks["STAUMIX"] = Block("STAUMIX")
   528     blocks["STAUMIX"].add_entry((1, 1,  math.cos(THETAL)))
   529     blocks["STAUMIX"].add_entry((1, 2, -math.sin(THETAL)))
   530     blocks["STAUMIX"].add_entry((2, 1,  math.sin(THETAL)))
   531     blocks["STAUMIX"].add_entry((2, 2,  math.cos(THETAL)))
   532     #
   533     ATSS, ABSS, ALSS = getnextvalidlineitems()
   534     blocks["AU"] = Block("AU")
   535     blocks["AU"].add_entry((3, 3, ATSS))
   536     blocks["AD"] = Block("AD")
   537     blocks["AD"].add_entry((3, 3, ABSS))
   538     blocks["AE"] = Block("AE")
   539     blocks["AE"].add_entry((3, 3, ALSS))
   540     #
   541     MUSS = getnextvalidlineitems()[0]
   542     blocks["MINPAR"].add_entry((4, MUSS))
   543     #
   544     return blocks, decays
   547 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   548 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   549 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   550 _PDGID2HERWIGID = {}
   551 _PDGID2HERWIGID[      -1] = 7
   552 _PDGID2HERWIGID[      -2] = 8
   553 _PDGID2HERWIGID[      -3] = 9
   554 _PDGID2HERWIGID[      -4] = 10
   555 _PDGID2HERWIGID[      -5] = 11
   556 _PDGID2HERWIGID[      -6] = 12
   557 _PDGID2HERWIGID[      21] = 13
   558 _PDGID2HERWIGID[      22] = 59
   559 _PDGID2HERWIGID[      11] = 121
   560 _PDGID2HERWIGID[      12] = 122
   561 _PDGID2HERWIGID[      13] = 123
   562 _PDGID2HERWIGID[      14] = 124
   563 _PDGID2HERWIGID[      15] = 125
   564 _PDGID2HERWIGID[      16] = 126
   565 _PDGID2HERWIGID[     -11] = 127
   566 _PDGID2HERWIGID[     -12] = 128
   567 _PDGID2HERWIGID[     -13] = 129
   568 _PDGID2HERWIGID[     -14] = 130
   569 _PDGID2HERWIGID[     -15] = 131
   570 _PDGID2HERWIGID[     -16] = 132
   571 _PDGID2HERWIGID[      24] = 198 ## W+
   572 _PDGID2HERWIGID[     -24] = 199 ## W-
   573 _PDGID2HERWIGID[      23] = 200 ## Z
   574 _PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
   575 _PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   576 _PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   577 _PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   578 _PDGID2HERWIGID[      37] = 206 ## HIGGS+
   579 _PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   580 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   581 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   582 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   583 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
   584 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   585 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
   586 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   587 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
   588 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   589 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
   590 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   591 _PDGID2HERWIGID[-1000006] = 412 ## SST1
   592 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   593 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   594 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   595 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   596 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   597 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   598 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   599 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   600 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   601 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   602 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
   603 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   604 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   605 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   606 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   607 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   608 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   609 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   610 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   611 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   612 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   613 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   614 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   615 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   616 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   617 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   618 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   619 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   620 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   621 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   622 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   623 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   624 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   625 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   626 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   627 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   628 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   629 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   630 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   631 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   632 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   633 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   634 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   635 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   636 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   637 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   639 def pdgid2herwigid(pdgid):
   640     """
   641     Convert a particle ID code in the standard PDG ID code definition into
   642     its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
   643     """
   644     return _PDGID2HERWIGID.get(pdgid, pdgid)
   649 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
   650     """
   651     Write an ISAWIG file from the supplied blocks and decays dicts.
   653     Other keyword parameters are passed to writeISAWIG.
   655     TODO: Handle RPV SUSY
   656     """
   657     f = open(isafilename, "w")
   658     f.write(writeISAWIG(blocks, decays, kwargs))
   659     f.close()
   662 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
   663     """
   664     Return an ISAWIG definition as a string, from the supplied blocks and decays dicts.
   666     ISAWIG parsing based on the HERWIG SUSY specification format, from
   667     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   669     If the ignorenobr parameter is True, do not write decay entries with a
   670     branching ratio of zero.
   671     """
   672     fmte = "%." + str(precision) + "e"
   674     masses = blocks["MASS"].entries
   676     ## Init output string
   677     out = ""
   679     ## First write out masses section:
   680     ##   Number of SUSY + top particles
   681     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   682     ##   repeated for each particle
   683     ## IDHW is the HERWIG identity code.
   684     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   685     massout = ""
   686     for pid in masses.keys():
   687         lifetime = -1
   688         try:
   689             width = decays[pid].totalwidth
   690             if width and width > 0:
   691                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   692         except:
   693             pass
   694         massout += ("%d " + fmte + " " + fmte + "\n") % (pdgid2herwigid(pid), masses[pid], lifetime)
   695     out += "%d\n" % massout.count("\n")
   696     out += massout
   698     assert(len(masses) == len(decays))
   700     ## Next each particles decay modes together with their branching ratios and matrix element codes
   701     ##   Number of decay modes for a given particle (IDK)
   702     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   703     ##     repeated for each mode.
   704     ##   Repeated for each particle.
   705     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   706     ## the decay mode. NME is a code for the matrix element to be used, either from the
   707     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   708     for i, pid in enumerate(decays.keys()):
   709         # if not decays.has_key(pid):
   710         #     continue
   711         hwid = pdgid2herwigid(pid)
   712         decayout = ""
   713         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   714         for i_d, d in enumerate(decays[pid].decays):
   715             ## Skip decay if it has no branching ratio
   716             if ignorenobr and d.br == 0:
   717                 continue
   719             ## Identify decay matrix element to use
   720             ## From std HW docs, or from this pair:
   721             ## Two new matrix element codes have been added for these new decays:
   722             ##    NME =	200 	3 body top quark via charged Higgs
   723             ##    	300 	3 body R-parity violating gaugino and gluino decays
   724             nme = 0
   725             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   726             if abs(pid) in (6, 12):
   727                 nme = 100
   728             ## Extra SUSY MEs
   729             if len(d.ids) == 3:
   730                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   731                 pass
   732             decayout += "%d " + fmte + " %d " % (hwid, d.br, nme)
   734             def is_quark(pid):
   735                 return (abs(pid) in range(1, 7))
   737             def is_lepton(pid):
   738                 return (abs(pid) in range(11, 17))
   740             def is_squark(pid):
   741                 if abs(pid) in range(1000001, 1000007):
   742                     return True
   743                 if abs(pid) in range(2000001, 2000007):
   744                     return True
   745                 return False
   747             def is_slepton(pid):
   748                 if abs(pid) in range(1000011, 1000017):
   749                     return True
   750                 if abs(pid) in range(2000011, 2000016, 2):
   751                     return True
   752                 return False
   754             def is_gaugino(pid):
   755                 if abs(pid) in range(1000022, 1000026):
   756                     return True
   757                 if abs(pid) in (1000035, 1000037):
   758                     return True
   759                 return False
   761             def is_susy(pid):
   762                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   764             absids = map(abs, d.ids)
   766             ## Order decay products as required by HERWIG
   767             ## Top
   768             if abs(pid) == 6:
   769                 def cmp_bottomlast(a, b):
   770                     """Comparison function which always puts b/bbar last"""
   771                     if abs(a) == 5:
   772                         return True
   773                     if abs(b) == 5:
   774                         return False
   775                     return cmp(a, b)
   776                 if len(absids) == 2:
   777                     ## 2 body mode, to Higgs: Higgs; Bottom
   778                     if (25 in absids or 26 in absids) and 5 in absids:
   779                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   780                 elif len(absids) == 3:
   781                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   782                     if 37 in absids or 23 in absids:
   783                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   784             ## Gluino
   785             elif abs(pid) == 1000021:
   786                 if len(absids) == 2:
   787                     ## 2 body mode
   788                     ## without gluon: any order
   789                     ## with gluon: gluon; colour neutral
   790                     if 21 in absids:
   791                         def cmp_gluonfirst(a, b):
   792                             """Comparison function which always puts gluon first"""
   793                             if a == 21:
   794                                 return False
   795                             if b == 21:
   796                                 return True
   797                             return cmp(a, b)
   798                         d.ids = sorted(d.ids, key=cmp_gluonfirst)
   799                 elif len(absids) == 3:
   800                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   801                     def cmp_quarkslast(a, b):
   802                         """Comparison function which always puts quarks last"""
   803                         if is_quark(a):
   804                             return True
   805                         if is_quark(b):
   806                             return False
   807                         return cmp(a, b)
   808                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   809             ## Squark/Slepton
   810             elif is_squark(pid) or is_slepton(pid):
   811                 def cmp_susy_quark_lepton(a, b):
   812                     if is_susy(a):
   813                         return False
   814                     if is_susy(b):
   815                         return True
   816                     if is_quark(a):
   817                         return False
   818                     if is_quark(b):
   819                         return True
   820                     return cmp(a, b)
   821                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   822                 ##                                                      Gluino       lepton
   823                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   824                 ## Squark
   825                 ##   2 body modes: Lepton Number Violated               quark     lepton
   826                 ##                 Baryon number violated               quark     quark
   827                 ## Slepton
   828                 ##   2 body modes: Lepton Number Violated               q or qbar
   829                 d.ids = sorted(d.ids, key=cmp_bottomlast)
   830             ## Higgs
   831             elif pid in (25, 26):
   832                 # TODO: Includes SUSY Higgses?
   833                 ## Higgs
   834                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   835                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   836                 ##   3 body modes:                                      colour neutral       q or qbar
   837                 if len(absids) == 3:
   838                     def cmp_quarkslast(a, b):
   839                         """Comparison function which always puts quarks last"""
   840                         if is_quark(a):
   841                             return True
   842                         if is_quark(b):
   843                             return False
   844                         return cmp(a, b)
   845                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   846             elif is_gaugino(pid):
   847                 # TODO: Is there actually anything to do here?
   848                 ## Gaugino
   849                 ##   2 body modes: Squark-quark                         q or sq
   850                 ##                 Slepton-lepton                       l or sl
   851                 ##
   852                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
   853                 ##                                                                           l or lbar
   854                 if len(absids) == 3:
   855                     def cmp_quarkslast(a, b):
   856                         """Comparison function which always puts quarks last"""
   857                         if is_quark(a):
   858                             return True
   859                         if is_quark(b):
   860                             return False
   861                         return cmp(a, b)
   862                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   864             # TODO: Gaugino/Gluino
   865             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
   867             ## Pad out IDs list with zeros
   868             ids = [0,0,0,0,0]
   869             for i, pid in enumerate(d.ids):
   870                 ids[i] = pid
   871             ids = map(str, ids)
   872             decayout += " ".join(ids) + "\n"
   873         decayout = "%d\n" % decayout.count("\n") + decayout
   874         out += decayout
   876     ## Now the SUSY parameters
   877     ## TANB, ALPHAH:
   878     out += (fmte + " " + fmte + "\n") % (blocks["MINPAR"].entries[3], blocks["ALPHA"].entries)
   879     ## Neutralino mixing matrix
   880     nmix = blocks["NMIX"].entries
   881     for i in xrange(1, 5):
   882         out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (nmix[i][1], nmix[i][2], nmix[i][3], nmix[i][4])
   883     ## Chargino mixing matrices V and U
   884     vmix = blocks["VMIX"].entries
   885     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (vmix[1][1], vmix[1][2], vmix[2][1], vmix[2][2])
   886     umix = blocks["UMIX"].entries
   887     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (umix[1][1], umix[1][2], umix[2][1], umix[2][2])
   888     ## THETAT,THETAB,THETAL
   889     import math
   890     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (math.acos(blocks["STOPMIX"].entries[1][1]),
   891                                                             math.acos(blocks["SBOTMIX"].entries[1][1]),
   892                                                             math.acos(blocks["STAUMIX"].entries[1][1]))
   893     ## ATSS,ABSS,ALSS
   894     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (blocks["AU"].entries[3][3],
   895                                                             blocks["AD"].entries[3][3],
   896                                                             blocks["AE"].entries[3][3])
   897     ## MUSS == sign(mu)
   898     out += "%f\n" % blocks["MINPAR"].entries[4]
   900     ## RPV SUSY
   901     isRPV = False
   902     out += "%d\n" % isRPV
   903     # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
   904     # TODO: Get the index orderings right
   905     # if isRPV: ...
   907     return out
   911 if __name__ == "__main__":
   912     import sys
   913     for a in sys.argv[1:]:
   914         if a.endswith(".isa"):
   915             blocks, decays = readISAWIGFile(a)
   916         else:
   917             blocks, decays = readSLHAFile(a)
   919         for bname, b in sorted(blocks.iteritems()):
   920             print b
   921             print
   923         print blocks.keys()
   925         print blocks["MASS"].entries[25]
   926         print
   928         for p in sorted(decays.values()):
   929             print p
   930             print
   932         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial