pyslha.py

Thu, 13 Oct 2011 17:10:35 +0200

author
Andy Buckley <andy@insectnation.org>
date
Thu, 13 Oct 2011 17:10:35 +0200
changeset 161
9306d1ad087f
parent 160
8acc44158b3f
child 162
e5b76d1917c0
permissions
-rw-r--r--

Restructuring the module a bit, and starting work on writing the ISAJET output format (as opposed to the ISAWIG output format) via two new writeISAJET and writeISAJETFile functions.

     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__()
   166 def readSLHA(spcstr, ignorenobr=False):
   167     """
   168     Read an SLHA definition from a string, returning dictionaries of blocks and
   169     decays.
   171     If the ignorenobr parameter is True, do not store decay entries with a
   172     branching ratio of zero.
   173     """
   174     blocks = {}
   175     decays = {}
   176     #
   177     import re
   178     currentblock = None
   179     currentdecay = None
   180     for line in spcstr.splitlines():
   181         ## Handle (ignore) comment lines
   182         if line.startswith("#"):
   183             continue
   184         if "#" in line:
   185             line = line[:line.index("#")]
   187         ## Handle BLOCK/DECAY start lines
   188         if line.upper().startswith("BLOCK"):
   189             #print line
   190             match = re.match(r"BLOCK\s+(\w+)(\s+Q\s*=\s*.+)?", line.upper())
   191             if not match:
   192                 continue
   193             blockname = match.group(1)
   194             qstr = match.group(2)
   195             if qstr is not None:
   196                 qstr = qstr[qstr.find("=")+1:].strip()
   197             currentblock = blockname
   198             currentdecay = None
   199             blocks[blockname] = Block(blockname, q=qstr)
   200         elif line.upper().startswith("DECAY"):
   201             match = re.match(r"DECAY\s+(\d+)\s+([\d\.E+-]+).*", line.upper())
   202             if not match:
   203                 continue
   204             pdgid = int(match.group(1))
   205             width = float(match.group(2))
   206             currentblock = "DECAY"
   207             currentdecay = pdgid
   208             decays[pdgid] = Particle(pdgid, width)
   209         else:
   210             ## In-block line
   211             if currentblock is not None:
   212                 items = line.split()
   213                 if len(items) < 1:
   214                     continue
   215                 if currentblock != "DECAY":
   216                     if len(items) < 2:
   217                         ## Treat the ALPHA block differently
   218                         blocks[currentblock].value = _autotype(items[0])
   219                         blocks[currentblock].entries = _autotype(items[0])
   220                     else:
   221                         blocks[currentblock].add_entry(items)
   222                 else:
   223                     br = float(items[0])
   224                     nda = int(items[1])
   225                     ids = map(int, items[2:])
   226                     if br > 0.0 or not ignorenobr:
   227                         decays[currentdecay].add_decay(br, nda, ids)
   229     ## Try to populate Particle masses from the MASS block
   230     # print blocks.keys()
   231     try:
   232         for pid in blocks["MASS"].entries.keys():
   233             if decays.has_key(pid):
   234                 decays[pid].mass = blocks["MASS"].entries[pid]
   235     except:
   236         raise Exception("No MASS block found, from which to populate particle masses")
   238     return blocks, decays
   243 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
   246 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
   247     """
   248     Return an SLHA definition as a string, from the supplied blocks and decays dicts.
   249     """
   250     fmte = "%." + str(precision) + "e"
   252     sep = "   "
   253     out = ""
   254     def dict_hier_strs(d, s=""):
   255         if type(d) is dict:
   256             for k, v in sorted(d.iteritems()):
   257                 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
   258                     yield s2
   259         else:
   260             yield s + sep + _autostr(d)
   261     ## Blocks
   262     for bname, b in sorted(blocks.iteritems()):
   263         namestr = b.name
   264         if b.q is not None:
   265             namestr += (" Q= " + fmte) % float(b.q)
   266         out += "BLOCK %s\n" % namestr
   267         for s in dict_hier_strs(b.entries):
   268             out += sep + s + "\n"
   269         out += "\n"
   270     ## Decays
   271     for pid, particle in sorted(decays.iteritems()):
   272         out += ("DECAY %d " + fmte + "\n") % (particle.pid, particle.totalwidth or -1)
   273         for d in sorted(particle.decays):
   274             if d.br > 0.0 or not ignorenobr:
   275                 products_str = "   ".join(map(str, d.ids))
   276                 out += sep + fmte % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
   277         out += "\n"
   278     return out
   282 ###############################################################################
   283 ## PDG <-> HERWIG particle ID code translations for ISAWIG handling
   285 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
   286 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   287 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   288 _HERWIGID2PDGID = {}
   289 _HERWIGID2PDGID[7]   = -1
   290 _HERWIGID2PDGID[8]   = -2
   291 _HERWIGID2PDGID[9]   = -3
   292 _HERWIGID2PDGID[10]  = -4
   293 _HERWIGID2PDGID[11]  = -5
   294 _HERWIGID2PDGID[12]  = -6
   295 _HERWIGID2PDGID[13]  =  21
   296 _HERWIGID2PDGID[59]  =  22
   297 _HERWIGID2PDGID[121] =  11
   298 _HERWIGID2PDGID[122] =  12
   299 _HERWIGID2PDGID[123] =  13
   300 _HERWIGID2PDGID[124] =  14
   301 _HERWIGID2PDGID[125] =  15
   302 _HERWIGID2PDGID[126] =  16
   303 _HERWIGID2PDGID[127] = -11
   304 _HERWIGID2PDGID[128] = -12
   305 _HERWIGID2PDGID[129] = -13
   306 _HERWIGID2PDGID[130] = -14
   307 _HERWIGID2PDGID[131] = -15
   308 _HERWIGID2PDGID[132] = -16
   309 _HERWIGID2PDGID[198] =  24 # W+
   310 _HERWIGID2PDGID[199] = -24 # W-
   311 _HERWIGID2PDGID[200] =  23 # Z0
   312 _HERWIGID2PDGID[201] =  25 ## SM HIGGS
   313 _HERWIGID2PDGID[203] =  25 ## HIGGSL0 (== PDG standard in this direction)
   314 _HERWIGID2PDGID[204] =  35 ## HIGGSH0
   315 _HERWIGID2PDGID[205] =  36 ## HIGGSA0
   316 _HERWIGID2PDGID[206] =  37 ## HIGGS+
   317 _HERWIGID2PDGID[207] = -37 ## HIGGS-
   318 _HERWIGID2PDGID[401] =  1000001 ## SSDLBR
   319 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
   320 _HERWIGID2PDGID[402] =  1000002 ## SSULBR
   321 _HERWIGID2PDGID[408] = -1000002 ## SSUL
   322 _HERWIGID2PDGID[403] =  1000003 ## SSSLBR
   323 _HERWIGID2PDGID[409] = -1000003 ## SSSL
   324 _HERWIGID2PDGID[404] =  1000004 ## SSCLBR
   325 _HERWIGID2PDGID[410] = -1000004 ## SSCL
   326 _HERWIGID2PDGID[405] =  1000005 ## SSB1BR
   327 _HERWIGID2PDGID[411] = -1000005 ## SSB1
   328 _HERWIGID2PDGID[406] =  1000006 ## SST1BR
   329 _HERWIGID2PDGID[412] = -1000006 ## SST1
   330 _HERWIGID2PDGID[413] =  2000001 ## SSDR
   331 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
   332 _HERWIGID2PDGID[414] =  2000002 ## SSUR
   333 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
   334 _HERWIGID2PDGID[415] =  2000003 ## SSSR
   335 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
   336 _HERWIGID2PDGID[416] =  2000004 ## SSCR
   337 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
   338 _HERWIGID2PDGID[417] =  2000005 ## SSB2
   339 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
   340 _HERWIGID2PDGID[418] =  2000006 ## SST2
   341 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
   342 _HERWIGID2PDGID[425] =  1000011 ## SSEL-
   343 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
   344 _HERWIGID2PDGID[426] =  1000012 ## SSNUEL
   345 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
   346 _HERWIGID2PDGID[427] =  1000013 ## SSMUL-
   347 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
   348 _HERWIGID2PDGID[428] =  1000014 ## SSNUMUL
   349 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
   350 _HERWIGID2PDGID[429] =  1000015 ## SSTAU1-
   351 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
   352 _HERWIGID2PDGID[430] =  1000016 ## SSNUTL
   353 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
   354 _HERWIGID2PDGID[437] =  2000011 ## SSEL-
   355 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
   356 _HERWIGID2PDGID[438] =  2000012 ## SSNUEL
   357 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
   358 _HERWIGID2PDGID[439] =  2000013 ## SSMUL-
   359 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
   360 _HERWIGID2PDGID[440] =  2000014 ## SSNUMUL
   361 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
   362 _HERWIGID2PDGID[441] =  2000015 ## SSTAU1-
   363 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
   364 _HERWIGID2PDGID[442] =  2000016 ## SSNUTL
   365 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
   366 _HERWIGID2PDGID[449] =  1000021 ## GLUINO
   367 _HERWIGID2PDGID[450] =  1000022 ## NTLINO1
   368 _HERWIGID2PDGID[451] =  1000023 ## NTLINO2
   369 _HERWIGID2PDGID[452] =  1000025 ## NTLINO3
   370 _HERWIGID2PDGID[453] =  1000035 ## NTLINO4
   371 _HERWIGID2PDGID[454] =  1000024 ## CHGINO1+
   372 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
   373 _HERWIGID2PDGID[455] =  1000037 ## CHGINO2+
   374 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
   375 _HERWIGID2PDGID[458] =  1000039 ## GRAVTINO
   377 def herwigid2pdgid(hwid):
   378     """
   379     Convert a particle ID code in the HERWIG internal IDHW format (as used by
   380     ISAWIG) into its equivalent in the standard PDG ID code definition.
   381     """
   382     return _HERWIGID2PDGID.get(hwid, hwid)
   385 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
   386 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
   387 ## + the IDPDG array and section 4.13 of the HERWIG manual.
   388 _PDGID2HERWIGID = {}
   389 _PDGID2HERWIGID[      -1] = 7
   390 _PDGID2HERWIGID[      -2] = 8
   391 _PDGID2HERWIGID[      -3] = 9
   392 _PDGID2HERWIGID[      -4] = 10
   393 _PDGID2HERWIGID[      -5] = 11
   394 _PDGID2HERWIGID[      -6] = 12
   395 _PDGID2HERWIGID[      21] = 13
   396 _PDGID2HERWIGID[      22] = 59
   397 _PDGID2HERWIGID[      11] = 121
   398 _PDGID2HERWIGID[      12] = 122
   399 _PDGID2HERWIGID[      13] = 123
   400 _PDGID2HERWIGID[      14] = 124
   401 _PDGID2HERWIGID[      15] = 125
   402 _PDGID2HERWIGID[      16] = 126
   403 _PDGID2HERWIGID[     -11] = 127
   404 _PDGID2HERWIGID[     -12] = 128
   405 _PDGID2HERWIGID[     -13] = 129
   406 _PDGID2HERWIGID[     -14] = 130
   407 _PDGID2HERWIGID[     -15] = 131
   408 _PDGID2HERWIGID[     -16] = 132
   409 _PDGID2HERWIGID[      24] = 198 ## W+
   410 _PDGID2HERWIGID[     -24] = 199 ## W-
   411 _PDGID2HERWIGID[      23] = 200 ## Z
   412 _PDGID2HERWIGID[      25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
   413 _PDGID2HERWIGID[      26] = 203 ## HIGGSL0
   414 _PDGID2HERWIGID[      35] = 204 ## HIGGSH0
   415 _PDGID2HERWIGID[      36] = 205 ## HIGGSA0
   416 _PDGID2HERWIGID[      37] = 206 ## HIGGS+
   417 _PDGID2HERWIGID[     -37] = 207 ## HIGGS-
   418 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
   419 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
   420 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
   421 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
   422 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
   423 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
   424 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
   425 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
   426 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
   427 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
   428 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
   429 _PDGID2HERWIGID[-1000006] = 412 ## SST1
   430 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
   431 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
   432 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
   433 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
   434 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
   435 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
   436 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
   437 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
   438 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
   439 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
   440 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
   441 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
   442 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
   443 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
   444 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
   445 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
   446 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
   447 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
   448 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
   449 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
   450 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
   451 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
   452 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
   453 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
   454 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
   455 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
   456 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
   457 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
   458 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
   459 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
   460 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
   461 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
   462 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
   463 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
   464 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
   465 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
   466 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
   467 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
   468 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
   469 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
   470 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
   471 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
   472 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
   473 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
   474 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
   475 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
   477 def pdgid2herwigid(pdgid):
   478     """
   479     Convert a particle ID code in the standard PDG ID code definition into
   480     its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
   481     """
   482     return _PDGID2HERWIGID.get(pdgid, pdgid)
   485 ###############################################################################
   486 ## ISAWIG format reading/writing
   489 def readISAWIG(isastr, ignorenobr=False):
   490     """
   491     Read a spectrum definition from a string in the ISAWIG format, returning
   492     dictionaries of blocks and decays. While this is not an SLHA format, it is
   493     informally supported as a useful mechanism for converting ISAWIG spectra to
   494     SLHA.
   496     ISAWIG parsing based on the HERWIG SUSY specification format, from
   497     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   499     If the ignorenobr parameter is True, do not store decay entries with a
   500     branching ratio of zero.
   501     """
   503     blocks = {}
   504     decays = {}
   505     LINES = isastr.splitlines()
   507     def getnextvalidline():
   508         while LINES:
   509             s = LINES.pop(0).strip()
   510             ## Return None if EOF reached
   511             if len(s) == 0:
   512                 continue
   513             ## Strip comments
   514             if "#" in s:
   515                 s = s[:s.index("#")].strip()
   516             ## Return if non-empty
   517             if len(s) > 0:
   518                 return s
   520     def getnextvalidlineitems():
   521         return map(_autotype, getnextvalidline().split())
   523     ## Populate MASS block and create decaying particle objects
   524     masses = Block("MASS")
   525     numentries = int(getnextvalidline())
   526     for i in xrange(numentries):
   527         hwid, mass, lifetime = getnextvalidlineitems()
   528         width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
   529         pdgid = herwigid2pdgid(hwid)
   530         masses.add_entry((pdgid, mass))
   531         decays[pdgid] = Particle(pdgid, width, mass)
   532         #print pdgid, mass, width
   533     blocks["MASS"] = masses
   535     ## Populate decays
   536     for n in xrange(numentries):
   537         numdecays = int(getnextvalidline())
   538         for d in xrange(numdecays):
   539             #print n, numentries-1, d, numdecays-1
   540             decayitems = getnextvalidlineitems()
   541             hwid = decayitems[0]
   542             pdgid = herwigid2pdgid(hwid)
   543             br = decayitems[1]
   544             nme = decayitems[2]
   545             daughter_hwids = decayitems[3:]
   546             daughter_pdgids = []
   547             for hw in daughter_hwids:
   548                 if hw != 0:
   549                     daughter_pdgids.append(herwigid2pdgid(hw))
   550             if not decays.has_key(pdgid):
   551                 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
   552                 decays[pdgid] = Particle(pdgid)
   553             decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
   556     ## Now the SUSY parameters
   557     TANB, ALPHAH = getnextvalidlineitems()
   558     blocks["MINPAR"] = Block("MINPAR")
   559     blocks["MINPAR"].add_entry((3, TANB))
   560     blocks["ALPHA"] = Block("ALPHA")
   561     blocks["ALPHA"].entries = ALPHAH
   562     #
   563     ## Neutralino mixing matrix
   564     blocks["NMIX"] = Block("NMIX")
   565     for i in xrange(1, 5):
   566         nmix_i = getnextvalidlineitems()
   567         for j, v in enumerate(nmix_i):
   568             blocks["NMIX"].add_entry((i, j+1, v))
   569     #
   570     ## Chargino mixing matrices V and U
   571     blocks["VMIX"] = Block("VMIX")
   572     vmix = getnextvalidlineitems()
   573     blocks["VMIX"].add_entry((1, 1, vmix[0]))
   574     blocks["VMIX"].add_entry((1, 2, vmix[1]))
   575     blocks["VMIX"].add_entry((2, 1, vmix[2]))
   576     blocks["VMIX"].add_entry((2, 2, vmix[3]))
   577     blocks["UMIX"] = Block("UMIX")
   578     umix = getnextvalidlineitems()
   579     blocks["UMIX"].add_entry((1, 1, umix[0]))
   580     blocks["UMIX"].add_entry((1, 2, umix[1]))
   581     blocks["UMIX"].add_entry((2, 1, umix[2]))
   582     blocks["UMIX"].add_entry((2, 2, umix[3]))
   583     #
   584     THETAT, THETAB, THETAL = getnextvalidlineitems()
   585     import math
   586     blocks["STOPMIX"] = Block("STOPMIX")
   587     blocks["STOPMIX"].add_entry((1, 1,  math.cos(THETAT)))
   588     blocks["STOPMIX"].add_entry((1, 2, -math.sin(THETAT)))
   589     blocks["STOPMIX"].add_entry((2, 1,  math.sin(THETAT)))
   590     blocks["STOPMIX"].add_entry((2, 2,  math.cos(THETAT)))
   591     blocks["SBOTMIX"] = Block("SBOTMIX")
   592     blocks["SBOTMIX"].add_entry((1, 1,  math.cos(THETAB)))
   593     blocks["SBOTMIX"].add_entry((1, 2, -math.sin(THETAB)))
   594     blocks["SBOTMIX"].add_entry((2, 1,  math.sin(THETAB)))
   595     blocks["SBOTMIX"].add_entry((2, 2,  math.cos(THETAB)))
   596     blocks["STAUMIX"] = Block("STAUMIX")
   597     blocks["STAUMIX"].add_entry((1, 1,  math.cos(THETAL)))
   598     blocks["STAUMIX"].add_entry((1, 2, -math.sin(THETAL)))
   599     blocks["STAUMIX"].add_entry((2, 1,  math.sin(THETAL)))
   600     blocks["STAUMIX"].add_entry((2, 2,  math.cos(THETAL)))
   601     #
   602     ATSS, ABSS, ALSS = getnextvalidlineitems()
   603     blocks["AU"] = Block("AU")
   604     blocks["AU"].add_entry((3, 3, ATSS))
   605     blocks["AD"] = Block("AD")
   606     blocks["AD"].add_entry((3, 3, ABSS))
   607     blocks["AE"] = Block("AE")
   608     blocks["AE"].add_entry((3, 3, ALSS))
   609     #
   610     MUSS = getnextvalidlineitems()[0]
   611     blocks["MINPAR"].add_entry((4, MUSS))
   612     #
   614     # TODO: Parse RPV boolean and couplings into SLHA2 blocks
   616     return blocks, decays
   619 def writeISAJET(blocks, decays, outname, ignorenobr=False, precision=8):
   620     """
   621     Return a SUSY spectrum definition in the format required for input by ISAJET,
   622     as a string, from the supplied blocks and decays dicts.
   624     The outname parameter specifies the desired output filename from ISAJET: this
   625     will appear in the first line of the return value.
   627     If the ignorenobr parameter is True, do not write decay entries with a
   628     branching ratio of zero.
   629     """
   630     fmte = "%." + str(precision) + "e"
   632     masses = blocks["MASS"].entries
   634     ## Init output string
   635     out = ""
   637     ## First line is the output name
   638     out += "'%s'" % outname + "\n"
   640     ## Next the top mass
   641     out += fmte % masses[6] + "\n"
   643     ## Next the top mass
   644     out += fmte % masses[6] + "\n"
   646     ## mSUGRA parameters (one line)
   647     # e.g. 1273.78,713.286,804.721,4.82337
   649     ## Masses and trilinear couplings (3 lines)
   650     # e.g. 1163.14,1114.15,1118.99,374.664,209.593
   651     # e.g. 1069.54,1112.7,919.908,374.556,209.381,-972.817,-326.745,-406.494
   652     # e.g. 1163.14,1114.15,1118.99,374.712,210.328
   654     ## RPV couplings (?? lines)
   655     # e.g. 232.615,445.477
   657     ## Etc ???!!!
   658     # e.g. /
   659     # e.g. n
   660     # e.g. y
   661     # e.g. y
   662     # e.g. 0.047441 3.80202e-23 0 0 0 2.17356e-22 0 0 5.23773e-09
   663     # e.g. y
   664     # 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
   665     # e.g. n
   666     # e.g. 'susy_RPV_stau_BC1scan_m560_tanb05.txt'
   668     return out
   671 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
   672     """
   673     Return a SUSY spectrum definition in the format produced by ISAWIG for inut to HERWIG
   674     as a string, from the supplied SLHA blocks and decays dicts.
   676     ISAWIG parsing based on the HERWIG SUSY specification format, from
   677     http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
   679     If the ignorenobr parameter is True, do not write decay entries with a
   680     branching ratio of zero.
   681     """
   682     fmte = "%." + str(precision) + "e"
   684     masses = blocks["MASS"].entries
   686     ## Init output string
   687     out = ""
   689     ## First write out masses section:
   690     ##   Number of SUSY + top particles
   691     ##   IDHW, RMASS(IDHW), RLTIM(IDHW)
   692     ##   repeated for each particle
   693     ## IDHW is the HERWIG identity code.
   694     ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
   695     massout = ""
   696     for pid in masses.keys():
   697         lifetime = -1
   698         try:
   699             width = decays[pid].totalwidth
   700             if width and width > 0:
   701                 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
   702         except:
   703             pass
   704         massout += ("%d " + fmte + " " + fmte + "\n") % (pdgid2herwigid(pid), masses[pid], lifetime)
   705     out += "%d\n" % massout.count("\n")
   706     out += massout
   708     assert(len(masses) == len(decays))
   710     ## Next each particles decay modes together with their branching ratios and matrix element codes
   711     ##   Number of decay modes for a given particle (IDK)
   712     ##     IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
   713     ##     repeated for each mode.
   714     ##   Repeated for each particle.
   715     ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
   716     ## the decay mode. NME is a code for the matrix element to be used, either from the
   717     ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
   718     for i, pid in enumerate(decays.keys()):
   719         # if not decays.has_key(pid):
   720         #     continue
   721         hwid = pdgid2herwigid(pid)
   722         decayout = ""
   723         #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
   724         for i_d, d in enumerate(decays[pid].decays):
   725             ## Skip decay if it has no branching ratio
   726             if ignorenobr and d.br == 0:
   727                 continue
   729             ## Identify decay matrix element to use
   730             ## From std HW docs, or from this pair:
   731             ## Two new matrix element codes have been added for these new decays:
   732             ##    NME =	200 	3 body top quark via charged Higgs
   733             ##    	300 	3 body R-parity violating gaugino and gluino decays
   734             nme = 0
   735             # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
   736             if abs(pid) in (6, 12):
   737                 nme = 100
   738             ## Extra SUSY MEs
   739             if len(d.ids) == 3:
   740                 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
   741                 pass
   742             decayout += "%d " + fmte + " %d " % (hwid, d.br, nme)
   744             def is_quark(pid):
   745                 return (abs(pid) in range(1, 7))
   747             def is_lepton(pid):
   748                 return (abs(pid) in range(11, 17))
   750             def is_squark(pid):
   751                 if abs(pid) in range(1000001, 1000007):
   752                     return True
   753                 if abs(pid) in range(2000001, 2000007):
   754                     return True
   755                 return False
   757             def is_slepton(pid):
   758                 if abs(pid) in range(1000011, 1000017):
   759                     return True
   760                 if abs(pid) in range(2000011, 2000016, 2):
   761                     return True
   762                 return False
   764             def is_gaugino(pid):
   765                 if abs(pid) in range(1000022, 1000026):
   766                     return True
   767                 if abs(pid) in (1000035, 1000037):
   768                     return True
   769                 return False
   771             def is_susy(pid):
   772                 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
   774             absids = map(abs, d.ids)
   776             ## Order decay products as required by HERWIG
   777             ## Top
   778             if abs(pid) == 6:
   779                 def cmp_bottomlast(a, b):
   780                     """Comparison function which always puts b/bbar last"""
   781                     if abs(a) == 5:
   782                         return True
   783                     if abs(b) == 5:
   784                         return False
   785                     return cmp(a, b)
   786                 if len(absids) == 2:
   787                     ## 2 body mode, to Higgs: Higgs; Bottom
   788                     if (25 in absids or 26 in absids) and 5 in absids:
   789                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   790                 elif len(absids) == 3:
   791                     ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
   792                     if 37 in absids or 23 in absids:
   793                         d.ids = sorted(d.ids, key=cmp_bottomlast)
   794             ## Gluino
   795             elif abs(pid) == 1000021:
   796                 if len(absids) == 2:
   797                     ## 2 body mode
   798                     ## without gluon: any order
   799                     ## with gluon: gluon; colour neutral
   800                     if 21 in absids:
   801                         def cmp_gluonfirst(a, b):
   802                             """Comparison function which always puts gluon first"""
   803                             if a == 21:
   804                                 return False
   805                             if b == 21:
   806                                 return True
   807                             return cmp(a, b)
   808                         d.ids = sorted(d.ids, key=cmp_gluonfirst)
   809                 elif len(absids) == 3:
   810                     ## 3-body modes, R-parity conserved: colour neutral; q or qbar
   811                     def cmp_quarkslast(a, b):
   812                         """Comparison function which always puts quarks last"""
   813                         if is_quark(a):
   814                             return True
   815                         if is_quark(b):
   816                             return False
   817                         return cmp(a, b)
   818                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   819             ## Squark/Slepton
   820             elif is_squark(pid) or is_slepton(pid):
   821                 def cmp_susy_quark_lepton(a, b):
   822                     if is_susy(a):
   823                         return False
   824                     if is_susy(b):
   825                         return True
   826                     if is_quark(a):
   827                         return False
   828                     if is_quark(b):
   829                         return True
   830                     return cmp(a, b)
   831                 ##   2 body modes: Gaugino/Gluino with Quark/Lepton     Gaugino      quark
   832                 ##                                                      Gluino       lepton
   833                 ##   3 body modes: Weak                                 sparticle    particles from W decay
   834                 ## Squark
   835                 ##   2 body modes: Lepton Number Violated               quark     lepton
   836                 ##                 Baryon number violated               quark     quark
   837                 ## Slepton
   838                 ##   2 body modes: Lepton Number Violated               q or qbar
   839                 d.ids = sorted(d.ids, key=cmp_bottomlast)
   840             ## Higgs
   841             elif pid in (25, 26):
   842                 # TODO: Includes SUSY Higgses?
   843                 ## Higgs
   844                 ##   2 body modes: (s)quark-(s)qbar                     (s)q or (s)qbar
   845                 ##                 (s)lepton-(s)lepton                  (s)l or (s)lbar
   846                 ##   3 body modes:                                      colour neutral       q or qbar
   847                 if len(absids) == 3:
   848                     def cmp_quarkslast(a, b):
   849                         """Comparison function which always puts quarks last"""
   850                         if is_quark(a):
   851                             return True
   852                         if is_quark(b):
   853                             return False
   854                         return cmp(a, b)
   855                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   856             elif is_gaugino(pid):
   857                 # TODO: Is there actually anything to do here?
   858                 ## Gaugino
   859                 ##   2 body modes: Squark-quark                         q or sq
   860                 ##                 Slepton-lepton                       l or sl
   861                 ##
   862                 ##   3 body modes: R-parity conserved                   colour neutral       q or qbar
   863                 ##                                                                           l or lbar
   864                 if len(absids) == 3:
   865                     def cmp_quarkslast(a, b):
   866                         """Comparison function which always puts quarks last"""
   867                         if is_quark(a):
   868                             return True
   869                         if is_quark(b):
   870                             return False
   871                         return cmp(a, b)
   872                     d.ids = sorted(d.ids, key=cmp_quarkslast)
   874             # TODO: Gaugino/Gluino
   875             ##   3 body modes:  R-parity violating:   Particles in the same order as the R-parity violating superpotential
   877             ## Pad out IDs list with zeros
   878             ids = [0,0,0,0,0]
   879             for i, pid in enumerate(d.ids):
   880                 ids[i] = pid
   881             ids = map(str, ids)
   882             decayout += " ".join(ids) + "\n"
   883         decayout = "%d\n" % decayout.count("\n") + decayout
   884         out += decayout
   886     ## Now the SUSY parameters
   887     ## TANB, ALPHAH:
   888     out += (fmte + " " + fmte + "\n") % (blocks["MINPAR"].entries[3], blocks["ALPHA"].entries)
   889     ## Neutralino mixing matrix
   890     nmix = blocks["NMIX"].entries
   891     for i in xrange(1, 5):
   892         out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (nmix[i][1], nmix[i][2], nmix[i][3], nmix[i][4])
   893     ## Chargino mixing matrices V and U
   894     vmix = blocks["VMIX"].entries
   895     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (vmix[1][1], vmix[1][2], vmix[2][1], vmix[2][2])
   896     umix = blocks["UMIX"].entries
   897     out += (fmte + " " + fmte + " " + fmte + " " + fmte + "\n") % (umix[1][1], umix[1][2], umix[2][1], umix[2][2])
   898     ## THETAT,THETAB,THETAL
   899     import math
   900     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (math.acos(blocks["STOPMIX"].entries[1][1]),
   901                                                             math.acos(blocks["SBOTMIX"].entries[1][1]),
   902                                                             math.acos(blocks["STAUMIX"].entries[1][1]))
   903     ## ATSS,ABSS,ALSS
   904     out += (fmte + " " + fmte + " " + fmte + " " + "\n") % (blocks["AU"].entries[3][3],
   905                                                             blocks["AD"].entries[3][3],
   906                                                             blocks["AE"].entries[3][3])
   907     ## MUSS == sign(mu)
   908     out += "%f\n" % blocks["MINPAR"].entries[4]
   910     ## RPV SUSY
   911     isRPV = False
   912     out += "%d\n" % isRPV
   913     # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
   914     # TODO: Get the index orderings right
   915     # if isRPV: ...
   917     return out
   920 ###############################################################################
   921 ## File-level functions
   924 def readSLHAFile(spcfilename, **kwargs):
   925     """
   926     Read an SLHA file, returning dictionaries of blocks and decays.
   928     Other keyword parameters are passed to readSLHA.
   929     """
   930     f = open(spcfilename, "r")
   931     rtn = readSLHA(f.read(), kwargs)
   932     f.close()
   933     return rtn
   936 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
   937     """
   938     Write an SLHA file from the supplied blocks and decays dicts.
   940     Other keyword parameters are passed to writeSLHA.
   941     """
   942     f = open(spcfilename, "w")
   943     f.write(writeSLHA(blocks, decays, kwargs))
   944     f.close()
   947 def readISAWIGFile(isafilename, **kwargs):
   948     """
   949     Read a spectrum definition from a file in the ISAWIG format, returning
   950     dictionaries of blocks and decays. While this is not an SLHA format, it is
   951     informally supported as a useful mechanism for converting ISAWIG spectra to
   952     SLHA.
   954     Other keyword parameters are passed to readSLHA.
   955     """
   956     f = open(isafilename, "r")
   957     rtn = readISAWIG(f.read(), kwargs)
   958     f.close()
   959     return rtn
   962 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
   963     """
   964     Write an ISAWIG file from the supplied blocks and decays dicts.
   966     Other keyword parameters are passed to writeISAWIG.
   967     """
   968     f = open(isafilename, "w")
   969     f.write(writeISAWIG(blocks, decays, kwargs))
   970     f.close()
   973 def writeISAJETFile(isafilename, blocks, decays, **kwargs):
   974     """
   975     Write an ISAJET file from the supplied blocks and decays dicts (see writeISAJET).
   977     Other keyword parameters are passed to writeISAJET.
   978     """
   979     f = open(isafilename, "w")
   980     f.write(writeISAWIG(blocks, decays, kwargs))
   981     f.close()
   985 ###############################################################################
   986 ## Main function for module testing
   989 if __name__ == "__main__":
   990     import sys
   991     for a in sys.argv[1:]:
   992         if a.endswith(".isa"):
   993             blocks, decays = readISAWIGFile(a)
   994         else:
   995             blocks, decays = readSLHAFile(a)
   997         for bname, b in sorted(blocks.iteritems()):
   998             print b
   999             print
  1001         print blocks.keys()
  1003         print blocks["MASS"].entries[25]
  1004         print
  1006         for p in sorted(decays.values()):
  1007             print p
  1008             print
  1010         print writeSLHA(blocks, decays, ignorenobr=True)

mercurial