pyslha.py

Sat, 24 Sep 2011 21:17:07 +0100

author
Andy Buckley <andy@insectnation.org>
date
Sat, 24 Sep 2011 21:17:07 +0100
changeset 154
da6a12eb28f1
parent 152
8ca4c23ba4a6
child 158
5cb1da0bddac
permissions
-rw-r--r--

Adding a max mass cutoff for slhaplot: don't show particle lines or associated decays for particles with masses greater than the cutoff. Version 1.2.8

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

mercurial