pyslha.py

Sun, 28 Apr 2013 21:13:25 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 21:13:25 +0200
changeset 199
80d4b063103a
parent 198
6a10f5f8a9b0
child 200
651b0fac2163
permissions
-rw-r--r--

Various fixes to ISAWIG output (and conversion to use new tuple-indexing on mixing matrices and single-value accessing on ALPHA).

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

mercurial