pyslha.py

Mon, 04 Oct 2010 09:29:35 +0100

author
Andy Buckley <andy@insectnation.org>
date
Mon, 04 Oct 2010 09:29:35 +0100
changeset 67
2733e0ddb7f9
parent 66
7c829ecb90e1
child 68
03bafdeeffa2
permissions
-rw-r--r--

Tidying sorting table

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@52 11 The current release supports SLHA version 1. Assistance with supporting version
andy@52 12 2 will be gladly accepted!
andy@52 13
andy@64 14 The plotting script currently provides output in a source form usable by the
andy@64 15 Rivet make-plots system, but alternative formats are in the pipeline. Again,
andy@64 16 contributions or suggestions of new formats are very welcome.
andy@64 17
andy@52 18 TODOs:
andy@64 19 * Identify HERWIG decay matrix element to use in ISAWIG
andy@64 20 * Order decay products as required by HERWIG in ISAWIG
andy@52 21 * What to do about unlisted PIDs?
andy@52 22 * Split writeSLHA into writeSLHA{Blocks,Decays}
andy@52 23 * Handle SLHA2
andy@52 24 * Handle RPV SUSY in ISAWIG
andy@8 25 """
andy@8 26
andy@8 27 from __future__ import with_statement
andy@26 28
andy@8 29 __author__ = "Andy Buckley <andy.buckley@cern.ch"
andy@64 30 __version__ = "0.4.0"
andy@8 31
andy@1 32
andy@4 33 def _autotype(var):
andy@30 34 """Automatically convert strings to numerical types if possible."""
andy@30 35 if type(var) is not str:
andy@4 36 return var
andy@36 37 if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
andy@4 38 return int(var)
andy@4 39 try:
andy@4 40 f = float(var)
andy@4 41 return f
andy@4 42 except ValueError:
andy@4 43 return var
andy@4 44
andy@30 45 def _autostr(var):
andy@30 46 """Automatically numerical types to the right sort of string."""
andy@30 47 if type(var) is float:
andy@30 48 return "%e" % var
andy@30 49 return str(var)
andy@30 50
andy@30 51
andy@4 52
andy@12 53 class Block(object):
andy@8 54 """
andy@8 55 Object representation of any BLOCK elements read from the SLHA file. Blocks
andy@8 56 have a name, may have an associated Q value, and then a collection of data
andy@8 57 entries, stored as a recursive dictionary. Types in the dictionary are
andy@8 58 numeric (int or float) when a cast from the string in the file has been
andy@8 59 possible.
andy@8 60 """
andy@5 61 def __init__(self, name, q=None):
andy@1 62 self.name = name
andy@1 63 self.entries = {}
andy@5 64 self.q = _autotype(q)
andy@1 65
andy@1 66 def add_entry(self, entry):
andy@1 67 #print entry
andy@1 68 nextparent = self.entries
andy@1 69 if len(entry) < 2:
andy@1 70 raise Exception("Block entries must be at least a 2-tuple")
andy@4 71 #print "in", entry
andy@4 72 entry = map(_autotype, entry)
andy@4 73 #print "out", entry
andy@1 74 for e in entry[:-2]:
andy@1 75 if e is not entry[-1]:
andy@1 76 nextparent = nextparent.setdefault(e, {})
andy@1 77 nextparent[entry[-2]] = entry[-1]
andy@1 78 #print self.entries
andy@1 79
andy@1 80 def __cmp__(self, other):
andy@31 81 return cmp(self.name, other.name)
andy@1 82
andy@1 83 def __str__(self):
andy@1 84 s = self.name
andy@5 85 if self.q is not None:
andy@5 86 s += " (Q=%s)" % self.q
andy@1 87 s += "\n"
andy@1 88 s += str(self.entries)
andy@1 89 return s
andy@1 90
andy@1 91
andy@12 92 class Decay(object):
andy@8 93 """
andy@8 94 Object representing a decay entry on a particle decribed by the SLHA file.
andy@8 95 'Decay' objects are not a direct representation of a DECAY block in an SLHA
andy@8 96 file... that role, somewhat confusingly, is taken by the Particle class.
andy@8 97
andy@8 98 Decay objects have three properties: a branching ratio, br, an nda number
andy@12 99 (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
andy@12 100 decay occurs. The PDG ID of the particle whose decay this represents may
andy@12 101 also be stored, but this is normally known via the Particle in which the
andy@12 102 decay is stored.
andy@8 103 """
andy@8 104 def __init__(self, br, nda, ids, parentid=None):
andy@8 105 self.parentid = parentid
andy@6 106 self.br = br
andy@6 107 self.nda = nda
andy@6 108 self.ids = ids
andy@29 109 assert(self.nda == len(self.ids))
andy@6 110
andy@6 111 def __cmp__(self, other):
andy@31 112 return cmp(other.br, self.br)
andy@6 113
andy@6 114 def __str__(self):
andy@12 115 return "%e %s" % (self.br, self.ids)
andy@6 116
andy@6 117
andy@12 118 class Particle(object):
andy@8 119 """
andy@8 120 Representation of a single, specific particle, decay block from an SLHA
andy@8 121 file. These objects are not themselves called 'Decay', since that concept
andy@8 122 applies more naturally to the various decays found inside this
andy@8 123 object. Particle classes store the PDG ID (pid) of the particle being
andy@8 124 represented, and optionally the mass (mass) and total decay width
andy@8 125 (totalwidth) of that particle in the SLHA scenario. Masses may also be found
andy@8 126 via the MASS block, from which the Particle.mass property is filled, if at
andy@8 127 all. They also store a list of Decay objects (decays) which are probably the
andy@8 128 item of most interest.
andy@8 129 """
andy@6 130 def __init__(self, pid, totalwidth=None, mass=None):
andy@6 131 self.pid = pid
andy@6 132 self.totalwidth = totalwidth
andy@6 133 self.mass = mass
andy@6 134 self.decays = []
andy@6 135
andy@6 136 def add_decay(self, br, nda, ids):
andy@6 137 self.decays.append(Decay(br, nda, ids))
andy@6 138 self.decays.sort()
andy@6 139
andy@6 140 def __cmp__(self, other):
andy@6 141 if abs(self.pid) == abs(other.pid):
andy@31 142 return cmp(self.pid, other.pid)
andy@31 143 return cmp(abs(self.pid), abs(other.pid))
andy@6 144
andy@6 145 def __str__(self):
andy@6 146 s = str(self.pid)
andy@7 147 if self.mass is not None:
andy@7 148 s += " : mass = %e GeV" % self.mass
andy@6 149 if self.totalwidth is not None:
andy@7 150 s += " : total width = %e GeV" % self.totalwidth
andy@6 151 for d in self.decays:
andy@12 152 if d.br > 0.0:
andy@12 153 s += "\n %s" % d
andy@6 154 return s
andy@1 155
andy@1 156
andy@31 157 def readSLHAFile(spcfilename, **kwargs):
andy@21 158 """
andy@21 159 Read an SLHA file, returning dictionaries of blocks and decays.
andy@31 160
andy@31 161 Other keyword parameters are passed to readSLHA.
andy@21 162 """
andy@21 163 with open(spcfilename, "r") as f:
andy@31 164 return readSLHA(f.read(), kwargs)
andy@21 165
andy@21 166
andy@31 167 def readSLHA(spcstr, ignorenobr=False):
andy@21 168 """
andy@31 169 Read an SLHA definition from a string, returning dictionaries of blocks and
andy@31 170 decays.
andy@31 171
andy@31 172 If the ignorenobr parameter is True, do not store decay entries with a
andy@31 173 branching ratio of zero.
andy@21 174 """
andy@1 175 blocks = {}
andy@1 176 decays = {}
andy@21 177 #
andy@34 178 import re
andy@21 179 currentblock = None
andy@21 180 currentdecay = None
andy@21 181 for line in spcstr.splitlines():
andy@21 182 ## Handle (ignore) comment lines
andy@21 183 if line.startswith("#"):
andy@21 184 continue
andy@21 185 if "#" in line:
andy@21 186 line = line[:line.index("#")]
andy@21 187
andy@21 188 ## Handle BLOCK/DECAY start lines
andy@21 189 if line.upper().startswith("BLOCK"):
andy@47 190 #print line
andy@47 191 match = re.match(r"BLOCK\s+(\w+)(\s+Q=\s*.+)?", line.upper())
andy@21 192 if not match:
andy@8 193 continue
andy@21 194 blockname = match.group(1)
andy@21 195 qstr = match.group(2)
andy@21 196 if qstr is not None:
andy@21 197 qstr = qstr[2:].strip()
andy@21 198 currentblock = blockname
andy@21 199 currentdecay = None
andy@21 200 blocks[blockname] = Block(blockname, q=qstr)
andy@21 201 elif line.upper().startswith("DECAY"):
andy@21 202 match = re.match(r"DECAY\s+(\d+)\s+([\d\.E+-]+).*", line.upper())
andy@21 203 if not match:
andy@21 204 continue
andy@21 205 pdgid = int(match.group(1))
andy@21 206 width = float(match.group(2))
andy@21 207 currentblock = "DECAY"
andy@21 208 currentdecay = pdgid
andy@21 209 decays[pdgid] = Particle(pdgid, width)
andy@21 210 else:
andy@21 211 ## In-block line
andy@21 212 if currentblock is not None:
andy@21 213 items = line.split()
andy@21 214 if len(items) < 1:
andy@6 215 continue
andy@21 216 if currentblock != "DECAY":
andy@21 217 if len(items) < 2:
andy@21 218 ## Treat the ALPHA block differently
andy@21 219 blocks[currentblock].value = _autotype(items[0])
andy@33 220 blocks[currentblock].entries = _autotype(items[0])
andy@8 221 else:
andy@21 222 blocks[currentblock].add_entry(items)
andy@21 223 else:
andy@21 224 br = float(items[0])
andy@21 225 nda = int(items[1])
andy@21 226 ids = map(int, items[2:])
andy@31 227 if br > 0.0 or not ignorenobr:
andy@31 228 decays[currentdecay].add_decay(br, nda, ids)
andy@1 229
andy@8 230 ## Try to populate Particle masses from the MASS block
andy@47 231 # print blocks.keys()
andy@47 232 try:
andy@47 233 for pid in blocks["MASS"].entries.keys():
andy@47 234 if decays.has_key(pid):
andy@47 235 decays[pid].mass = blocks["MASS"].entries[pid]
andy@47 236 except:
andy@47 237 raise Exception("No MASS block found, from which to populate particle masses")
andy@8 238
andy@1 239 return blocks, decays
andy@1 240
andy@1 241
andy@33 242 def readISAWIGFile(isafilename, **kwargs):
andy@33 243 """
andy@33 244 Read a spectrum definition from a file in the ISAWIG format, returning
andy@33 245 dictionaries of blocks and decays. While this is not an SLHA format, it is
andy@33 246 informally supported as a useful mechanism for converting ISAWIG spectra to
andy@33 247 SLHA.
andy@33 248
andy@33 249 Other keyword parameters are passed to readSLHA.
andy@33 250 """
andy@33 251 with open(isafilename, "r") as f:
andy@33 252 return readISAWIG(f.read(), kwargs)
andy@33 253
andy@33 254
andy@33 255 def readISAWIG(isastr, ignorenobr=False):
andy@33 256 """
andy@33 257 Read a spectrum definition from a string in the ISAWIG format, returning
andy@33 258 dictionaries of blocks and decays. While this is not an SLHA format, it is
andy@33 259 informally supported as a useful mechanism for converting ISAWIG spectra to
andy@33 260 SLHA.
andy@33 261
andy@33 262 ISAWIG parsing based on the HERWIG SUSY specification format, from
andy@33 263 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
andy@33 264
andy@33 265 If the ignorenobr parameter is True, do not store decay entries with a
andy@33 266 branching ratio of zero.
andy@33 267 """
andy@33 268
andy@33 269 ## PDG MC ID codes mapped to HERWIG SUSY ID codes, based on
andy@33 270 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
andy@33 271 HERWIGID2PDGID = {}
andy@33 272 HERWIGID2PDGID[203] = 25 ## HIGGSL0
andy@33 273 HERWIGID2PDGID[204] = 35 ## HIGGSH0
andy@33 274 HERWIGID2PDGID[205] = 36 ## HIGGSA0
andy@33 275 HERWIGID2PDGID[206] = 37 ## HIGGS+
andy@33 276 HERWIGID2PDGID[207] = -37 ## HIGGS-
andy@33 277 HERWIGID2PDGID[401] = 1000001 ## SSDLBR
andy@33 278 HERWIGID2PDGID[407] = -1000001 ## SSDLBR
andy@33 279 HERWIGID2PDGID[402] = 1000002 ## SSULBR
andy@33 280 HERWIGID2PDGID[408] = -1000002 ## SSUL
andy@33 281 HERWIGID2PDGID[403] = 1000003 ## SSSLBR
andy@33 282 HERWIGID2PDGID[409] = -1000003 ## SSSL
andy@33 283 HERWIGID2PDGID[404] = 1000004 ## SSCLBR
andy@33 284 HERWIGID2PDGID[410] = -1000004 ## SSCL
andy@33 285 HERWIGID2PDGID[405] = 1000005 ## SSB1BR
andy@33 286 HERWIGID2PDGID[411] = -1000005 ## SSB1
andy@33 287 HERWIGID2PDGID[406] = 1000006 ## SST1BR
andy@33 288 HERWIGID2PDGID[412] = -1000006 ## SST1
andy@33 289 HERWIGID2PDGID[413] = 2000001 ## SSDR
andy@33 290 HERWIGID2PDGID[419] = -2000001 ## SSDRBR
andy@33 291 HERWIGID2PDGID[414] = 2000002 ## SSUR
andy@33 292 HERWIGID2PDGID[420] = -2000002 ## SSURBR
andy@33 293 HERWIGID2PDGID[415] = 2000003 ## SSSR
andy@33 294 HERWIGID2PDGID[421] = -2000003 ## SSSRBR
andy@33 295 HERWIGID2PDGID[416] = 2000004 ## SSCR
andy@33 296 HERWIGID2PDGID[422] = -2000004 ## SSCRBR
andy@33 297 HERWIGID2PDGID[417] = 2000005 ## SSB2
andy@33 298 HERWIGID2PDGID[423] = -2000005 ## SSB2BR
andy@33 299 HERWIGID2PDGID[418] = 2000006 ## SST2
andy@33 300 HERWIGID2PDGID[424] = -2000006 ## SST2BR
andy@33 301 HERWIGID2PDGID[425] = 1000011 ## SSEL-
andy@33 302 HERWIGID2PDGID[431] = -1000011 ## SSEL+
andy@33 303 HERWIGID2PDGID[426] = 1000012 ## SSNUEL
andy@33 304 HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
andy@33 305 HERWIGID2PDGID[427] = 1000013 ## SSMUL-
andy@33 306 HERWIGID2PDGID[433] = -1000013 ## SSMUL+
andy@33 307 HERWIGID2PDGID[428] = 1000014 ## SSNUMUL
andy@33 308 HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
andy@33 309 HERWIGID2PDGID[429] = 1000015 ## SSTAU1-
andy@33 310 HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
andy@33 311 HERWIGID2PDGID[430] = 1000016 ## SSNUTL
andy@33 312 HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
andy@33 313 HERWIGID2PDGID[437] = 2000011 ## SSEL-
andy@33 314 HERWIGID2PDGID[443] = -2000011 ## SSEL+
andy@33 315 HERWIGID2PDGID[438] = 2000012 ## SSNUEL
andy@33 316 HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
andy@33 317 HERWIGID2PDGID[439] = 2000013 ## SSMUL-
andy@33 318 HERWIGID2PDGID[445] = -2000013 ## SSMUL+
andy@33 319 HERWIGID2PDGID[440] = 2000014 ## SSNUMUL
andy@33 320 HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
andy@33 321 HERWIGID2PDGID[441] = 2000015 ## SSTAU1-
andy@33 322 HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
andy@33 323 HERWIGID2PDGID[442] = 2000016 ## SSNUTL
andy@33 324 HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
andy@33 325 HERWIGID2PDGID[449] = 1000021 ## GLUINO
andy@33 326 HERWIGID2PDGID[450] = 1000022 ## NTLINO1
andy@33 327 HERWIGID2PDGID[451] = 1000023 ## NTLINO2
andy@33 328 HERWIGID2PDGID[452] = 1000025 ## NTLINO3
andy@33 329 HERWIGID2PDGID[453] = 1000035 ## NTLINO4
andy@33 330 HERWIGID2PDGID[454] = 1000024 ## CHGINO1+
andy@33 331 HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
andy@33 332 HERWIGID2PDGID[455] = 1000037 ## CHGINO2+
andy@33 333 HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
andy@33 334 HERWIGID2PDGID[458] = 1000039 ## GRAVTINO
andy@33 335
andy@33 336 blocks = {}
andy@33 337 decays = {}
andy@35 338 LINES = isastr.splitlines()
andy@33 339
andy@33 340 def getnextvalidline():
andy@35 341 while LINES:
andy@35 342 s = LINES.pop(0).strip()
andy@33 343 ## Return None if EOF reached
andy@33 344 if len(s) == 0:
andy@33 345 continue
andy@33 346 ## Strip comments
andy@33 347 if "#" in s:
andy@33 348 s = s[:s.index("#")].strip()
andy@33 349 ## Return if non-empty
andy@33 350 if len(s) > 0:
andy@33 351 return s
andy@33 352
andy@34 353 def getnextvalidlineitems():
andy@34 354 return map(_autotype, getnextvalidline().split())
andy@34 355
andy@34 356 ## Populate MASS block and create decaying particle objects
andy@35 357 masses = Block("MASS")
andy@33 358 numentries = int(getnextvalidline())
andy@33 359 for i in xrange(numentries):
andy@34 360 hwid, mass, lifetime = getnextvalidlineitems()
andy@34 361 width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
andy@34 362 pdgid = HERWIGID2PDGID.get(hwid, hwid)
andy@34 363 masses.add_entry((pdgid, mass))
andy@34 364 decays[pdgid] = Particle(pdgid, width, mass)
andy@34 365 #print pdgid, mass, width
andy@34 366 blocks["MASS"] = masses
andy@33 367
andy@34 368 ## Populate decays
andy@34 369 for n in xrange(numentries):
andy@34 370 numdecays = int(getnextvalidline())
andy@34 371 for d in xrange(numdecays):
andy@40 372 #print n, numentries-1, d, numdecays-1
andy@34 373 decayitems = getnextvalidlineitems()
andy@34 374 hwid = decayitems[0]
andy@34 375 pdgid = HERWIGID2PDGID.get(hwid, hwid)
andy@34 376 br = decayitems[1]
andy@34 377 nme = decayitems[2]
andy@34 378 daughter_hwids = decayitems[3:]
andy@34 379 daughter_pdgids = []
andy@34 380 for hw in daughter_hwids:
andy@34 381 if hw != 0:
andy@34 382 daughter_pdgids.append(HERWIGID2PDGID.get(hw, hw))
andy@35 383 if not decays.has_key(pdgid):
andy@40 384 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
andy@38 385 decays[pdgid] = Particle(pdgid)
andy@38 386 decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
andy@33 387
andy@33 388
andy@34 389 ## Now the SUSY parameters
andy@34 390 TANB, ALPHAH = getnextvalidlineitems()
andy@34 391 blocks["MINPAR"] = Block("MINPAR")
andy@34 392 blocks["MINPAR"].add_entry((3, TANB))
andy@34 393 blocks["ALPHA"] = Block("ALPHA")
andy@34 394 blocks["ALPHA"].entries = ALPHAH
andy@34 395 #
andy@34 396 ## Neutralino mixing matrix
andy@34 397 blocks["NMIX"] = Block("NMIX")
andy@34 398 for i in xrange(1, 5):
andy@34 399 nmix_i = getnextvalidlineitems()
andy@34 400 for j, v in enumerate(nmix_i):
andy@34 401 blocks["NMIX"].add_entry((i, j+1, v))
andy@34 402 #
andy@34 403 ## Chargino mixing matrices V and U
andy@34 404 blocks["VMIX"] = Block("VMIX")
andy@34 405 vmix = getnextvalidlineitems()
andy@34 406 blocks["VMIX"].add_entry((1, 1, vmix[0]))
andy@34 407 blocks["VMIX"].add_entry((1, 2, vmix[1]))
andy@34 408 blocks["VMIX"].add_entry((2, 1, vmix[2]))
andy@34 409 blocks["VMIX"].add_entry((2, 2, vmix[3]))
andy@34 410 blocks["UMIX"] = Block("UMIX")
andy@34 411 umix = getnextvalidlineitems()
andy@34 412 blocks["UMIX"].add_entry((1, 1, umix[0]))
andy@34 413 blocks["UMIX"].add_entry((1, 2, umix[1]))
andy@34 414 blocks["UMIX"].add_entry((2, 1, umix[2]))
andy@34 415 blocks["UMIX"].add_entry((2, 2, umix[3]))
andy@34 416 #
andy@34 417 THETAT, THETAB, THETAL = getnextvalidlineitems()
andy@34 418 import math
andy@34 419 blocks["STOPMIX"] = Block("STOPMIX")
andy@34 420 blocks["STOPMIX"].add_entry((1, 1, math.cos(THETAT)))
andy@34 421 blocks["STOPMIX"].add_entry((1, 2, -math.sin(THETAT)))
andy@34 422 blocks["STOPMIX"].add_entry((2, 1, math.sin(THETAT)))
andy@34 423 blocks["STOPMIX"].add_entry((2, 2, math.cos(THETAT)))
andy@34 424 blocks["SBOTMIX"] = Block("SBOTMIX")
andy@34 425 blocks["SBOTMIX"].add_entry((1, 1, math.cos(THETAB)))
andy@34 426 blocks["SBOTMIX"].add_entry((1, 2, -math.sin(THETAB)))
andy@34 427 blocks["SBOTMIX"].add_entry((2, 1, math.sin(THETAB)))
andy@34 428 blocks["SBOTMIX"].add_entry((2, 2, math.cos(THETAB)))
andy@34 429 blocks["STAUMIX"] = Block("STAUMIX")
andy@34 430 blocks["STAUMIX"].add_entry((1, 1, math.cos(THETAL)))
andy@34 431 blocks["STAUMIX"].add_entry((1, 2, -math.sin(THETAL)))
andy@34 432 blocks["STAUMIX"].add_entry((2, 1, math.sin(THETAL)))
andy@34 433 blocks["STAUMIX"].add_entry((2, 2, math.cos(THETAL)))
andy@34 434 #
andy@34 435 ATSS, ABSS, ALSS = getnextvalidlineitems()
andy@34 436 blocks["AU"] = Block("AU")
andy@34 437 blocks["AU"].add_entry((3, 3, ATSS))
andy@34 438 blocks["AD"] = Block("AD")
andy@34 439 blocks["AD"].add_entry((3, 3, ABSS))
andy@34 440 blocks["AE"] = Block("AE")
andy@34 441 blocks["AE"].add_entry((3, 3, ALSS))
andy@34 442 #
andy@34 443 MUSS = getnextvalidlineitems()[0]
andy@34 444 blocks["MINPAR"].add_entry((4, MUSS))
andy@34 445 #
andy@34 446 return blocks, decays
andy@33 447
andy@21 448
andy@31 449 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
andy@29 450 """
andy@29 451 Write an SLHA file from the supplied blocks and decays dicts.
andy@31 452
andy@31 453 Other keyword parameters are passed to writeSLHA.
andy@29 454 """
andy@29 455 with open(spcfilename, "w") as f:
andy@31 456 f.write(writeSLHA(blocks, decays, kwargs))
andy@29 457
andy@29 458
andy@33 459 ## TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
andy@33 460
andy@31 461 def writeSLHA(blocks, decays, ignorenobr=False):
andy@29 462 """
andy@29 463 Return an SLHA definition as a string, from the supplied blocks and decays dicts.
andy@29 464 """
andy@31 465 sep = " "
andy@29 466 out = ""
andy@30 467 def dict_hier_strs(d, s=""):
andy@30 468 if type(d) is dict:
andy@30 469 for k, v in sorted(d.iteritems()):
andy@31 470 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
andy@30 471 yield s2
andy@30 472 else:
andy@31 473 yield s + sep + _autostr(d)
andy@30 474 ## Blocks
andy@30 475 for bname, b in sorted(blocks.iteritems()):
andy@29 476 namestr = b.name
andy@29 477 if b.q is not None:
andy@29 478 namestr += " Q= %e" % b.q
andy@29 479 out += "BLOCK %s\n" % namestr
andy@30 480 for s in dict_hier_strs(b.entries):
andy@31 481 out += sep + s + "\n"
andy@29 482 out += "\n"
andy@30 483 ## Decays
andy@30 484 for pid, particle in sorted(decays.iteritems()):
andy@40 485 out += "DECAY %d %e\n" % (particle.pid, particle.totalwidth or -1)
andy@30 486 for d in sorted(particle.decays):
andy@31 487 if d.br > 0.0 or not ignorenobr:
andy@31 488 products_str = " ".join(map(str, d.ids))
andy@31 489 out += sep + "%e" % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
andy@29 490 out += "\n"
andy@29 491 return out
andy@29 492
andy@29 493
andy@29 494
andy@48 495 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
andy@48 496 """
andy@48 497 Write an ISAWIG file from the supplied blocks and decays dicts.
andy@48 498
andy@48 499 Other keyword parameters are passed to writeISAWIG.
andy@48 500
andy@48 501 TODO: Handle RPV SUSY
andy@48 502 """
andy@48 503 with open(isafilename, "w") as f:
andy@48 504 f.write(writeISAWIG(blocks, decays, kwargs))
andy@48 505
andy@48 506
andy@48 507 def writeISAWIG(blocks, decays, ignorenobr=False):
andy@48 508 """
andy@48 509 Return an ISAWIG definition as a string, from the supplied blocks and decays dicts.
andy@48 510
andy@48 511 ISAWIG parsing based on the HERWIG SUSY specification format, from
andy@48 512 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
andy@48 513
andy@48 514 If the ignorenobr parameter is True, do not write decay entries with a
andy@48 515 branching ratio of zero.
andy@48 516 """
andy@48 517 ## PDG MC ID codes mapped to HERWIG SUSY ID codes, based on
andy@48 518 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
andy@48 519 PDGID2HERWIGID = {}
andy@48 520 PDGID2HERWIGID[ 25] = 203 ## HIGGSL0 (ADDED)
andy@48 521 PDGID2HERWIGID[ 26] = 203 ## HIGGSL0
andy@48 522 PDGID2HERWIGID[ 35] = 204 ## HIGGSH0
andy@48 523 PDGID2HERWIGID[ 36] = 205 ## HIGGSA0
andy@48 524 PDGID2HERWIGID[ 37] = 206 ## HIGGS+
andy@48 525 PDGID2HERWIGID[ -37] = 207 ## HIGGS-
andy@48 526 PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
andy@48 527 PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
andy@48 528 PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
andy@48 529 PDGID2HERWIGID[-1000002] = 408 ## SSUL
andy@48 530 PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
andy@48 531 PDGID2HERWIGID[-1000003] = 409 ## SSSL
andy@48 532 PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
andy@48 533 PDGID2HERWIGID[-1000004] = 410 ## SSCL
andy@48 534 PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
andy@48 535 PDGID2HERWIGID[-1000005] = 411 ## SSB1
andy@48 536 PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
andy@48 537 PDGID2HERWIGID[-1000006] = 412 ## SST1
andy@48 538 PDGID2HERWIGID[ 2000001] = 413 ## SSDR
andy@48 539 PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
andy@48 540 PDGID2HERWIGID[ 2000002] = 414 ## SSUR
andy@48 541 PDGID2HERWIGID[-2000002] = 420 ## SSURBR
andy@48 542 PDGID2HERWIGID[ 2000003] = 415 ## SSSR
andy@48 543 PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
andy@48 544 PDGID2HERWIGID[ 2000004] = 416 ## SSCR
andy@48 545 PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
andy@48 546 PDGID2HERWIGID[ 2000005] = 417 ## SSB2
andy@48 547 PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
andy@48 548 PDGID2HERWIGID[ 2000006] = 418 ## SST2
andy@48 549 PDGID2HERWIGID[-2000006] = 424 ## SST2BR
andy@48 550 PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
andy@48 551 PDGID2HERWIGID[-1000011] = 431 ## SSEL+
andy@48 552 PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
andy@48 553 PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
andy@48 554 PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
andy@48 555 PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
andy@48 556 PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
andy@48 557 PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
andy@48 558 PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
andy@48 559 PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
andy@48 560 PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
andy@48 561 PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
andy@48 562 PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
andy@48 563 PDGID2HERWIGID[-2000011] = 443 ## SSEL+
andy@48 564 PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
andy@48 565 PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
andy@48 566 PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
andy@48 567 PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
andy@48 568 PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
andy@48 569 PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
andy@48 570 PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
andy@48 571 PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
andy@48 572 PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
andy@48 573 PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
andy@48 574 PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
andy@48 575 PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
andy@48 576 PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
andy@48 577 PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
andy@48 578 PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
andy@48 579 PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
andy@48 580 PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
andy@48 581 PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
andy@48 582 PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
andy@48 583 PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
andy@48 584
andy@48 585 masses = blocks["MASS"].entries
andy@48 586
andy@48 587 ## Init output string
andy@48 588 out = ""
andy@48 589
andy@48 590 ## First write out masses section:
andy@48 591 ## Number of SUSY + top particles
andy@48 592 ## IDHW, RMASS(IDHW), RLTIM(IDHW)
andy@48 593 ## repeated for each particle
andy@48 594 ## IDHW is the HERWIG identity code.
andy@48 595 ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
andy@48 596 massout = ""
andy@48 597 for pid in masses.keys():
andy@48 598 lifetime = -1
andy@48 599 try:
andy@48 600 width = decays[pid].totalwidth
andy@48 601 if width and width > 0:
andy@48 602 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
andy@48 603 except:
andy@48 604 pass
andy@48 605 massout += "%d %e %e\n" % (PDGID2HERWIGID.get(pid, pid), masses[pid], lifetime)
andy@48 606 out += "%d\n" % massout.count("\n")
andy@48 607 out += massout
andy@48 608
andy@48 609 assert(len(masses) == len(decays))
andy@48 610
andy@48 611 ## Next each particles decay modes together with their branching ratios and matrix element codes
andy@48 612 ## Number of decay modes for a given particle (IDK)
andy@48 613 ## IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
andy@48 614 ## repeated for each mode.
andy@48 615 ## Repeated for each particle.
andy@48 616 ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
andy@48 617 ## the decay mode. NME is a code for the matrix element to be used, either from the
andy@48 618 ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
andy@48 619 for i, pid in enumerate(decays.keys()):
andy@48 620 # if not decays.has_key(pid):
andy@48 621 # continue
andy@48 622 hwid = PDGID2HERWIGID.get(pid, pid)
andy@48 623 decayout = ""
andy@48 624 #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
andy@48 625 for i_d, d in enumerate(decays[pid].decays):
andy@48 626 ## Skip decay if it has no branching ratio
andy@48 627 if ignorenobr and d.br == 0:
andy@48 628 continue
andy@48 629 ## TODO: Identify decay matrix element to use
andy@59 630 ## From std HW docs, or from this pair:
andy@59 631 ## Two new matrix element codes have been added for these new decays:
andy@59 632 ## NME = 200 3 body top quark via charged Higgs
andy@59 633 ## 300 3 body R-parity violating gaugino and gluino decays
andy@59 634 ## How to determine which to use?
andy@48 635 #decayout += str(i_d + 1) + " "
andy@48 636 decayout += "%d %e 0 " % (hwid, d.br)
andy@66 637 ## Order decay products as required by HERWIG
andy@66 638 ## Top
andy@66 639 absids = map(abs, d.ids)
andy@66 640 if abs(pid) == 6:
andy@66 641 def cmp_bottomlast(a, b):
andy@66 642 """Comparison function which always puts b/bbar last"""
andy@66 643 if abs(a) == 5:
andy@66 644 return True
andy@66 645 elif abs(b) == 5:
andy@66 646 return False
andy@66 647 return cmp(a, b)
andy@66 648 if len(absids) == 2:
andy@66 649 ## 2 body mode, to Higgs: Higgs; Bottom
andy@66 650 if (25 in absids or 26 in absids) and 5 in absids:
andy@66 651 d.ids = sorted(d.ids, key=cmp_bottomlast)
andy@66 652 elif len(absids) == 3:
andy@66 653 ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
andy@66 654 if 37 in absids or 23 in absids:
andy@66 655 d.ids = sorted(d.ids, key=cmp_bottomlast)
andy@66 656 ## Gluino
andy@66 657 elif abs(pid) == 1000021:
andy@66 658 if len(absids) == 2:
andy@66 659 ## 2 body mode
andy@66 660 ## without gluon: any order
andy@66 661 ## with gluon: gluon; colour neutral
andy@66 662 if 21 in absids:
andy@66 663 def cmp_gluonfirst(a, b):
andy@66 664 """Comparison function which always puts gluon first"""
andy@66 665 if a == 21:
andy@66 666 return False
andy@66 667 elif b == 21:
andy@66 668 return True
andy@66 669 return cmp(a, b)
andy@66 670 d.ids = sorted(d.ids, key=cmp_gluonfirst)
andy@66 671 elif len(absids) == 3:
andy@66 672 ## 3-body modes, R-parity conserved: colour neutral; q or qbar
andy@66 673 pass
andy@66 674
andy@66 675 ## TODO: More sorting!
andy@67 676 ## Squark/Slepton
andy@67 677 ## 2 body modes: Gaugino/Gluino with Quark/Lepton Gaugino quark
andy@67 678 ## Gluino lepton
andy@67 679 ## 3 body modes: Weak sparticle particles from W decay
andy@67 680 ## Squark
andy@67 681 ## 2 body modes: Lepton Number Violated quark lepton
andy@67 682 ## Baryon number violated quark quark
andy@66 683
andy@67 684 ## TODO: Meaning of alignments below? Re-check web table
andy@67 685
andy@67 686 ## Slepton
andy@67 687 ## 2 body modes: Lepton Number Violated q or qbar
andy@67 688 ## Higgs
andy@67 689 ## 2 body modes: (s)quark-(s)qbar (s)q or (s)qbar
andy@67 690 ## (s)lepton-(s)lepton (s)l or (s)lbar
andy@67 691 ## 3 body modes: Colour neutral q or qbar
andy@67 692 ## l or lbar
andy@67 693 ## Gaugino
andy@67 694 ## 2 body modes: Squark-quark q or sq
andy@67 695 ## Slepton-lepton l or sl
andy@67 696 ##
andy@67 697 ## 3 body modes: R-parity conserved colour neutral q or qbar
andy@67 698 ## l or lbar
andy@67 699
andy@67 700 ## TODO: Gaugino/Gluino
andy@67 701 ## 3 body modes: R-parity violating: Particles in the same order as the R-parity violating superpotential
andy@66 702
andy@66 703 ## Pad out IDs list with zeros
andy@48 704 ids = [0,0,0,0,0]
andy@48 705 for i, pid in enumerate(d.ids):
andy@48 706 ids[i] = pid
andy@48 707 ids = map(str, ids)
andy@48 708 decayout += " ".join(ids) + "\n"
andy@48 709 decayout = "%d\n" % decayout.count("\n") + decayout
andy@48 710 out += decayout
andy@48 711
andy@48 712 ## Now the SUSY parameters
andy@48 713 ## TANB, ALPHAH:
andy@48 714 out += "%e %e\n" % (blocks["MINPAR"].entries[3], blocks["ALPHA"].entries)
andy@48 715 ## Neutralino mixing matrix
andy@48 716 nmix = blocks["NMIX"].entries
andy@48 717 for i in xrange(1, 5):
andy@48 718 out += "%e %e %e %e\n" % (nmix[i][1], nmix[i][2], nmix[i][3], nmix[i][4])
andy@48 719 ## Chargino mixing matrices V and U
andy@48 720 vmix = blocks["VMIX"].entries
andy@48 721 out += "%e %e %e %e\n" % (vmix[1][1], vmix[1][2], vmix[2][1], vmix[2][2])
andy@48 722 umix = blocks["UMIX"].entries
andy@48 723 out += "%e %e %e %e\n" % (umix[1][1], umix[1][2], umix[2][1], umix[2][2])
andy@48 724 # THETAT,THETAB,THETAL
andy@48 725 import math
andy@48 726 out += "%e %e %e\n" % (math.acos(blocks["STOPMIX"].entries[1][1]),
andy@48 727 math.acos(blocks["SBOTMIX"].entries[1][1]),
andy@48 728 math.acos(blocks["STAUMIX"].entries[1][1]))
andy@48 729 # ATSS,ABSS,ALSS
andy@48 730 out += "%e %e %e\n" % (blocks["AU"].entries[3][3],
andy@48 731 blocks["AD"].entries[3][3],
andy@48 732 blocks["AE"].entries[3][3])
andy@48 733 # MUSS == sign(mu)
andy@48 734 out += "%f\n" % blocks["MINPAR"].entries[4]
andy@48 735
andy@48 736 ## TODO: Handle RPV SUSY
andy@48 737
andy@48 738 return out
andy@48 739
andy@48 740
andy@48 741
andy@1 742 if __name__ == "__main__":
andy@1 743 import sys
andy@1 744 for a in sys.argv[1:]:
andy@35 745 if a.endswith(".isa"):
andy@35 746 blocks, decays = readISAWIGFile(a)
andy@35 747 else:
andy@35 748 blocks, decays = readSLHAFile(a)
andy@3 749
andy@5 750 for bname, b in sorted(blocks.iteritems()):
andy@5 751 print b
andy@5 752 print
andy@3 753
andy@47 754 print blocks.keys()
andy@47 755
andy@4 756 print blocks["MASS"].entries[25]
andy@6 757 print
andy@6 758
andy@6 759 for p in sorted(decays.values()):
andy@6 760 print p
andy@6 761 print
andy@29 762
andy@31 763 print writeSLHA(blocks, decays, ignorenobr=True)

mercurial