pyslha.py

Sun, 28 Apr 2013 21:56:51 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 21:56:51 +0200
changeset 202
d3d069f4549a
parent 201
3ac555efdbba
child 203
ce90a0dace07
permissions
-rw-r--r--

Removing inappropriate uses of add_entry

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

mercurial