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