pyslha.py@269bc4611cdd
pyslha.py
Mon, 15 Jul 2013 16:50:38 +0200
- author
- Andy Buckley <andy@insectnation.org>
- date
- Mon, 15 Jul 2013 16:50:38 +0200
- changeset 244
- 269bc4611cdd
- parent 242
- 41ec8fb352b3
- child 245
- ca52fe2c0d71
- permissions
- -rw-r--r--
Add credit to ChangeLog
1 #! /usr/bin/env python
3 """\
4 A simple but flexible handler of the SUSY Les Houches Accord (SLHA) data format.
6 pyslha is a parser/writer module for particle physics SUSY Les Houches Accord
7 (SLHA) supersymmetric spectrum/decay files, and a collection of scripts which
8 use the interface, e.g. for conversion to and from the legacy ISAWIG format, or
9 to plot the mass spectrum and decay chains.
11 The current release supports SLHA version 1, and as far as I'm aware is also
12 fully compatible with SLHA2: the block structures are read and accessed
13 generically. If you have any problems, please provide an example input file and
14 I'll happily investigate.
16 The plotting script provides output in PDF, EPS and PNG via LaTeX and the TikZ
17 graphics package, and as LaTeX/TikZ source for direct embedding into documents or
18 user-tweaking of the generated output.
20 Users of version 1.x should note that the interface has changed a little in
21 version 2.0.0 and onward, in particular in the interface of the Block objects,
22 which are now more dict-like: entries can be added and accessed via the usual
23 square-brackets indexing operators, including for multiple indices as is common
24 for mixing matrices e.g. NMIX[1,2] as opposed to the old NMIX.entries[1][2]
25 way. This does break backward compatibility but is a big improvement both for
26 internal code sanity and usability. The Block interface also now supplies
27 dict-like has_key(), keys(), and items() methods, as well as more specialist
28 value(), set_value() and is_single_valued() methods for improved access to ALPHA
29 and any other unindexed blocks.
31 If you use PySLHA, for either model data handling or spectrum visualisation,
32 please cite the paper: http://arxiv.org/abs/1305.4194
34 TODOs:
36 For 2.1.4:
37 * In set_value, if first item is non-int, treat as None-indexed
38 * Refine value string heuristic for strings with ints in them?
40 For 2.2.0:
41 * Add Sphinx docs.
42 * Preserve comments from read -> write (needs full-line/inline comment
43 separation). Can use separate comment dicts in Block and Decay, and
44 attach a multiline .comment attr to the returned/written dicts.
46 For 3.0.0:
47 * Add a new Doc type which encapsulates the blocks, decays, xsections,
48 etc. data sections and will act as a more robust single return value from
49 the read commands / argument to the write commands. Can also handle
50 document-level header comments this way and improve the __repr__ resulting
51 from an interactive read() call without a target variable (i.e. hide the
52 uninteresting OrderedDict boilerplate).
54 Later/maybe:
55 * Identify HERWIG decay matrix element to use in ISAWIG
56 * Handle RPV SUSY in ISAWIG
57 """
59 __author__ = "Andy Buckley <andy.buckley@cern.ch>"
60 __version__ = "2.1.3"
64 ###############################################################################
65 ## Private utility functions
67 def _mkdict():
68 """Try to return an empty ordered dict, but fall back to normal dict if necessary"""
69 try:
70 from collections import OrderedDict
71 return OrderedDict()
72 except:
73 try:
74 from ordereddict import OrderedDict
75 return OrderedDict()
76 except:
77 return dict()
79 def _autotype(var):
80 """Automatically convert strings to numerical types if possible."""
81 if type(var) is not str:
82 return var
83 if var.isdigit() or (var.startswith("-") and var[1:].isdigit()):
84 return int(var)
85 try:
86 f = float(var)
87 return f
88 except ValueError:
89 return var
91 def _autostr(var, precision=8):
92 """Automatically format numerical types as the right sort of string."""
93 if type(var) is float:
94 return ("% ." + str(precision) + "e") % var
95 elif not hasattr(var, "__iter__"):
96 return str(var)
97 else:
98 return ",".join(_autostr(subval) for subval in var)
100 def _autotuple(a):
101 """Automatically convert the supplied iterable to a scalar or tuple as appropriate."""
102 if len(a) == 1:
103 return a[0]
104 return tuple(a)
108 ###############################################################################
109 ## Exceptions
111 class AccessError(Exception):
112 "Exception object to be raised when a SLHA block is accessed in an invalid way"
113 def __init__(self, errmsg):
114 self.msg = errmsg
115 def __str__(self):
116 return self.msg
118 class ParseError(Exception):
119 "Exception object to be raised when a spectrum file/string is malformed"
120 def __init__(self, errmsg):
121 self.msg = errmsg
122 def __str__(self):
123 return self.msg
127 ###############################################################################
128 ## The data block, decay and particle classes
130 class Block(object):
131 """
132 Object representation of any BLOCK elements read from an SLHA file.
134 Blocks have a name, may have an associated Q value, and contain a collection
135 of data entries, each indexed by one or more keys. Entries in the dictionary
136 are stored as numeric types (int or float) when a cast from the string in
137 the file has been possible.
139 Block is closely related to a Python dict (and, in fact, is implemented via
140 an OrderedDict when possible). The preferred methods of entry access use the
141 dict-like [] operator for getting and setting, and the keys() and items()
142 methods for iteration. Purely iterating over the object behaves like keys(),
143 as for an ordinary dict.
145 Multiple (integer) indices are possible, especially for entries in mixing
146 matrix blocks. These are now implemented in the natural way, e.g. for access
147 to the (1,2) element of a mixing matrix block, use bmix[1,2] = 0.123 and
148 print bmix[1,2]. The value() and set_value() functions behave
149 similarly. Multi-element values are also permitted.
151 It is possible, although not usual, to store unindexed values in a
152 block. This is only supported when that entry is the only one in the block,
153 and it is stored in the normal fashion but with None as the lookup key. The
154 value() method may be used without a key argument to retrieve this value, if
155 the is_single_valued() method returns True, and similarly the set_value()
156 method may be used to set it if only one argument is supplied and the object
157 is compatible.
158 """
159 def __init__(self, name, q=None, entries=None):
160 self.name = name
161 self.entries = _mkdict()
162 if entries is not None:
163 self.entries.update(entries)
164 self.q = _autotype(q)
166 # TODO: Rename? To what?
167 def add_entry(self, args):
168 """Add an entry to the block from an iterable (i.e. list or tuple) of
169 strings, or from a whitespace-separated string.
171 This method is just for convenience: it splits the single string
172 argument if necessary and converts the list of strings into numeric
173 types when possible. For the treatment of the resulting iterable see the
174 set_value method.
175 """
176 ## If the argument is a single string, split it and proceed
177 if type(args) is str:
178 args = args.split()
179 ## Check that the arg is an iterable
180 if not hasattr(args, "__iter__"):
181 raise AccessError("Block entries must be iterable")
182 ## Auto-convert the types in the list
183 args = map(_autotype, args)
184 ## Re-join consecutive strings into single entries
185 i = 0
186 while i < len(args)-1:
187 if type(args[i]) is str and type(args[i+1]) is str:
188 args[i] += " " + args[i+1]
189 del args[i+1]
190 continue
191 i += 1
192 ## Add the entry to the map, with appropriate indices
193 self.set_value(*args)
195 def is_single_valued(self):
196 """Return true if there is only one entry, and it has no index: the
197 'value()' attribute may be used in that case without an argument."""
198 return len(self.entries) == 1 and self.entries.keys()[0] is None
200 def value(self, key=None):
201 """Get a value from the block with the supplied key.
203 If no key is given, then the block must contain only one non-indexed
204 value otherwise an AccessError exception will be raised.\
205 """
206 if key == None and not self.is_single_valued():
207 raise AccessError("Tried to access unique value of multi-value block")
208 return self.entries[key]
210 def set_value(self, *args):
211 """Set a value in the block via supplied key/value arguments.
213 Indexing is determined automatically: any leading integers will be
214 treated as a multi-dimensional index, with the remaining entries being a
215 (potentially multi-dimensional) value. If all N args are ints, then the
216 first N-1 are treated as the index and the Nth as the value.
218 If there is only one arg it will be treated as the value of a
219 single-valued block. In this case the block must already contain at most
220 one non-indexed value otherwise an AccessError exception will be
221 raised.\
222 """
223 if len(args) == 0:
224 raise AccessError("set_value() called without arguments")
225 elif len(args) == 1:
226 if len(self.entries) > 0 and not self.is_single_valued():
227 raise AccessError("Tried to set a unique value on a multi-value block")
228 self.entries[None] = args[0]
229 else:
230 ## Find the first non-integer -- all previous items are indices
231 i_first_nonint = -1
232 for i, x in enumerate(args):
233 if type(x) is not int:
234 i_first_nonint = i
235 break
236 if i_first_nonint == 0:
237 raise AccessError("Attempted to set a block entry with a non-integer(s) index")
238 else:
239 self.entries[_autotuple(args[:i_first_nonint])] = _autotuple(args[i_first_nonint:])
241 def has_key(self, key):
242 """Does the block have the given key?"""
243 return self.entries.has_key(key)
245 def keys(self):
246 """Access the block item keys."""
247 return self.entries.keys()
249 def values(self):
250 """Access the block item values."""
251 return self.entries.values()
253 def items(self, key=None):
254 """Access the block items as (key,value) tuples.
256 Note: The Python 3 dict attribute 'items()' is used rather than the
257 'old' Python 2 'iteritems()' name for forward-looking compatibility.\
258 """
259 return self.entries.iteritems()
261 def __len__(self):
262 return len(self.entries)
264 def __iter(self):
265 return self.entries.__iter__()
267 def __getitem__(self, key):
268 return self.entries[key]
270 def __setitem__(self, key, value):
271 if key is not None and type(key) is not int and not all(type(x) is int for x in key):
272 raise AccessError("Attempted to set a block entry with a non-integer(s) index")
273 self.entries[key] = value
275 def __cmp__(self, other):
276 return cmp(self.name, other.name) and cmp(self.entries, other.entries)
278 def __repr__(self):
279 s = self.name
280 if self.q is not None:
281 s += " (Q=%s)" % self.q
282 s += " { " + "; ".join(_autostr(k) + " : " + _autostr(v) for k, v in self.items()) + " }"
283 return s
286 class Decay(object):
287 """
288 Object representing a decay entry on a particle decribed by the SLHA file.
289 'Decay' objects are not a direct representation of a DECAY block in an SLHA
290 file... that role, somewhat confusingly, is taken by the Particle class.
292 Decay objects have three properties: a branching ratio, br, an nda number
293 (number of daughters == len(ids)), and a tuple of PDG PIDs to which the
294 decay occurs. The PDG ID of the particle whose decay this represents may
295 also be stored, but this is normally known via the Particle in which the
296 decay is stored.
297 """
298 def __init__(self, br, nda, ids, parentid=None):
299 self.parentid = parentid
300 self.br = br
301 self.nda = nda
302 self.ids = ids
303 assert(self.nda == len(self.ids))
305 def __cmp__(self, other):
306 return cmp(other.br, self.br)
308 def __repr__(self):
309 return "% .8e %s" % (self.br, self.ids)
312 class Particle(object):
313 """
314 Representation of a single, specific particle, decay block from an SLHA
315 file. These objects are not themselves called 'Decay', since that concept
316 applies more naturally to the various decays found inside this
317 object. Particle classes store the PDG ID (pid) of the particle being
318 represented, and optionally the mass (mass) and total decay width
319 (totalwidth) of that particle in the SLHA scenario. Masses may also be found
320 via the MASS block, from which the Particle.mass property is filled, if at
321 all. They also store a list of Decay objects (decays) which are probably the
322 item of most interest.
323 """
324 def __init__(self, pid, totalwidth=None, mass=None):
325 self.pid = pid
326 self.totalwidth = totalwidth
327 self.mass = mass
328 self.decays = []
330 def add_decay(self, br, nda, ids):
331 self.decays.append(Decay(br, nda, ids))
332 self.decays.sort()
334 def __cmp__(self, other):
335 if abs(self.pid) == abs(other.pid):
336 return cmp(self.pid, other.pid)
337 return cmp(abs(self.pid), abs(other.pid))
339 def __repr__(self):
340 s = str(self.pid)
341 if self.mass is not None:
342 s += " : mass = %.8e GeV" % self.mass
343 if self.totalwidth is not None:
344 s += " : total width = %.8e GeV" % self.totalwidth
345 for d in self.decays:
346 if d.br > 0.0:
347 s += "\n %s" % d
348 return s
352 ###############################################################################
353 ## SLHA parsing and writing functions
355 def readSLHA(spcstr, ignorenobr=False, ignorenomass=False):
356 """
357 Read an SLHA definition from a string, returning dictionaries of blocks and
358 decays.
360 If the ignorenobr parameter is True, do not store decay entries with a
361 branching ratio of zero.
363 If the ignorenomass parameter is True, parse file even if mass block is
364 absent in the file (default is to raise a ParseError).
365 """
366 blocks = _mkdict()
367 decays = _mkdict()
368 #
369 import re
370 currentblock = None
371 currentdecay = None
372 for line in spcstr.splitlines():
373 ## Handle (ignore) comment lines
374 # TODO: Store block/entry comments
375 if line.startswith("#"):
376 continue
377 if "#" in line:
378 line = line[:line.index("#")]
379 ## Ignore empty lines (after comment removal and whitespace trimming)
380 if not line.strip():
381 continue
383 ## Section header lines start with a non-whitespace character, data lines have a whitespace indent
384 # TODO: Are tabs also allowed for indents? Check the SLHA standard.
385 if not line.startswith(" "):
386 # TODO: Should we now strip the line to remove any trailing whitespace?
388 ## Handle BLOCK start lines
389 if line.upper().startswith("BLOCK"):
390 #print line
391 match = re.match(r"BLOCK\s+(\w+)(\s+Q\s*=\s*.+)?", line.upper())
392 if not match:
393 continue
394 blockname = match.group(1)
395 qstr = match.group(2)
396 if qstr is not None:
397 qstr = qstr[qstr.find("=")+1:].strip()
398 currentblock = blockname
399 currentdecay = None
400 blocks[blockname] = Block(blockname, q=qstr)
401 ## Handle DECAY start lines
402 elif line.upper().startswith("DECAY"):
403 match = re.match(r"DECAY\s+(-?\d+)\s+([\d\.E+-]+|NAN).*", line.upper())
404 if not match:
405 continue
406 pdgid = int(match.group(1))
407 width = float(match.group(2)) if match.group(2) != "NAN" else None
408 currentblock = "DECAY"
409 currentdecay = pdgid
410 decays[pdgid] = Particle(pdgid, width)
411 ## Handle unknown section type start lines (and continue ignoring until a non-header line is found)
412 elif type(_autotype(line.split()[0])) is str:
413 import sys
414 sys.stderr.write("Ignoring unknown section type: %s\n" % line.split()[0])
415 currentblock = None
416 currentdecay = None
418 ## This non-empty line starts with an indent, hence must be an in-section data line
419 else:
420 # TODO: Should we now strip the line to remove the indent (and any trailing whitespace)?
421 if currentblock is not None:
422 items = line.split()
423 if len(items) < 1:
424 continue
425 if currentblock != "DECAY":
426 blocks[currentblock].add_entry(items)
427 # TODO: Add handling of XSECTION if/when standardised
428 else:
429 br = float(items[0]) if items[0].upper() != "NAN" else None
430 nda = int(items[1])
431 ids = map(int, items[2:])
432 if br > 0.0 or not ignorenobr: # br == None is < 0
433 decays[currentdecay].add_decay(br, nda, ids)
435 ## Try to populate Particle masses from the MASS block
436 # print blocks.keys()
437 try:
438 for pid in blocks["MASS"].keys():
439 if decays.has_key(pid):
440 decays[pid].mass = blocks["MASS"][pid]
441 except:
442 if not ignorenomass:
443 raise ParseError("No MASS block found: cannot populate particle masses")
445 return blocks, decays
449 def writeSLHABlocks(blocks, precision=8):
450 """Return an SLHA definition as a string, from the supplied blocks dict."""
451 sep = 3 * " "
452 blockstrs = []
453 for bname, b in blocks.iteritems():
454 namestr = b.name
455 if b.q is not None:
456 namestr += " Q= " + _autostr(float(b.q), precision)
457 blockstr = "BLOCK %s\n" % namestr
458 entrystrs = []
459 for k, v in b.items():
460 entrystr = ""
461 if type(k) == tuple:
462 entrystr += sep.join(_autostr(i, precision) for i in k)
463 elif k is not None:
464 entrystr += _autostr(k, precision)
465 entrystr += sep + _autostr(v, precision)
466 entrystrs.append(entrystr)
467 blockstr += "\n".join(entrystrs)
468 blockstrs.append(blockstr)
469 return "\n\n".join(blockstrs)
473 def writeSLHADecays(decays, ignorenobr=False, precision=8):
474 """Return an SLHA decay definition as a string, from the supplied decays dict."""
475 sep = 3 * " "
476 blockstrs = []
477 for pid, particle in decays.iteritems():
478 blockstr = ("DECAY %d " % particle.pid) + _autostr(particle.totalwidth or -1, precision) + "\n"
479 decaystrs = []
480 for d in particle.decays:
481 if d.br > 0.0 or not ignorenobr:
482 products_str = sep.join("% d" % i for i in d.ids)
483 decaystr = sep + _autostr(d.br, precision) + sep + str(len(d.ids)) + sep + products_str
484 decaystrs.append(decaystr)
485 blockstr += "\n".join(decaystrs)
486 blockstrs.append(blockstr)
487 return "\n\n".join(blockstrs)
491 def writeSLHA(blocks, decays, ignorenobr=False, precision=8):
492 """Return an SLHA definition as a string, from the supplied blocks and decays dicts."""
493 ss = [x for x in (writeSLHABlocks(blocks, precision), writeSLHADecays(decays, ignorenobr, precision)) if x]
494 return "\n\n".join(ss)
498 ###############################################################################
499 ## PDG <-> HERWIG particle ID code translations for ISAWIG handling
501 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
502 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
503 ## + the IDPDG array and section 4.13 of the HERWIG manual.
504 _HERWIGID2PDGID = {}
505 _HERWIGID2PDGID[7] = -1
506 _HERWIGID2PDGID[8] = -2
507 _HERWIGID2PDGID[9] = -3
508 _HERWIGID2PDGID[10] = -4
509 _HERWIGID2PDGID[11] = -5
510 _HERWIGID2PDGID[12] = -6
511 _HERWIGID2PDGID[13] = 21
512 _HERWIGID2PDGID[59] = 22
513 _HERWIGID2PDGID[121] = 11
514 _HERWIGID2PDGID[122] = 12
515 _HERWIGID2PDGID[123] = 13
516 _HERWIGID2PDGID[124] = 14
517 _HERWIGID2PDGID[125] = 15
518 _HERWIGID2PDGID[126] = 16
519 _HERWIGID2PDGID[127] = -11
520 _HERWIGID2PDGID[128] = -12
521 _HERWIGID2PDGID[129] = -13
522 _HERWIGID2PDGID[130] = -14
523 _HERWIGID2PDGID[131] = -15
524 _HERWIGID2PDGID[132] = -16
525 _HERWIGID2PDGID[198] = 24 # W+
526 _HERWIGID2PDGID[199] = -24 # W-
527 _HERWIGID2PDGID[200] = 23 # Z0
528 _HERWIGID2PDGID[201] = 25 ## SM HIGGS
529 _HERWIGID2PDGID[203] = 25 ## HIGGSL0 (== PDG standard in this direction)
530 _HERWIGID2PDGID[204] = 35 ## HIGGSH0
531 _HERWIGID2PDGID[205] = 36 ## HIGGSA0
532 _HERWIGID2PDGID[206] = 37 ## HIGGS+
533 _HERWIGID2PDGID[207] = -37 ## HIGGS-
534 _HERWIGID2PDGID[401] = 1000001 ## SSDLBR
535 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
536 _HERWIGID2PDGID[402] = 1000002 ## SSULBR
537 _HERWIGID2PDGID[408] = -1000002 ## SSUL
538 _HERWIGID2PDGID[403] = 1000003 ## SSSLBR
539 _HERWIGID2PDGID[409] = -1000003 ## SSSL
540 _HERWIGID2PDGID[404] = 1000004 ## SSCLBR
541 _HERWIGID2PDGID[410] = -1000004 ## SSCL
542 _HERWIGID2PDGID[405] = 1000005 ## SSB1BR
543 _HERWIGID2PDGID[411] = -1000005 ## SSB1
544 _HERWIGID2PDGID[406] = 1000006 ## SST1BR
545 _HERWIGID2PDGID[412] = -1000006 ## SST1
546 _HERWIGID2PDGID[413] = 2000001 ## SSDR
547 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
548 _HERWIGID2PDGID[414] = 2000002 ## SSUR
549 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
550 _HERWIGID2PDGID[415] = 2000003 ## SSSR
551 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
552 _HERWIGID2PDGID[416] = 2000004 ## SSCR
553 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
554 _HERWIGID2PDGID[417] = 2000005 ## SSB2
555 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
556 _HERWIGID2PDGID[418] = 2000006 ## SST2
557 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
558 _HERWIGID2PDGID[425] = 1000011 ## SSEL-
559 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
560 _HERWIGID2PDGID[426] = 1000012 ## SSNUEL
561 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
562 _HERWIGID2PDGID[427] = 1000013 ## SSMUL-
563 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
564 _HERWIGID2PDGID[428] = 1000014 ## SSNUMUL
565 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
566 _HERWIGID2PDGID[429] = 1000015 ## SSTAU1-
567 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
568 _HERWIGID2PDGID[430] = 1000016 ## SSNUTL
569 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
570 _HERWIGID2PDGID[437] = 2000011 ## SSEL-
571 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
572 _HERWIGID2PDGID[438] = 2000012 ## SSNUEL
573 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
574 _HERWIGID2PDGID[439] = 2000013 ## SSMUL-
575 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
576 _HERWIGID2PDGID[440] = 2000014 ## SSNUMUL
577 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
578 _HERWIGID2PDGID[441] = 2000015 ## SSTAU1-
579 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
580 _HERWIGID2PDGID[442] = 2000016 ## SSNUTL
581 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
582 _HERWIGID2PDGID[449] = 1000021 ## GLUINO
583 _HERWIGID2PDGID[450] = 1000022 ## NTLINO1
584 _HERWIGID2PDGID[451] = 1000023 ## NTLINO2
585 _HERWIGID2PDGID[452] = 1000025 ## NTLINO3
586 _HERWIGID2PDGID[453] = 1000035 ## NTLINO4
587 _HERWIGID2PDGID[454] = 1000024 ## CHGINO1+
588 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
589 _HERWIGID2PDGID[455] = 1000037 ## CHGINO2+
590 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
591 _HERWIGID2PDGID[458] = 1000039 ## GRAVTINO
593 def herwigid2pdgid(hwid):
594 """
595 Convert a particle ID code in the HERWIG internal IDHW format (as used by
596 ISAWIG) into its equivalent in the standard PDG ID code definition.
597 """
598 return _HERWIGID2PDGID.get(hwid, hwid)
601 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
602 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
603 ## + the IDPDG array and section 4.13 of the HERWIG manual.
604 _PDGID2HERWIGID = {}
605 _PDGID2HERWIGID[ -1] = 7
606 _PDGID2HERWIGID[ -2] = 8
607 _PDGID2HERWIGID[ -3] = 9
608 _PDGID2HERWIGID[ -4] = 10
609 _PDGID2HERWIGID[ -5] = 11
610 _PDGID2HERWIGID[ -6] = 12
611 _PDGID2HERWIGID[ 21] = 13
612 _PDGID2HERWIGID[ 22] = 59
613 _PDGID2HERWIGID[ 11] = 121
614 _PDGID2HERWIGID[ 12] = 122
615 _PDGID2HERWIGID[ 13] = 123
616 _PDGID2HERWIGID[ 14] = 124
617 _PDGID2HERWIGID[ 15] = 125
618 _PDGID2HERWIGID[ 16] = 126
619 _PDGID2HERWIGID[ -11] = 127
620 _PDGID2HERWIGID[ -12] = 128
621 _PDGID2HERWIGID[ -13] = 129
622 _PDGID2HERWIGID[ -14] = 130
623 _PDGID2HERWIGID[ -15] = 131
624 _PDGID2HERWIGID[ -16] = 132
625 _PDGID2HERWIGID[ 24] = 198 ## W+
626 _PDGID2HERWIGID[ -24] = 199 ## W-
627 _PDGID2HERWIGID[ 23] = 200 ## Z
628 _PDGID2HERWIGID[ 25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW) # TODO: should be 201?
629 _PDGID2HERWIGID[ 26] = 203 ## HIGGSL0
630 _PDGID2HERWIGID[ 35] = 204 ## HIGGSH0
631 _PDGID2HERWIGID[ 36] = 205 ## HIGGSA0
632 _PDGID2HERWIGID[ 37] = 206 ## HIGGS+
633 _PDGID2HERWIGID[ -37] = 207 ## HIGGS-
634 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
635 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
636 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
637 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
638 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
639 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
640 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
641 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
642 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
643 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
644 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
645 _PDGID2HERWIGID[-1000006] = 412 ## SST1
646 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
647 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
648 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
649 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
650 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
651 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
652 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
653 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
654 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
655 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
656 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
657 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
658 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
659 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
660 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
661 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
662 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
663 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
664 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
665 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
666 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
667 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
668 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
669 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
670 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
671 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
672 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
673 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
674 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
675 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
676 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
677 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
678 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
679 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
680 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
681 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
682 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
683 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
684 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
685 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
686 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
687 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
688 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
689 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
690 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
691 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
693 def pdgid2herwigid(pdgid):
694 """
695 Convert a particle ID code in the standard PDG ID code definition into
696 its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
697 """
698 return _PDGID2HERWIGID.get(pdgid, pdgid)
701 ###############################################################################
702 ## ISAWIG format reading/writing
704 def readISAWIG(isastr, ignorenobr=False):
705 """
706 Read a spectrum definition from a string in the ISAWIG format, returning
707 dictionaries of blocks and decays. While this is not an SLHA format, it is
708 informally supported as a useful mechanism for converting ISAWIG spectra to
709 SLHA.
711 ISAWIG parsing based on the HERWIG SUSY specification format, from
712 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
714 If the ignorenobr parameter is True, do not store decay entries with a
715 branching ratio of zero.
716 """
718 blocks = _mkdict()
719 decays = _mkdict()
720 LINES = isastr.splitlines()
722 def getnextvalidline():
723 while LINES:
724 s = LINES.pop(0).strip()
725 # print "*", s, "*"
726 ## Return None if EOF reached
727 if len(s) == 0:
728 continue
729 ## Strip comments
730 if "#" in s:
731 s = s[:s.index("#")].strip()
732 ## Return if non-empty
733 if len(s) > 0:
734 return s
736 def getnextvalidlineitems():
737 return map(_autotype, getnextvalidline().split())
739 ## Populate MASS block and create decaying particle objects
740 masses = Block("MASS")
741 numentries = int(getnextvalidline())
742 for i in xrange(numentries):
743 hwid, mass, lifetime = getnextvalidlineitems()
744 width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
745 pdgid = herwigid2pdgid(hwid)
746 masses[pdgid] = mass
747 decays[pdgid] = Particle(pdgid, width, mass)
748 #print pdgid, mass, width
749 blocks["MASS"] = masses
751 ## Populate decays
752 for n in xrange(numentries):
753 numdecays = int(getnextvalidline())
754 for d in xrange(numdecays):
755 #print n, numentries-1, d, numdecays-1
756 decayitems = getnextvalidlineitems()
757 hwid = decayitems[0]
758 pdgid = herwigid2pdgid(hwid)
759 br = decayitems[1]
760 nme = decayitems[2]
761 daughter_hwids = decayitems[3:]
762 daughter_pdgids = []
763 for hw in daughter_hwids:
764 if hw != 0:
765 daughter_pdgids.append(herwigid2pdgid(hw))
766 if not decays.has_key(pdgid):
767 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
768 decays[pdgid] = Particle(pdgid)
769 decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
772 ## Now the SUSY parameters
773 TANB, ALPHAH = getnextvalidlineitems()
774 blocks["MINPAR"] = Block("MINPAR")
775 blocks["MINPAR"][3] = TANB
776 blocks["ALPHA"] = Block("ALPHA")
777 blocks["ALPHA"].set_value(ALPHAH)
778 #
779 ## Neutralino mixing matrix
780 blocks["NMIX"] = Block("NMIX")
781 for i in xrange(1, 5):
782 nmix_i = getnextvalidlineitems()
783 for j, v in enumerate(nmix_i):
784 blocks["NMIX"][i, j+1] = v
785 #
786 ## Chargino mixing matrices V and U
787 blocks["VMIX"] = Block("VMIX")
788 vmix = getnextvalidlineitems()
789 blocks["VMIX"][1, 1] = vmix[0]
790 blocks["VMIX"][1, 2] = vmix[1]
791 blocks["VMIX"][2, 1] = vmix[2]
792 blocks["VMIX"][2, 2] = vmix[3]
793 blocks["UMIX"] = Block("UMIX")
794 umix = getnextvalidlineitems()
795 blocks["UMIX"][1, 1] = umix[0]
796 blocks["UMIX"][1, 2] = umix[1]
797 blocks["UMIX"][2, 1] = umix[2]
798 blocks["UMIX"][2, 2] = umix[3]
799 #
800 THETAT, THETAB, THETAL = getnextvalidlineitems()
801 import math
802 blocks["STOPMIX"] = Block("STOPMIX")
803 blocks["STOPMIX"][1, 1] = math.cos(THETAT)
804 blocks["STOPMIX"][1, 2] = -math.sin(THETAT)
805 blocks["STOPMIX"][2, 1] = math.sin(THETAT)
806 blocks["STOPMIX"][2, 2] = math.cos(THETAT)
807 blocks["SBOTMIX"] = Block("SBOTMIX")
808 blocks["SBOTMIX"][1, 1] = math.cos(THETAB)
809 blocks["SBOTMIX"][1, 2] = -math.sin(THETAB)
810 blocks["SBOTMIX"][2, 1] = math.sin(THETAB)
811 blocks["SBOTMIX"][2, 2] = math.cos(THETAB)
812 blocks["STAUMIX"] = Block("STAUMIX")
813 blocks["STAUMIX"][1, 1] = math.cos(THETAL)
814 blocks["STAUMIX"][1, 2] = -math.sin(THETAL)
815 blocks["STAUMIX"][2, 1] = math.sin(THETAL)
816 blocks["STAUMIX"][2, 2] = math.cos(THETAL)
817 #
818 ATSS, ABSS, ALSS = getnextvalidlineitems()
819 blocks["AU"] = Block("AU")
820 blocks["AU"][3, 3] = ATSS
821 blocks["AD"] = Block("AD")
822 blocks["AD"][3, 3] = ABSS
823 blocks["AE"] = Block("AE")
824 blocks["AE"][3, 3] = ALSS
825 #
826 MUSS = getnextvalidlineitems()[0]
827 blocks["MINPAR"][4] = MUSS
828 #
830 # TODO: Parse RPV boolean and couplings into SLHA2 blocks
832 return blocks, decays
835 def writeISAWIG(blocks, decays, ignorenobr=False, precision=8):
836 """
837 Return a SUSY spectrum definition in the format produced by ISAWIG for inut to HERWIG
838 as a string, from the supplied SLHA blocks and decays dicts.
840 ISAWIG parsing based on the HERWIG SUSY specification format, from
841 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
843 If the ignorenobr parameter is True, do not write decay entries with a
844 branching ratio of zero.
845 """
846 masses = blocks["MASS"]
848 ## Init output string
849 out = ""
851 ## First write out masses section:
852 ## Number of SUSY + top particles
853 ## IDHW, RMASS(IDHW), RLTIM(IDHW)
854 ## repeated for each particle
855 ## IDHW is the HERWIG identity code.
856 ## RMASS and RTLIM are the mass in GeV, and lifetime in seconds respectively.
857 massout = ""
858 for pid in masses.keys():
859 lifetime = -1
860 try:
861 width = decays[pid].totalwidth
862 if width and width > 0:
863 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
864 except:
865 pass
866 massout += ("%d " % pdgid2herwigid(pid)) + _autostr(masses[pid], precision) + " " + _autostr(lifetime, precision) + "\n"
867 out += "%d\n" % massout.count("\n")
868 out += massout
870 assert(len(masses) == len(decays))
872 ## Next each particles decay modes together with their branching ratios and matrix element codes
873 ## Number of decay modes for a given particle (IDK)
874 ## IDK(*), BRFRAC(*), NME(*) & IDKPRD(1-5,*)
875 ## repeated for each mode.
876 ## Repeated for each particle.
877 ## IDK is the HERWIG code for the decaying particle, BRFRAC is the branching ratio of
878 ## the decay mode. NME is a code for the matrix element to be used, either from the
879 ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
880 for i, pid in enumerate(decays.keys()):
881 # if not decays.has_key(pid):
882 # continue
883 hwid = pdgid2herwigid(pid)
884 decayout = ""
885 #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
886 for i_d, d in enumerate(decays[pid].decays):
887 ## Skip decay if it has no branching ratio
888 if ignorenobr and d.br == 0:
889 continue
891 ## Identify decay matrix element to use
892 ## From std HW docs, or from this pair:
893 ## Two new matrix element codes have been added for these new decays:
894 ## NME = 200 3 body top quark via charged Higgs
895 ## 300 3 body R-parity violating gaugino and gluino decays
896 nme = 0
897 # TODO: Get correct condition for using ME 100... this guessed from some ISAWIG output
898 if abs(pid) in (6, 12):
899 nme = 100
900 ## Extra SUSY MEs
901 if len(d.ids) == 3:
902 # TODO: How to determine the conditions for using 200 and 300 MEs? Enumeration of affected decays?
903 pass
904 decayout += ("%d " % hwid) + _autostr(d.br, precision) + (" %d " % nme)
906 def is_quark(pid):
907 return (abs(pid) in range(1, 7))
909 def is_lepton(pid):
910 return (abs(pid) in range(11, 17))
912 def is_squark(pid):
913 if abs(pid) in range(1000001, 1000007):
914 return True
915 if abs(pid) in range(2000001, 2000007):
916 return True
917 return False
919 def is_slepton(pid):
920 if abs(pid) in range(1000011, 1000017):
921 return True
922 if abs(pid) in range(2000011, 2000016, 2):
923 return True
924 return False
926 def is_gaugino(pid):
927 if abs(pid) in range(1000022, 1000026):
928 return True
929 if abs(pid) in (1000035, 1000037):
930 return True
931 return False
933 def is_susy(pid):
934 return (is_squark(pid) or is_slepton(pid) or is_gaugino(pid) or pid == 1000021)
936 absids = map(abs, d.ids)
938 ## Order decay products as required by HERWIG
939 ## Top
940 if abs(pid) == 6:
941 def cmp_bottomlast(a, b):
942 """Comparison function which always puts b/bbar last"""
943 if abs(a) == 5:
944 return True
945 if abs(b) == 5:
946 return False
947 return cmp(a, b)
948 if len(absids) == 2:
949 ## 2 body mode, to Higgs: Higgs; Bottom
950 if (25 in absids or 26 in absids) and 5 in absids:
951 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
952 elif len(absids) == 3:
953 ## 3 body mode, via charged Higgs/W: quarks or leptons from W/Higgs; Bottom
954 if 37 in absids or 23 in absids:
955 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
956 ## Gluino
957 elif abs(pid) == 1000021:
958 if len(absids) == 2:
959 ## 2 body mode
960 ## without gluon: any order
961 ## with gluon: gluon; colour neutral
962 if 21 in absids:
963 def cmp_gluonfirst(a, b):
964 """Comparison function which always puts gluon first"""
965 if a == 21:
966 return False
967 if b == 21:
968 return True
969 return cmp(a, b)
970 d.ids = sorted(d.ids, cmp=cmp_gluonfirst)
971 elif len(absids) == 3:
972 ## 3-body modes, R-parity conserved: colour neutral; q or qbar
973 def cmp_quarkslast(a, b):
974 """Comparison function which always puts quarks last"""
975 if is_quark(a):
976 return True
977 if is_quark(b):
978 return False
979 return cmp(a, b)
980 d.ids = sorted(d.ids, cmp=cmp_quarkslast)
981 ## Squark/Slepton
982 elif is_squark(pid) or is_slepton(pid):
983 def cmp_susy_quark_lepton(a, b):
984 if is_susy(a):
985 return False
986 if is_susy(b):
987 return True
988 if is_quark(a):
989 return False
990 if is_quark(b):
991 return True
992 return cmp(a, b)
993 ## 2 body modes: Gaugino/Gluino with Quark/Lepton Gaugino quark
994 ## Gluino lepton
995 ## 3 body modes: Weak sparticle particles from W decay
996 ## Squark
997 ## 2 body modes: Lepton Number Violated quark lepton
998 ## Baryon number violated quark quark
999 ## Slepton
1000 ## 2 body modes: Lepton Number Violated q or qbar
1001 d.ids = sorted(d.ids, cmp=cmp_bottomlast)
1002 ## Higgs
1003 elif pid in (25, 26):
1004 # TODO: Includes SUSY Higgses?
1005 ## Higgs
1006 ## 2 body modes: (s)quark-(s)qbar (s)q or (s)qbar
1007 ## (s)lepton-(s)lepton (s)l or (s)lbar
1008 ## 3 body modes: colour neutral q or qbar
1009 if len(absids) == 3:
1010 def cmp_quarkslast(a, b):
1011 """Comparison function which always puts quarks last"""
1012 if is_quark(a):
1013 return True
1014 if is_quark(b):
1015 return False
1016 return cmp(a, b)
1017 d.ids = sorted(d.ids, cmp=cmp_quarkslast)
1018 elif is_gaugino(pid):
1019 # TODO: Is there actually anything to do here?
1020 ## Gaugino
1021 ## 2 body modes: Squark-quark q or sq
1022 ## Slepton-lepton l or sl
1023 ##
1024 ## 3 body modes: R-parity conserved colour neutral q or qbar
1025 ## l or lbar
1026 if len(absids) == 3:
1027 def cmp_quarkslast(a, b):
1028 """Comparison function which always puts quarks last"""
1029 if is_quark(a):
1030 return True
1031 if is_quark(b):
1032 return False
1033 return cmp(a, b)
1034 d.ids = sorted(d.ids, cmp=cmp_quarkslast)
1036 # TODO: Gaugino/Gluino
1037 ## 3 body modes: R-parity violating: Particles in the same order as the R-parity violating superpotential
1039 ## Pad out IDs list with zeros
1040 ids = [0,0,0,0,0]
1041 for i, pid in enumerate(d.ids):
1042 ids[i] = pid
1043 ids = map(str, ids)
1044 decayout += " ".join(ids) + "\n"
1045 decayout = "%d\n" % decayout.count("\n") + decayout
1046 out += decayout
1048 ## Now the SUSY parameters
1049 ## TANB, ALPHAH:
1050 out += _autostr(blocks["MINPAR"][3], precision) + " " + _autostr(blocks["ALPHA"].value(), precision) + "\n"
1051 ## Neutralino mixing matrix
1052 nmix = blocks["NMIX"]
1053 for i in xrange(1, 5):
1054 out += _autostr(nmix[i,1], precision) + " " + \
1055 _autostr(nmix[i,2], precision) + " " + \
1056 _autostr(nmix[i,3], precision) + " " + \
1057 _autostr(nmix[i,4], precision) + "\n"
1058 ## Chargino mixing matrices V and U
1059 vmix = blocks["VMIX"]
1060 out += _autostr(vmix[1,1], precision) + " " + \
1061 _autostr(vmix[1,2], precision) + " " + \
1062 _autostr(vmix[2,1], precision) + " " + \
1063 _autostr(vmix[2,2], precision) + "\n"
1064 umix = blocks["UMIX"]
1065 out += _autostr(umix[1,1], precision) + " " + \
1066 _autostr(umix[1,2], precision) + " " + \
1067 _autostr(umix[2,1], precision) + " " + \
1068 _autostr(umix[2,2], precision) + "\n"
1069 ## THETAT,THETAB,THETAL
1070 import math
1071 out += _autostr(math.acos(blocks["STOPMIX"][1,1]), precision) + " " + \
1072 _autostr(math.acos(blocks["SBOTMIX"][1,1]), precision) + " " + \
1073 _autostr(math.acos(blocks["STAUMIX"][1,1]), precision) + "\n"
1074 ## ATSS,ABSS,ALSS
1075 out += _autostr(blocks["AU"][3,3], precision) + " " + \
1076 _autostr(blocks["AD"][3,3], precision) + " " + \
1077 _autostr(blocks["AE"][3,3], precision) + "\n"
1078 ## MUSS == sign(mu)
1079 out += "%f\n" % blocks["MINPAR"][4]
1081 ## RPV SUSY
1082 isRPV = False
1083 out += "%d\n" % isRPV
1084 # TODO: Write RPV couplings if RPV is True (lambda1,2,3; 27 params in each, sci format.
1085 # TODO: Get the index orderings right
1086 # if isRPV: ...
1088 return out
1092 ###############################################################################
1093 ## File-level read/write functions
1095 def readSLHAFile(spcfilename, **kwargs):
1096 """
1097 Read an SLHA file, returning dictionaries of blocks and decays.
1099 Other keyword parameters are passed to readSLHA.
1100 """
1101 f = open(spcfilename, "r")
1102 rtn = readSLHA(f.read(), **kwargs)
1103 f.close()
1104 return rtn
1107 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
1108 """
1109 Write an SLHA file from the supplied blocks and decays dicts.
1111 Other keyword parameters are passed to writeSLHA.
1112 """
1113 f = open(spcfilename, "w")
1114 f.write(writeSLHA(blocks, decays, **kwargs))
1115 f.close()
1118 def readISAWIGFile(isafilename, **kwargs):
1119 """
1120 Read a spectrum definition from a file in the ISAWIG format, returning
1121 dictionaries of blocks and decays. While this is not an SLHA format, it is
1122 informally supported as a useful mechanism for converting ISAWIG spectra to
1123 SLHA.
1125 Other keyword parameters are passed to readSLHA.
1126 """
1127 f = open(isafilename, "r")
1128 rtn = readISAWIG(f.read(), **kwargs)
1129 f.close()
1130 return rtn
1133 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
1134 """
1135 Write an ISAWIG file from the supplied blocks and decays dicts.
1137 Other keyword parameters are passed to writeISAWIG.
1138 """
1139 f = open(isafilename, "w")
1140 f.write(writeISAWIG(blocks, decays, **kwargs))
1141 f.close()
1145 ###############################################################################
1146 ## Main function for module testing
1148 if __name__ == "__main__":
1149 import sys
1150 for a in sys.argv[1:]:
1151 if a.endswith(".isa"):
1152 blocks, decays = readISAWIGFile(a)
1153 else:
1154 blocks, decays = readSLHAFile(a)
1156 for bname, b in sorted(blocks.iteritems()):
1157 print b
1158 print
1160 print blocks.keys()
1162 print blocks["MASS"][25]
1163 print
1165 for p in sorted(decays.values()):
1166 print p
1167 print
1169 print writeSLHA(blocks, decays, ignorenobr=True)
