pyslha.py

Sun, 28 Apr 2013 22:13:19 +0200

author
Andy Buckley <andy@insectnation.org>
date
Sun, 28 Apr 2013 22:13:19 +0200
changeset 203
ce90a0dace07
parent 202
d3d069f4549a
child 204
bef82eaef56e
permissions
-rw-r--r--

Adding Block.set_value(*args) and Block documentation.

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

mercurial