pyslha.py

Mon, 08 Apr 2013 17:24:14 +0200

author
Andy Buckley <andy@insectnation.org>
date
Mon, 08 Apr 2013 17:24:14 +0200
changeset 185
6cae67cd3cc1
parent 183
5dcd3f1622be
child 188
e706e3b2a647
permissions
-rw-r--r--

pyslha.py: Fix handling of particles with NaN widths (and decay channels with NaN BRs). Both are filled into the blocks in readSLHA(File) now, with None as the value for the invalid width/BR: test for correctness with e.g. 'myparticle.width is not None'. BR == NaN decays will be added if the ignorenobr arg == False (the default). Thanks to Lukas Vanelderen for the report.

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

mercurial