pyslha.py

Mon, 15 Jul 2013 16:50:38 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 15 Jul 2013 16:50:38 +0200
changeset 244
269bc4611cdd
parent 242
41ec8fb352b3
child 245
ca52fe2c0d71
permissions
-rw-r--r--

Add credit to ChangeLog

andy@1 1 #! /usr/bin/env python
andy@1 2
andy@8 3 """\
andy@208 4 A simple but flexible handler of the SUSY Les Houches Accord (SLHA) data format.
andy@25 5
andy@52 6 pyslha is a parser/writer module for particle physics SUSY Les Houches Accord
andy@52 7 (SLHA) supersymmetric spectrum/decay files, and a collection of scripts which
andy@52 8 use the interface, e.g. for conversion to and from the legacy ISAWIG format, or
andy@52 9 to plot the mass spectrum and decay chains.
andy@52 10
andy@192 11 The current release supports SLHA version 1, and as far as I'm aware is also
andy@145 12 fully compatible with SLHA2: the block structures are read and accessed
andy@192 13 generically. If you have any problems, please provide an example input file and
andy@208 14 I'll happily investigate.
andy@52 15
andy@131 16 The plotting script provides output in PDF, EPS and PNG via LaTeX and the TikZ
andy@131 17 graphics package, and as LaTeX/TikZ source for direct embedding into documents or
andy@131 18 user-tweaking of the generated output.
andy@64 19
andy@208 20 Users of version 1.x should note that the interface has changed a little in
andy@208 21 version 2.0.0 and onward, in particular in the interface of the Block objects,
andy@208 22 which are now more dict-like: entries can be added and accessed via the usual
andy@208 23 square-brackets indexing operators, including for multiple indices as is common
andy@208 24 for mixing matrices e.g. NMIX[1,2] as opposed to the old NMIX.entries[1][2]
andy@208 25 way. This does break backward compatibility but is a big improvement both for
andy@208 26 internal code sanity and usability. The Block interface also now supplies
andy@208 27 dict-like has_key(), keys(), and items() methods, as well as more specialist
andy@208 28 value(), set_value() and is_single_valued() methods for improved access to ALPHA
andy@208 29 and any other unindexed blocks.
andy@208 30
andy@236 31 If you use PySLHA, for either model data handling or spectrum visualisation,
andy@236 32 please cite the paper: http://arxiv.org/abs/1305.4194
andy@208 33
andy@52 34 TODOs:
andy@192 35
andy@239 36 For 2.1.4:
andy@235 37 * In set_value, if first item is non-int, treat as None-indexed
andy@235 38 * Refine value string heuristic for strings with ints in them?
andy@231 39
andy@229 40 For 2.2.0:
andy@225 41 * Add Sphinx docs.
andy@212 42 * Preserve comments from read -> write (needs full-line/inline comment
andy@223 43 separation). Can use separate comment dicts in Block and Decay, and
andy@212 44 attach a multiline .comment attr to the returned/written dicts.
andy@198 45
andy@239 46 For 3.0.0:
andy@239 47 * Add a new Doc type which encapsulates the blocks, decays, xsections,
andy@239 48 etc. data sections and will act as a more robust single return value from
andy@239 49 the read commands / argument to the write commands. Can also handle
andy@239 50 document-level header comments this way and improve the __repr__ resulting
andy@239 51 from an interactive read() call without a target variable (i.e. hide the
andy@239 52 uninteresting OrderedDict boilerplate).
andy@239 53
andy@208 54 Later/maybe:
andy@208 55 * Identify HERWIG decay matrix element to use in ISAWIG
andy@208 56 * Handle RPV SUSY in ISAWIG
andy@8 57 """
andy@8 58
andy@223 59 __author__ = "Andy Buckley <andy.buckley@cern.ch>"
andy@242 60 __version__ = "2.1.3"
andy@8 61
andy@1 62
andy@223 63
andy@223 64 ###############################################################################
andy@223 65 ## Private utility functions
andy@223 66
andy@189 67 def _mkdict():
andy@189 68 """Try to return an empty ordered dict, but fall back to normal dict if necessary"""
andy@189 69 try:
andy@189 70 from collections import OrderedDict
andy@189 71 return OrderedDict()
andy@189 72 except:
andy@189 73 try:
andy@189 74 from ordereddict import OrderedDict
andy@189 75 return OrderedDict()
andy@189 76 except:
andy@189 77 return dict()
andy@189 78
andy@4 79 def _autotype(var):
andy@30 80 """Automatically convert strings to numerical types if possible."""
andy@30 81 if type(var) is not str:
andy@4 82 return var
andy@36 83 if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
andy@4 84 return int(var)
andy@4 85 try:
andy@4 86 f = float(var)
andy@4 87 return f
andy@4 88 except ValueError:
andy@4 89 return var
andy@4 90
andy@147 91 def _autostr(var, precision=8):
andy@192 92 """Automatically format numerical types as the right sort of string."""
andy@30 93 if type(var) is float:
andy@213 94 return ("% ." + str(precision) + "e") % var
andy@231 95 elif not hasattr(var, "__iter__"):
andy@231 96 return str(var)
andy@231 97 else:
andy@231 98 return ",".join(_autostr(subval) for subval in var)
andy@30 99
andy@221 100 def _autotuple(a):
andy@221 101 """Automatically convert the supplied iterable to a scalar or tuple as appropriate."""
andy@221 102 if len(a) == 1:
andy@221 103 return a[0]
andy@221 104 return tuple(a)
andy@221 105
andy@30 106
andy@223 107
andy@223 108 ###############################################################################
andy@223 109 ## Exceptions
andy@223 110
andy@193 111 class AccessError(Exception):
andy@193 112 "Exception object to be raised when a SLHA block is accessed in an invalid way"
andy@193 113 def __init__(self, errmsg):
andy@193 114 self.msg = errmsg
andy@193 115 def __str__(self):
andy@193 116 return self.msg
andy@193 117
andy@166 118 class ParseError(Exception):
andy@166 119 "Exception object to be raised when a spectrum file/string is malformed"
andy@166 120 def __init__(self, errmsg):
andy@166 121 self.msg = errmsg
andy@168 122 def __str__(self):
andy@166 123 return self.msg
andy@166 124
andy@166 125
andy@4 126
andy@223 127 ###############################################################################
andy@223 128 ## The data block, decay and particle classes
andy@223 129
andy@12 130 class Block(object):
andy@8 131 """
andy@203 132 Object representation of any BLOCK elements read from an SLHA file.
andy@203 133
andy@203 134 Blocks have a name, may have an associated Q value, and contain a collection
andy@203 135 of data entries, each indexed by one or more keys. Entries in the dictionary
andy@203 136 are stored as numeric types (int or float) when a cast from the string in
andy@203 137 the file has been possible.
andy@203 138
andy@203 139 Block is closely related to a Python dict (and, in fact, is implemented via
andy@203 140 an OrderedDict when possible). The preferred methods of entry access use the
andy@203 141 dict-like [] operator for getting and setting, and the keys() and items()
andy@203 142 methods for iteration. Purely iterating over the object behaves like keys(),
andy@203 143 as for an ordinary dict.
andy@203 144
andy@223 145 Multiple (integer) indices are possible, especially for entries in mixing
andy@223 146 matrix blocks. These are now implemented in the natural way, e.g. for access
andy@223 147 to the (1,2) element of a mixing matrix block, use bmix[1,2] = 0.123 and
andy@223 148 print bmix[1,2]. The value() and set_value() functions behave
andy@223 149 similarly. Multi-element values are also permitted.
andy@203 150
andy@203 151 It is possible, although not usual, to store unindexed values in a
andy@203 152 block. This is only supported when that entry is the only one in the block,
andy@203 153 and it is stored in the normal fashion but with None as the lookup key. The
andy@203 154 value() method may be used without a key argument to retrieve this value, if
andy@203 155 the is_single_valued() method returns True, and similarly the set_value()
andy@203 156 method may be used to set it if only one argument is supplied and the object
andy@203 157 is compatible.
andy@8 158 """
andy@204 159 def __init__(self, name, q=None, entries=None):
andy@1 160 self.name = name
andy@189 161 self.entries = _mkdict()
andy@204 162 if entries is not None:
andy@204 163 self.entries.update(entries)
andy@5 164 self.q = _autotype(q)
andy@1 165
andy@219 166 # TODO: Rename? To what?
andy@219 167 def add_entry(self, args):
andy@219 168 """Add an entry to the block from an iterable (i.e. list or tuple) of
andy@219 169 strings, or from a whitespace-separated string.
andy@203 170
andy@219 171 This method is just for convenience: it splits the single string
andy@219 172 argument if necessary and converts the list of strings into numeric
andy@223 173 types when possible. For the treatment of the resulting iterable see the
andy@223 174 set_value method.
andy@196 175 """
andy@203 176 ## If the argument is a single string, split it and proceed
andy@219 177 if type(args) is str:
andy@219 178 args = args.split()
andy@203 179 ## Check that the arg is an iterable
andy@219 180 if not hasattr(args, "__iter__"):
andy@196 181 raise AccessError("Block entries must be iterable")
andy@203 182 ## Auto-convert the types in the list
andy@219 183 args = map(_autotype, args)
andy@232 184 ## Re-join consecutive strings into single entries
andy@232 185 i = 0
andy@232 186 while i < len(args)-1:
andy@232 187 if type(args[i]) is str and type(args[i+1]) is str:
andy@232 188 args[i] += " " + args[i+1]
andy@232 189 del args[i+1]
andy@232 190 continue
andy@232 191 i += 1
andy@203 192 ## Add the entry to the map, with appropriate indices
andy@221 193 self.set_value(*args)
andy@192 194
andy@196 195 def is_single_valued(self):
andy@196 196 """Return true if there is only one entry, and it has no index: the
andy@196 197 'value()' attribute may be used in that case without an argument."""
andy@199 198 return len(self.entries) == 1 and self.entries.keys()[0] is None
andy@196 199
andy@196 200 def value(self, key=None):
andy@196 201 """Get a value from the block with the supplied key.
andy@196 202
andy@196 203 If no key is given, then the block must contain only one non-indexed
andy@196 204 value otherwise an AccessError exception will be raised.\
andy@196 205 """
andy@196 206 if key == None and not self.is_single_valued():
andy@196 207 raise AccessError("Tried to access unique value of multi-value block")
andy@196 208 return self.entries[key]
andy@196 209
andy@203 210 def set_value(self, *args):
andy@219 211 """Set a value in the block via supplied key/value arguments.
andy@203 212
andy@218 213 Indexing is determined automatically: any leading integers will be
andy@218 214 treated as a multi-dimensional index, with the remaining entries being a
andy@219 215 (potentially multi-dimensional) value. If all N args are ints, then the
andy@219 216 first N-1 are treated as the index and the Nth as the value.
andy@218 217
andy@219 218 If there is only one arg it will be treated as the value of a
andy@219 219 single-valued block. In this case the block must already contain at most
andy@219 220 one non-indexed value otherwise an AccessError exception will be
andy@219 221 raised.\
andy@203 222 """
andy@203 223 if len(args) == 0:
andy@203 224 raise AccessError("set_value() called without arguments")
andy@203 225 elif len(args) == 1:
andy@203 226 if len(self.entries) > 0 and not self.is_single_valued():
andy@203 227 raise AccessError("Tried to set a unique value on a multi-value block")
andy@203 228 self.entries[None] = args[0]
andy@203 229 else:
andy@218 230 ## Find the first non-integer -- all previous items are indices
andy@222 231 i_first_nonint = -1
andy@221 232 for i, x in enumerate(args):
andy@218 233 if type(x) is not int:
andy@218 234 i_first_nonint = i
andy@218 235 break
andy@218 236 if i_first_nonint == 0:
andy@221 237 raise AccessError("Attempted to set a block entry with a non-integer(s) index")
andy@218 238 else:
andy@222 239 self.entries[_autotuple(args[:i_first_nonint])] = _autotuple(args[i_first_nonint:])
andy@203 240
andy@205 241 def has_key(self, key):
andy@205 242 """Does the block have the given key?"""
andy@205 243 return self.entries.has_key(key)
andy@205 244
andy@200 245 def keys(self):
andy@200 246 """Access the block item keys."""
andy@201 247 return self.entries.keys()
andy@200 248
andy@231 249 def values(self):
andy@231 250 """Access the block item values."""
andy@231 251 return self.entries.values()
andy@231 252
andy@195 253 def items(self, key=None):
andy@196 254 """Access the block items as (key,value) tuples.
andy@196 255
andy@196 256 Note: The Python 3 dict attribute 'items()' is used rather than the
andy@196 257 'old' Python 2 'iteritems()' name for forward-looking compatibility.\
andy@196 258 """
andy@193 259 return self.entries.iteritems()
andy@193 260
andy@193 261 def __len__(self):
andy@193 262 return len(self.entries)
andy@193 263
andy@196 264 def __iter(self):
andy@196 265 return self.entries.__iter__()
andy@196 266
andy@192 267 def __getitem__(self, key):
andy@192 268 return self.entries[key]
andy@1 269
andy@201 270 def __setitem__(self, key, value):
andy@222 271 if key is not None and type(key) is not int and not all(type(x) is int for x in key):
andy@218 272 raise AccessError("Attempted to set a block entry with a non-integer(s) index")
andy@201 273 self.entries[key] = value
andy@201 274
andy@1 275 def __cmp__(self, other):
andy@196 276 return cmp(self.name, other.name) and cmp(self.entries, other.entries)
andy@1 277
andy@193 278 def __repr__(self):
andy@1 279 s = self.name
andy@5 280 if self.q is not None:
andy@5 281 s += " (Q=%s)" % self.q
andy@231 282 s += " { " + "; ".join(_autostr(k) + " : " + _autostr(v) for k, v in self.items()) + " }"
andy@1 283 return s
andy@1 284
andy@1 285
andy@12 286 class Decay(object):
andy@8 287 """
andy@8 288 Object representing a decay entry on a particle decribed by the SLHA file.
andy@8 289 'Decay' objects are not a direct representation of a DECAY block in an SLHA
andy@8 290 file... that role, somewhat confusingly, is taken by the Particle class.
andy@8 291
andy@8 292 Decay objects have three properties: a branching ratio, br, an nda number
andy@12 293 (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
andy@12 294 decay occurs. The PDG ID of the particle whose decay this represents may
andy@12 295 also be stored, but this is normally known via the Particle in which the
andy@12 296 decay is stored.
andy@8 297 """
andy@8 298 def __init__(self, br, nda, ids, parentid=None):
andy@8 299 self.parentid = parentid
andy@6 300 self.br = br
andy@6 301 self.nda = nda
andy@6 302 self.ids = ids
andy@29 303 assert(self.nda == len(self.ids))
andy@6 304
andy@6 305 def __cmp__(self, other):
andy@31 306 return cmp(other.br, self.br)
andy@6 307
andy@215 308 def __repr__(self):
andy@233 309 return "% .8e %s" % (self.br, self.ids)
andy@6 310
andy@6 311
andy@12 312 class Particle(object):
andy@8 313 """
andy@8 314 Representation of a single, specific particle, decay block from an SLHA
andy@8 315 file. These objects are not themselves called 'Decay', since that concept
andy@8 316 applies more naturally to the various decays found inside this
andy@8 317 object. Particle classes store the PDG ID (pid) of the particle being
andy@8 318 represented, and optionally the mass (mass) and total decay width
andy@8 319 (totalwidth) of that particle in the SLHA scenario. Masses may also be found
andy@8 320 via the MASS block, from which the Particle.mass property is filled, if at
andy@8 321 all. They also store a list of Decay objects (decays) which are probably the
andy@8 322 item of most interest.
andy@8 323 """
andy@6 324 def __init__(self, pid, totalwidth=None, mass=None):
andy@6 325 self.pid = pid
andy@6 326 self.totalwidth = totalwidth
andy@6 327 self.mass = mass
andy@6 328 self.decays = []
andy@6 329
andy@6 330 def add_decay(self, br, nda, ids):
andy@6 331 self.decays.append(Decay(br, nda, ids))
andy@6 332 self.decays.sort()
andy@6 333
andy@6 334 def __cmp__(self, other):
andy@6 335 if abs(self.pid) == abs(other.pid):
andy@31 336 return cmp(self.pid, other.pid)
andy@31 337 return cmp(abs(self.pid), abs(other.pid))
andy@6 338
andy@215 339 def __repr__(self):
andy@6 340 s = str(self.pid)
andy@7 341 if self.mass is not None:
andy@147 342 s += " : mass = %.8e GeV" % self.mass
andy@6 343 if self.totalwidth is not None:
andy@147 344 s += " : total width = %.8e GeV" % self.totalwidth
andy@6 345 for d in self.decays:
andy@12 346 if d.br > 0.0:
andy@12 347 s += "\n %s" % d
andy@6 348 return s
andy@1 349
andy@21 350
andy@223 351
andy@223 352 ###############################################################################
andy@223 353 ## SLHA parsing and writing functions
andy@223 354
andy@241 355 def readSLHA(spcstr, ignorenobr=False, ignorenomass=False):
andy@21 356 """
andy@31 357 Read an SLHA definition from a string, returning dictionaries of blocks and
andy@31 358 decays.
andy@31 359
andy@31 360 If the ignorenobr parameter is True, do not store decay entries with a
andy@31 361 branching ratio of zero.
andy@241 362
andy@241 363 If the ignorenomass parameter is True, parse file even if mass block is
andy@241 364 absent in the file (default is to raise a ParseError).
andy@21 365 """
andy@189 366 blocks = _mkdict()
andy@189 367 decays = _mkdict()
andy@21 368 #
andy@34 369 import re
andy@21 370 currentblock = None
andy@21 371 currentdecay = None
andy@21 372 for line in spcstr.splitlines():
andy@21 373 ## Handle (ignore) comment lines
andy@192 374 # TODO: Store block/entry comments
andy@21 375 if line.startswith("#"):
andy@21 376 continue
andy@21 377 if "#" in line:
andy@21 378 line = line[:line.index("#")]
andy@238 379 ## Ignore empty lines (after comment removal and whitespace trimming)
andy@238 380 if not line.strip():
andy@238 381 continue
andy@21 382
andy@238 383 ## Section header lines start with a non-whitespace character, data lines have a whitespace indent
andy@238 384 # TODO: Are tabs also allowed for indents? Check the SLHA standard.
andy@238 385 if not line.startswith(" "):
andy@238 386 # TODO: Should we now strip the line to remove any trailing whitespace?
andy@238 387
andy@238 388 ## Handle BLOCK start lines
andy@238 389 if line.upper().startswith("BLOCK"):
andy@238 390 #print line
andy@238 391 match = re.match(r"BLOCK\s+(\w+)(\s+Q\s*=\s*.+)?", line.upper())
andy@238 392 if not match:
andy@238 393 continue
andy@238 394 blockname = match.group(1)
andy@238 395 qstr = match.group(2)
andy@238 396 if qstr is not None:
andy@238 397 qstr = qstr[qstr.find("=")+1:].strip()
andy@238 398 currentblock = blockname
andy@238 399 currentdecay = None
andy@238 400 blocks[blockname] = Block(blockname, q=qstr)
andy@238 401 ## Handle DECAY start lines
andy@238 402 elif line.upper().startswith("DECAY"):
andy@238 403 match = re.match(r"DECAY\s+(-?\d+)\s+([\d\.E+-]+|NAN).*", line.upper())
andy@238 404 if not match:
andy@238 405 continue
andy@238 406 pdgid = int(match.group(1))
andy@238 407 width = float(match.group(2)) if match.group(2) != "NAN" else None
andy@238 408 currentblock = "DECAY"
andy@238 409 currentdecay = pdgid
andy@238 410 decays[pdgid] = Particle(pdgid, width)
andy@238 411 ## Handle unknown section type start lines (and continue ignoring until a non-header line is found)
andy@238 412 elif type(_autotype(line.split()[0])) is str:
andy@238 413 import sys
andy@238 414 sys.stderr.write("Ignoring unknown section type: %s\n" % line.split()[0])
andy@238 415 currentblock = None
andy@238 416 currentdecay = None
andy@238 417
andy@238 418 ## This non-empty line starts with an indent, hence must be an in-section data line
andy@21 419 else:
andy@238 420 # TODO: Should we now strip the line to remove the indent (and any trailing whitespace)?
andy@21 421 if currentblock is not None:
andy@21 422 items = line.split()
andy@21 423 if len(items) < 1:
andy@6 424 continue
andy@21 425 if currentblock != "DECAY":
andy@200 426 blocks[currentblock].add_entry(items)
andy@238 427 # TODO: Add handling of XSECTION if/when standardised
andy@21 428 else:
andy@185 429 br = float(items[0]) if items[0].upper() != "NAN" else None
andy@21 430 nda = int(items[1])
andy@21 431 ids = map(int, items[2:])
andy@185 432 if br > 0.0 or not ignorenobr: # br == None is < 0
andy@31 433 decays[currentdecay].add_decay(br, nda, ids)
andy@1 434
andy@8 435 ## Try to populate Particle masses from the MASS block
andy@47 436 # print blocks.keys()
andy@47 437 try:
andy@201 438 for pid in blocks["MASS"].keys():
andy@47 439 if decays.has_key(pid):
andy@201 440 decays[pid].mass = blocks["MASS"][pid]
andy@47 441 except:
andy@241 442 if not ignorenomass:
andy@241 443 raise ParseError("No MASS block found: cannot populate particle masses")
andy@8 444
andy@1 445 return blocks, decays
andy@1 446
andy@1 447
andy@124 448
andy@210 449 def writeSLHABlocks(blocks, precision=8):
andy@210 450 """Return an SLHA definition as a string, from the supplied blocks dict."""
andy@211 451 sep = 3 * " "
andy@192 452 blockstrs = []
andy@189 453 for bname, b in blocks.iteritems():
andy@124 454 namestr = b.name
andy@124 455 if b.q is not None:
andy@211 456 namestr += " Q= " + _autostr(float(b.q), precision)
andy@192 457 blockstr = "BLOCK %s\n" % namestr
andy@192 458 entrystrs = []
andy@201 459 for k, v in b.items():
andy@192 460 entrystr = ""
andy@192 461 if type(k) == tuple:
andy@211 462 entrystr += sep.join(_autostr(i, precision) for i in k)
andy@194 463 elif k is not None:
andy@211 464 entrystr += _autostr(k, precision)
andy@211 465 entrystr += sep + _autostr(v, precision)
andy@192 466 entrystrs.append(entrystr)
andy@192 467 blockstr += "\n".join(entrystrs)
andy@192 468 blockstrs.append(blockstr)
andy@210 469 return "\n\n".join(blockstrs)
andy@210 470
andy@210 471
andy@210 472
andy@210 473 def writeSLHADecays(decays, ignorenobr=False, precision=8):
andy@210 474 """Return an SLHA decay definition as a string, from the supplied decays dict."""
andy@211 475 sep = 3 * " "
andy@210 476 blockstrs = []
andy@189 477 for pid, particle in decays.iteritems():
andy@211 478 blockstr = ("DECAY %d " % particle.pid) + _autostr(particle.totalwidth or -1, precision) + "\n"
andy@192 479 decaystrs = []
andy@189 480 for d in particle.decays:
andy@124 481 if d.br > 0.0 or not ignorenobr:
andy@233 482 products_str = sep.join("% d" % i for i in d.ids)
andy@211 483 decaystr = sep + _autostr(d.br, precision) + sep + str(len(d.ids)) + sep + products_str
andy@192 484 decaystrs.append(decaystr)
andy@192 485 blockstr += "\n".join(decaystrs)
andy@192 486 blockstrs.append(blockstr)
andy@192 487 return "\n\n".join(blockstrs)
andy@124 488
andy@124 489
andy@124 490
andy@210 491 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
andy@210 492 """Return an SLHA definition as a string, from the supplied blocks and decays dicts."""
andy@214 493 ss = [x for x in (writeSLHABlocks(blocks, precision), writeSLHADecays(decays, ignorenobr, precision)) if x]
andy@214 494 return "\n\n".join(ss)
andy@210 495
andy@210 496
andy@210 497
andy@124 498 ###############################################################################
andy@161 499 ## PDG <-> HERWIG particle ID code translations for ISAWIG handling
andy@124 500
andy@124 501 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
andy@124 502 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
andy@124 503 ## + the IDPDG array and section 4.13 of the HERWIG manual.
andy@124 504 _HERWIGID2PDGID = {}
andy@124 505 _HERWIGID2PDGID[7] = -1
andy@124 506 _HERWIGID2PDGID[8] = -2
andy@124 507 _HERWIGID2PDGID[9] = -3
andy@124 508 _HERWIGID2PDGID[10] = -4
andy@124 509 _HERWIGID2PDGID[11] = -5
andy@124 510 _HERWIGID2PDGID[12] = -6
andy@124 511 _HERWIGID2PDGID[13] = 21
andy@124 512 _HERWIGID2PDGID[59] = 22
andy@124 513 _HERWIGID2PDGID[121] = 11
andy@124 514 _HERWIGID2PDGID[122] = 12
andy@124 515 _HERWIGID2PDGID[123] = 13
andy@124 516 _HERWIGID2PDGID[124] = 14
andy@124 517 _HERWIGID2PDGID[125] = 15
andy@124 518 _HERWIGID2PDGID[126] = 16
andy@124 519 _HERWIGID2PDGID[127] = -11
andy@124 520 _HERWIGID2PDGID[128] = -12
andy@124 521 _HERWIGID2PDGID[129] = -13
andy@124 522 _HERWIGID2PDGID[130] = -14
andy@124 523 _HERWIGID2PDGID[131] = -15
andy@124 524 _HERWIGID2PDGID[132] = -16
andy@139 525 _HERWIGID2PDGID[198] = 24 # W+
andy@139 526 _HERWIGID2PDGID[199] = -24 # W-
andy@139 527 _HERWIGID2PDGID[200] = 23 # Z0
andy@139 528 _HERWIGID2PDGID[201] = 25 ## SM HIGGS
andy@124 529 _HERWIGID2PDGID[203] = 25 ## HIGGSL0 (== PDG standard in this direction)
andy@124 530 _HERWIGID2PDGID[204] = 35 ## HIGGSH0
andy@124 531 _HERWIGID2PDGID[205] = 36 ## HIGGSA0
andy@124 532 _HERWIGID2PDGID[206] = 37 ## HIGGS+
andy@124 533 _HERWIGID2PDGID[207] = -37 ## HIGGS-
andy@124 534 _HERWIGID2PDGID[401] = 1000001 ## SSDLBR
andy@124 535 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
andy@124 536 _HERWIGID2PDGID[402] = 1000002 ## SSULBR
andy@124 537 _HERWIGID2PDGID[408] = -1000002 ## SSUL
andy@124 538 _HERWIGID2PDGID[403] = 1000003 ## SSSLBR
andy@124 539 _HERWIGID2PDGID[409] = -1000003 ## SSSL
andy@124 540 _HERWIGID2PDGID[404] = 1000004 ## SSCLBR
andy@124 541 _HERWIGID2PDGID[410] = -1000004 ## SSCL
andy@124 542 _HERWIGID2PDGID[405] = 1000005 ## SSB1BR
andy@124 543 _HERWIGID2PDGID[411] = -1000005 ## SSB1
andy@124 544 _HERWIGID2PDGID[406] = 1000006 ## SST1BR
andy@124 545 _HERWIGID2PDGID[412] = -1000006 ## SST1
andy@124 546 _HERWIGID2PDGID[413] = 2000001 ## SSDR
andy@124 547 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
andy@124 548 _HERWIGID2PDGID[414] = 2000002 ## SSUR
andy@124 549 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
andy@124 550 _HERWIGID2PDGID[415] = 2000003 ## SSSR
andy@124 551 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
andy@124 552 _HERWIGID2PDGID[416] = 2000004 ## SSCR
andy@124 553 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
andy@124 554 _HERWIGID2PDGID[417] = 2000005 ## SSB2
andy@124 555 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
andy@124 556 _HERWIGID2PDGID[418] = 2000006 ## SST2
andy@124 557 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
andy@124 558 _HERWIGID2PDGID[425] = 1000011 ## SSEL-
andy@124 559 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
andy@124 560 _HERWIGID2PDGID[426] = 1000012 ## SSNUEL
andy@124 561 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
andy@124 562 _HERWIGID2PDGID[427] = 1000013 ## SSMUL-
andy@124 563 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
andy@124 564 _HERWIGID2PDGID[428] = 1000014 ## SSNUMUL
andy@124 565 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
andy@124 566 _HERWIGID2PDGID[429] = 1000015 ## SSTAU1-
andy@124 567 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
andy@124 568 _HERWIGID2PDGID[430] = 1000016 ## SSNUTL
andy@124 569 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
andy@124 570 _HERWIGID2PDGID[437] = 2000011 ## SSEL-
andy@124 571 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
andy@124 572 _HERWIGID2PDGID[438] = 2000012 ## SSNUEL
andy@124 573 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
andy@124 574 _HERWIGID2PDGID[439] = 2000013 ## SSMUL-
andy@124 575 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
andy@124 576 _HERWIGID2PDGID[440] = 2000014 ## SSNUMUL
andy@124 577 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
andy@124 578 _HERWIGID2PDGID[441] = 2000015 ## SSTAU1-
andy@124 579 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
andy@124 580 _HERWIGID2PDGID[442] = 2000016 ## SSNUTL
andy@124 581 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
andy@124 582 _HERWIGID2PDGID[449] = 1000021 ## GLUINO
andy@124 583 _HERWIGID2PDGID[450] = 1000022 ## NTLINO1
andy@124 584 _HERWIGID2PDGID[451] = 1000023 ## NTLINO2
andy@124 585 _HERWIGID2PDGID[452] = 1000025 ## NTLINO3
andy@124 586 _HERWIGID2PDGID[453] = 1000035 ## NTLINO4
andy@124 587 _HERWIGID2PDGID[454] = 1000024 ## CHGINO1+
andy@124 588 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
andy@124 589 _HERWIGID2PDGID[455] = 1000037 ## CHGINO2+
andy@124 590 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
andy@124 591 _HERWIGID2PDGID[458] = 1000039 ## GRAVTINO
andy@124 592
andy@124 593 def herwigid2pdgid(hwid):
andy@124 594 """
andy@124 595 Convert a particle ID code in the HERWIG internal IDHW format (as used by
andy@124 596 ISAWIG) into its equivalent in the standard PDG ID code definition.
andy@124 597 """
andy@124 598 return _HERWIGID2PDGID.get(hwid, hwid)
andy@124 599
andy@120 600
andy@124 601 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
andy@124 602 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
andy@124 603 ## + the IDPDG array and section 4.13 of the HERWIG manual.
andy@124 604 _PDGID2HERWIGID = {}
andy@124 605 _PDGID2HERWIGID[ -1] = 7
andy@124 606 _PDGID2HERWIGID[ -2] = 8
andy@124 607 _PDGID2HERWIGID[ -3] = 9
andy@124 608 _PDGID2HERWIGID[ -4] = 10
andy@124 609 _PDGID2HERWIGID[ -5] = 11
andy@124 610 _PDGID2HERWIGID[ -6] = 12
andy@124 611 _PDGID2HERWIGID[ 21] = 13
andy@124 612 _PDGID2HERWIGID[ 22] = 59
andy@124 613 _PDGID2HERWIGID[ 11] = 121
andy@124 614 _PDGID2HERWIGID[ 12] = 122
andy@124 615 _PDGID2HERWIGID[ 13] = 123
andy@124 616 _PDGID2HERWIGID[ 14] = 124
andy@124 617 _PDGID2HERWIGID[ 15] = 125
andy@124 618 _PDGID2HERWIGID[ 16] = 126
andy@124 619 _PDGID2HERWIGID[ -11] = 127
andy@124 620 _PDGID2HERWIGID[ -12] = 128
andy@124 621 _PDGID2HERWIGID[ -13] = 129
andy@124 622 _PDGID2HERWIGID[ -14] = 130
andy@124 623 _PDGID2HERWIGID[ -15] = 131
andy@124 624 _PDGID2HERWIGID[ -16] = 132
andy@139 625 _PDGID2HERWIGID[ 24] = 198 ## W+
andy@139 626 _PDGID2HERWIGID[ -24] = 199 ## W-
andy@139 627 _PDGID2HERWIGID[ 23] = 200 ## Z
andy@139 628 _PDGID2HERWIGID[ 25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
andy@124 629 _PDGID2HERWIGID[ 26] = 203 ## HIGGSL0
andy@124 630 _PDGID2HERWIGID[ 35] = 204 ## HIGGSH0
andy@124 631 _PDGID2HERWIGID[ 36] = 205 ## HIGGSA0
andy@124 632 _PDGID2HERWIGID[ 37] = 206 ## HIGGS+
andy@124 633 _PDGID2HERWIGID[ -37] = 207 ## HIGGS-
andy@124 634 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
andy@124 635 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
andy@124 636 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
andy@124 637 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
andy@124 638 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
andy@124 639 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
andy@124 640 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
andy@124 641 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
andy@124 642 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
andy@124 643 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
andy@124 644 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
andy@124 645 _PDGID2HERWIGID[-1000006] = 412 ## SST1
andy@124 646 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
andy@124 647 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
andy@124 648 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
andy@124 649 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
andy@124 650 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
andy@124 651 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
andy@124 652 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
andy@124 653 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
andy@124 654 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
andy@124 655 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
andy@124 656 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
andy@124 657 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
andy@124 658 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
andy@124 659 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
andy@124 660 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
andy@124 661 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
andy@124 662 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
andy@124 663 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
andy@124 664 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
andy@124 665 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
andy@124 666 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
andy@124 667 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
andy@124 668 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
andy@124 669 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
andy@124 670 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
andy@124 671 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
andy@124 672 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
andy@124 673 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
andy@124 674 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
andy@124 675 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
andy@124 676 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
andy@124 677 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
andy@124 678 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
andy@124 679 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
andy@124 680 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
andy@124 681 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
andy@124 682 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
andy@124 683 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
andy@124 684 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
andy@124 685 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
andy@124 686 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
andy@124 687 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
andy@124 688 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
andy@124 689 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
andy@124 690 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
andy@124 691 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
andy@124 692
andy@124 693 def pdgid2herwigid(pdgid):
andy@29 694 """
andy@124 695 Convert a particle ID code in the standard PDG ID code definition into
andy@124 696 its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
andy@124 697 """
andy@124 698 return _PDGID2HERWIGID.get(pdgid, pdgid)
andy@31 699
andy@29 700
andy@161 701 ###############################################################################
andy@161 702 ## ISAWIG format reading/writing
andy@29 703
andy@161 704 def readISAWIG(isastr, ignorenobr=False):
andy@48 705 """
andy@161 706 Read a spectrum definition from a string in the ISAWIG format, returning
andy@161 707 dictionaries of blocks and decays. While this is not an SLHA format, it is
andy@161 708 informally supported as a useful mechanism for converting ISAWIG spectra to
andy@161 709 SLHA.
andy@48 710
andy@161 711 ISAWIG parsing based on the HERWIG SUSY specification format, from
andy@161 712 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
andy@48 713
andy@161 714 If the ignorenobr parameter is True, do not store decay entries with a
andy@161 715 branching ratio of zero.
andy@48 716 """
andy@161 717
andy@189 718 blocks = _mkdict()
andy@189 719 decays = _mkdict()
andy@161 720 LINES = isastr.splitlines()
andy@161 721
andy@161 722 def getnextvalidline():
andy@161 723 while LINES:
andy@161 724 s = LINES.pop(0).strip()
andy@207 725 # print "*", s, "*"
andy@161 726 ## Return None if EOF reached
andy@161 727 if len(s) == 0:
andy@161 728 continue
andy@161 729 ## Strip comments
andy@161 730 if "#" in s:
andy@161 731 s = s[:s.index("#")].strip()
andy@161 732 ## Return if non-empty
andy@161 733 if len(s) > 0:
andy@161 734 return s
andy@161 735
andy@161 736 def getnextvalidlineitems():
andy@161 737 return map(_autotype, getnextvalidline().split())
andy@161 738
andy@161 739 ## Populate MASS block and create decaying particle objects
andy@161 740 masses = Block("MASS")
andy@161 741 numentries = int(getnextvalidline())
andy@161 742 for i in xrange(numentries):
andy@161 743 hwid, mass, lifetime = getnextvalidlineitems()
andy@161 744 width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
andy@161 745 pdgid = herwigid2pdgid(hwid)
andy@202 746 masses[pdgid] = mass
andy@161 747 decays[pdgid] = Particle(pdgid, width, mass)
andy@161 748 #print pdgid, mass, width
andy@161 749 blocks["MASS"] = masses
andy@161 750
andy@161 751 ## Populate decays
andy@161 752 for n in xrange(numentries):
andy@161 753 numdecays = int(getnextvalidline())
andy@161 754 for d in xrange(numdecays):
andy@161 755 #print n, numentries-1, d, numdecays-1
andy@161 756 decayitems = getnextvalidlineitems()
andy@161 757 hwid = decayitems[0]
andy@161 758 pdgid = herwigid2pdgid(hwid)
andy@161 759 br = decayitems[1]
andy@161 760 nme = decayitems[2]
andy@161 761 daughter_hwids = decayitems[3:]
andy@161 762 daughter_pdgids = []
andy@161 763 for hw in daughter_hwids:
andy@161 764 if hw != 0:
andy@161 765 daughter_pdgids.append(herwigid2pdgid(hw))
andy@161 766 if not decays.has_key(pdgid):
andy@161 767 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
andy@161 768 decays[pdgid] = Particle(pdgid)
andy@161 769 decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
andy@161 770
andy@161 771
andy@161 772 ## Now the SUSY parameters
andy@161 773 TANB, ALPHAH = getnextvalidlineitems()
andy@161 774 blocks["MINPAR"] = Block("MINPAR")
andy@202 775 blocks["MINPAR"][3] = TANB
andy@161 776 blocks["ALPHA"] = Block("ALPHA")
andy@222 777 blocks["ALPHA"].set_value(ALPHAH)
andy@161 778 #
andy@161 779 ## Neutralino mixing matrix
andy@161 780 blocks["NMIX"] = Block("NMIX")
andy@161 781 for i in xrange(1, 5):
andy@161 782 nmix_i = getnextvalidlineitems()
andy@161 783 for j, v in enumerate(nmix_i):
andy@202 784 blocks["NMIX"][i, j+1] = v
andy@161 785 #
andy@161 786 ## Chargino mixing matrices V and U
andy@161 787 blocks["VMIX"] = Block("VMIX")
andy@161 788 vmix = getnextvalidlineitems()
andy@202 789 blocks["VMIX"][1, 1] = vmix[0]
andy@202 790 blocks["VMIX"][1, 2] = vmix[1]
andy@202 791 blocks["VMIX"][2, 1] = vmix[2]
andy@202 792 blocks["VMIX"][2, 2] = vmix[3]
andy@161 793 blocks["UMIX"] = Block("UMIX")
andy@161 794 umix = getnextvalidlineitems()
andy@202 795 blocks["UMIX"][1, 1] = umix[0]
andy@202 796 blocks["UMIX"][1, 2] = umix[1]
andy@202 797 blocks["UMIX"][2, 1] = umix[2]
andy@202 798 blocks["UMIX"][2, 2] = umix[3]
andy@161 799 #
andy@161 800 THETAT, THETAB, THETAL = getnextvalidlineitems()
andy@161 801 import math
andy@161 802 blocks["STOPMIX"] = Block("STOPMIX")
andy@202 803 blocks["STOPMIX"][1, 1] = math.cos(THETAT)
andy@202 804 blocks["STOPMIX"][1, 2] = -math.sin(THETAT)
andy@202 805 blocks["STOPMIX"][2, 1] = math.sin(THETAT)
andy@202 806 blocks["STOPMIX"][2, 2] = math.cos(THETAT)
andy@161 807 blocks["SBOTMIX"] = Block("SBOTMIX")
andy@202 808 blocks["SBOTMIX"][1, 1] = math.cos(THETAB)
andy@202 809 blocks["SBOTMIX"][1, 2] = -math.sin(THETAB)
andy@202 810 blocks["SBOTMIX"][2, 1] = math.sin(THETAB)
andy@202 811 blocks["SBOTMIX"][2, 2] = math.cos(THETAB)
andy@161 812 blocks["STAUMIX"] = Block("STAUMIX")
andy@202 813 blocks["STAUMIX"][1, 1] = math.cos(THETAL)
andy@202 814 blocks["STAUMIX"][1, 2] = -math.sin(THETAL)
andy@202 815 blocks["STAUMIX"][2, 1] = math.sin(THETAL)
andy@202 816 blocks["STAUMIX"][2, 2] = math.cos(THETAL)
andy@161 817 #
andy@161 818 ATSS, ABSS, ALSS = getnextvalidlineitems()
andy@161 819 blocks["AU"] = Block("AU")
andy@205 820 blocks["AU"][3, 3] = ATSS
andy@161 821 blocks["AD"] = Block("AD")
andy@205 822 blocks["AD"][3, 3] = ABSS
andy@161 823 blocks["AE"] = Block("AE")
andy@205 824 blocks["AE"][3, 3] = ALSS
andy@161 825 #
andy@161 826 MUSS = getnextvalidlineitems()[0]
andy@202 827 blocks["MINPAR"][4] = MUSS
andy@161 828 #
andy@161 829
andy@161 830 # TODO: Parse RPV boolean and couplings into SLHA2 blocks
andy@161 831
andy@161 832 return blocks, decays
andy@161 833
andy@161 834
andy@147 835 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
andy@48 836 """
andy@161 837 Return a SUSY spectrum definition in the format produced by ISAWIG for inut to HERWIG
andy@161 838 as a string, from the supplied SLHA blocks and decays dicts.
andy@48 839
andy@48 840 ISAWIG parsing based on the HERWIG SUSY specification format, from
andy@48 841 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
andy@48 842
andy@48 843 If the ignorenobr parameter is True, do not write decay entries with a
andy@48 844 branching ratio of zero.
andy@48 845 """
andy@201 846 masses = blocks["MASS"]
andy@48 847
andy@48 848 ## Init output string
andy@48 849 out = ""
andy@48 850
andy@48 851 ## First write out masses section:
andy@48 852 ## Number of SUSY + top particles
andy@48 853 ## IDHW, RMASS(IDHW), RLTIM(IDHW)
andy@48 854 ## repeated for each particle
andy@48 855 ## IDHW is the HERWIG identity code.
andy@48 856 ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
andy@48 857 massout = ""
andy@48 858 for pid in masses.keys():
andy@48 859 lifetime = -1
andy@48 860 try:
andy@48 861 width = decays[pid].totalwidth
andy@48 862 if width and width > 0:
andy@48 863 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
andy@48 864 except:
andy@48 865 pass
andy@211 866 massout += ("%d " % pdgid2herwigid(pid)) + _autostr(masses[pid], precision) + " " + _autostr(lifetime, precision) + "\n"
andy@48 867 out += "%d\n" % massout.count("\n")
andy@48 868 out += massout
andy@48 869
andy@48 870 assert(len(masses) == len(decays))
andy@48 871
andy@48 872 ## Next each particles decay modes together with their branching ratios and matrix element codes
andy@48 873 ## Number of decay modes for a given particle (IDK)
andy@48 874 ## IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
andy@48 875 ## repeated for each mode.
andy@48 876 ## Repeated for each particle.
andy@48 877 ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
andy@48 878 ## the decay mode. NME is a code for the matrix element to be used, either from the
andy@48 879 ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
andy@48 880 for i, pid in enumerate(decays.keys()):
andy@48 881 # if not decays.has_key(pid):
andy@48 882 # continue
andy@124 883 hwid = pdgid2herwigid(pid)
andy@48 884 decayout = ""
andy@48 885 #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
andy@48 886 for i_d, d in enumerate(decays[pid].decays):
andy@48 887 ## Skip decay if it has no branching ratio
andy@48 888 if ignorenobr and d.br == 0:
andy@48 889 continue
andy@71 890
andy@71 891 ## Identify decay matrix element to use
andy@59 892 ## From std HW docs, or from this pair:
andy@59 893 ## Two new matrix element codes have been added for these new decays:
andy@59 894 ## NME = 200 3 body top quark via charged Higgs
andy@59 895 ## 300 3 body R-parity violating gaugino and gluino decays
andy@71 896 nme = 0
andy@71 897 # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
andy@71 898 if abs(pid) in (6, 12):
andy@71 899 nme = 100
andy@71 900 ## Extra SUSY MEs
andy@71 901 if len(d.ids) == 3:
andy@71 902 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
andy@71 903 pass
andy@211 904 decayout += ("%d " % hwid) + _autostr(d.br, precision) + (" %d " % nme)
andy@71 905
andy@71 906 def is_quark(pid):
andy@71 907 return (abs(pid) in range(1, 7))
andy@71 908
andy@71 909 def is_lepton(pid):
andy@71 910 return (abs(pid) in range(11, 17))
andy@71 911
andy@71 912 def is_squark(pid):
andy@71 913 if abs(pid) in range(1000001, 1000007):
andy@71 914 return True
andy@71 915 if abs(pid) in range(2000001, 2000007):
andy@71 916 return True
andy@71 917 return False
andy@71 918
andy@71 919 def is_slepton(pid):
andy@71 920 if abs(pid) in range(1000011, 1000017):
andy@71 921 return True
andy@71 922 if abs(pid) in range(2000011, 2000016, 2):
andy@71 923 return True
andy@71 924 return False
andy@71 925
andy@71 926 def is_gaugino(pid):
andy@71 927 if abs(pid) in range(1000022, 1000026):
andy@71 928 return True
andy@71 929 if abs(pid) in (1000035, 1000037):
andy@71 930 return True
andy@71 931 return False
andy@71 932
andy@71 933 def is_susy(pid):
andy@71 934 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
andy@71 935
andy@71 936 absids = map(abs, d.ids)
andy@71 937
andy@66 938 ## Order decay products as required by HERWIG
andy@66 939 ## Top
andy@66 940 if abs(pid) == 6:
andy@66 941 def cmp_bottomlast(a, b):
andy@66 942 """Comparison function which always puts b/bbar last"""
andy@66 943 if abs(a) == 5:
andy@66 944 return True
andy@71 945 if abs(b) == 5:
andy@66 946 return False
andy@66 947 return cmp(a, b)
andy@66 948 if len(absids) == 2:
andy@66 949 ## 2 body mode, to Higgs: Higgs; Bottom
andy@66 950 if (25 in absids or 26 in absids) and 5 in absids:
andy@199 951 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
andy@66 952 elif len(absids) == 3:
andy@66 953 ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
andy@66 954 if 37 in absids or 23 in absids:
andy@199 955 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
andy@66 956 ## Gluino
andy@66 957 elif abs(pid) == 1000021:
andy@66 958 if len(absids) == 2:
andy@66 959 ## 2 body mode
andy@66 960 ## without gluon: any order
andy@66 961 ## with gluon: gluon; colour neutral
andy@66 962 if 21 in absids:
andy@66 963 def cmp_gluonfirst(a, b):
andy@66 964 """Comparison function which always puts gluon first"""
andy@66 965 if a == 21:
andy@66 966 return False
andy@71 967 if b == 21:
andy@66 968 return True
andy@66 969 return cmp(a, b)
andy@199 970 d.ids = sorted(d.ids, cmp=cmp_gluonfirst)
andy@66 971 elif len(absids) == 3:
andy@66 972 ## 3-body modes, R-parity conserved: colour neutral; q or qbar
andy@70 973 def cmp_quarkslast(a, b):
andy@70 974 """Comparison function which always puts quarks last"""
andy@71 975 if is_quark(a):
andy@70 976 return True
andy@71 977 if is_quark(b):
andy@71 978 return False
andy@71 979 return cmp(a, b)
andy@199 980 d.ids = sorted(d.ids, cmp=cmp_quarkslast)
andy@71 981 ## Squark/Slepton
andy@71 982 elif is_squark(pid) or is_slepton(pid):
andy@73 983 def cmp_susy_quark_lepton(a, b):
andy@71 984 if is_susy(a):
andy@71 985 return False
andy@71 986 if is_susy(b):
andy@71 987 return True
andy@71 988 if is_quark(a):
andy@71 989 return False
andy@71 990 if is_quark(b):
andy@71 991 return True
andy@71 992 return cmp(a, b)
andy@71 993 ## 2 body modes: Gaugino/Gluino with Quark/Lepton Gaugino quark
andy@71 994 ## Gluino lepton
andy@71 995 ## 3 body modes: Weak sparticle particles from W decay
andy@71 996 ## Squark
andy@71 997 ## 2 body modes: Lepton Number Violated quark lepton
andy@71 998 ## Baryon number violated quark quark
andy@71 999 ## Slepton
andy@71 1000 ## 2 body modes: Lepton Number Violated q or qbar
andy@199 1001 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
andy@71 1002 ## Higgs
andy@71 1003 elif pid in (25, 26):
andy@71 1004 # TODO: Includes SUSY Higgses?
andy@71 1005 ## Higgs
andy@71 1006 ## 2 body modes: (s)quark-(s)qbar (s)q or (s)qbar
andy@71 1007 ## (s)lepton-(s)lepton (s)l or (s)lbar
andy@71 1008 ## 3 body modes: colour neutral q or qbar
andy@71 1009 if len(absids) == 3:
andy@71 1010 def cmp_quarkslast(a, b):
andy@71 1011 """Comparison function which always puts quarks last"""
andy@71 1012 if is_quark(a):
andy@71 1013 return True
andy@71 1014 if is_quark(b):
andy@71 1015 return False
andy@71 1016 return cmp(a, b)
andy@199 1017 d.ids = sorted(d.ids, cmp=cmp_quarkslast)
andy@71 1018 elif is_gaugino(pid):
andy@71 1019 # TODO: Is there actually anything to do here?
andy@71 1020 ## Gaugino
andy@71 1021 ## 2 body modes: Squark-quark q or sq
andy@71 1022 ## Slepton-lepton l or sl
andy@71 1023 ##
andy@71 1024 ## 3 body modes: R-parity conserved colour neutral q or qbar
andy@71 1025 ## l or lbar
andy@71 1026 if len(absids) == 3:
andy@71 1027 def cmp_quarkslast(a, b):
andy@71 1028 """Comparison function which always puts quarks last"""
andy@71 1029 if is_quark(a):
andy@71 1030 return True
andy@71 1031 if is_quark(b):
andy@70 1032 return False
andy@70 1033 return cmp(a, b)
andy@199 1034 d.ids = sorted(d.ids, cmp=cmp_quarkslast)
andy@66 1035
andy@71 1036 # TODO: Gaugino/Gluino
andy@67 1037 ## 3 body modes: R-parity violating: Particles in the same order as the R-parity violating superpotential
andy@66 1038
andy@66 1039 ## Pad out IDs list with zeros
andy@48 1040 ids = [0,0,0,0,0]
andy@48 1041 for i, pid in enumerate(d.ids):
andy@48 1042 ids[i] = pid
andy@48 1043 ids = map(str, ids)
andy@48 1044 decayout += " ".join(ids) + "\n"
andy@48 1045 decayout = "%d\n" % decayout.count("\n") + decayout
andy@48 1046 out += decayout
andy@48 1047
andy@48 1048 ## Now the SUSY parameters
andy@48 1049 ## TANB, ALPHAH:
andy@211 1050 out += _autostr(blocks["MINPAR"][3], precision) + " " + _autostr(blocks["ALPHA"].value(), precision) + "\n"
andy@48 1051 ## Neutralino mixing matrix
andy@200 1052 nmix = blocks["NMIX"]
andy@48 1053 for i in xrange(1, 5):
andy@211 1054 out += _autostr(nmix[i,1], precision) + " " + \
andy@211 1055 _autostr(nmix[i,2], precision) + " " + \
andy@211 1056 _autostr(nmix[i,3], precision) + " " + \
andy@211 1057 _autostr(nmix[i,4], precision) + "\n"
andy@48 1058 ## Chargino mixing matrices V and U
andy@200 1059 vmix = blocks["VMIX"]
andy@211 1060 out += _autostr(vmix[1,1], precision) + " " + \
andy@211 1061 _autostr(vmix[1,2], precision) + " " + \
andy@211 1062 _autostr(vmix[2,1], precision) + " " + \
andy@211 1063 _autostr(vmix[2,2], precision) + "\n"
andy@200 1064 umix = blocks["UMIX"]
andy@211 1065 out += _autostr(umix[1,1], precision) + " " + \
andy@211 1066 _autostr(umix[1,2], precision) + " " + \
andy@211 1067 _autostr(umix[2,1], precision) + " " + \
andy@211 1068 _autostr(umix[2,2], precision) + "\n"
andy@160 1069 ## THETAT,THETAB,THETAL
andy@48 1070 import math
andy@211 1071 out += _autostr(math.acos(blocks["STOPMIX"][1,1]), precision) + " " + \
andy@211 1072 _autostr(math.acos(blocks["SBOTMIX"][1,1]), precision) + " " + \
andy@211 1073 _autostr(math.acos(blocks["STAUMIX"][1,1]), precision) + "\n"
andy@160 1074 ## ATSS,ABSS,ALSS
andy@211 1075 out += _autostr(blocks["AU"][3,3], precision) + " " + \
andy@211 1076 _autostr(blocks["AD"][3,3], precision) + " " + \
andy@211 1077 _autostr(blocks["AE"][3,3], precision) + "\n"
andy@160 1078 ## MUSS == sign(mu)
andy@201 1079 out += "%f\n" % blocks["MINPAR"][4]
andy@48 1080
andy@160 1081 ## RPV SUSY
andy@160 1082 isRPV = False
andy@160 1083 out += "%d\n" % isRPV
andy@160 1084 # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
andy@160 1085 # TODO: Get the index orderings right
andy@160 1086 # if isRPV: ...
andy@48 1087
andy@48 1088 return out
andy@48 1089
andy@48 1090
andy@223 1091
andy@161 1092 ###############################################################################
andy@223 1093 ## File-level read/write functions
andy@161 1094
andy@161 1095 def readSLHAFile(spcfilename, **kwargs):
andy@161 1096 """
andy@161 1097 Read an SLHA file, returning dictionaries of blocks and decays.
andy@161 1098
andy@161 1099 Other keyword parameters are passed to readSLHA.
andy@161 1100 """
andy@161 1101 f = open(spcfilename, "r")
andy@240 1102 rtn = readSLHA(f.read(), **kwargs)
andy@161 1103 f.close()
andy@161 1104 return rtn
andy@161 1105
andy@161 1106
andy@161 1107 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
andy@161 1108 """
andy@161 1109 Write an SLHA file from the supplied blocks and decays dicts.
andy@161 1110
andy@161 1111 Other keyword parameters are passed to writeSLHA.
andy@161 1112 """
andy@161 1113 f = open(spcfilename, "w")
andy@240 1114 f.write(writeSLHA(blocks, decays, **kwargs))
andy@161 1115 f.close()
andy@161 1116
andy@161 1117
andy@161 1118 def readISAWIGFile(isafilename, **kwargs):
andy@161 1119 """
andy@161 1120 Read a spectrum definition from a file in the ISAWIG format, returning
andy@161 1121 dictionaries of blocks and decays. While this is not an SLHA format, it is
andy@161 1122 informally supported as a useful mechanism for converting ISAWIG spectra to
andy@161 1123 SLHA.
andy@161 1124
andy@161 1125 Other keyword parameters are passed to readSLHA.
andy@161 1126 """
andy@161 1127 f = open(isafilename, "r")
andy@240 1128 rtn = readISAWIG(f.read(), **kwargs)
andy@161 1129 f.close()
andy@161 1130 return rtn
andy@161 1131
andy@161 1132
andy@161 1133 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
andy@161 1134 """
andy@161 1135 Write an ISAWIG file from the supplied blocks and decays dicts.
andy@161 1136
andy@161 1137 Other keyword parameters are passed to writeISAWIG.
andy@161 1138 """
andy@161 1139 f = open(isafilename, "w")
andy@240 1140 f.write(writeISAWIG(blocks, decays, **kwargs))
andy@161 1141 f.close()
andy@161 1142
andy@161 1143
andy@161 1144
andy@161 1145 ###############################################################################
andy@161 1146 ## Main function for module testing
andy@161 1147
andy@1 1148 if __name__ == "__main__":
andy@1 1149 import sys
andy@1 1150 for a in sys.argv[1:]:
andy@35 1151 if a.endswith(".isa"):
andy@35 1152 blocks, decays = readISAWIGFile(a)
andy@35 1153 else:
andy@35 1154 blocks, decays = readSLHAFile(a)
andy@3 1155
andy@5 1156 for bname, b in sorted(blocks.iteritems()):
andy@5 1157 print b
andy@5 1158 print
andy@3 1159
andy@47 1160 print blocks.keys()
andy@47 1161
andy@201 1162 print blocks["MASS"][25]
andy@6 1163 print
andy@6 1164
andy@6 1165 for p in sorted(decays.values()):
andy@6 1166 print p
andy@6 1167 print
andy@29 1168
andy@31 1169 print writeSLHA(blocks, decays, ignorenobr=True)

mercurial