pyslha.py

Sat, 03 Mar 2012 20:38:59 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sat, 03 Mar 2012 20:38:59 +0100
changeset 170
50af1eaa6dfb
parent 168
9736876969a4
child 171
741d03a41ba3
permissions
-rw-r--r--

Correct version

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

mercurial