pyslha.py

Fri, 26 Apr 2013 20:57:42 +0200

author
Andy Buckley <andy@insectnation.org>
date
Fri, 26 Apr 2013 20:57:42 +0200
changeset 192
453a523cba25
parent 189
60c73b489420
child 193
40d024dac179
permissions
-rw-r--r--

Converting entries storage to avoid recursion, by using tuples as dict keys if necessary. Adding direct entry access via the [] operator to the Block class's entries. Improved SLHA output formatting.

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

mercurial