pyslha.py

Sat, 03 Mar 2012 00:28:48 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sat, 03 Mar 2012 00:28:48 +0100
changeset 162
e5b76d1917c0
parent 161
9306d1ad087f
child 164
fa5d5692a68b
permissions
-rw-r--r--

Converting rendering system to use the separate tex2pix library (which itself is based on the slhaplot rendering experience). Test for presence of the mathpazo package before including it in the default preamble. Version 1.3.0

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

mercurial