pyslha.py

changeset 124
7bd52be093b2
parent 123
c85e29bc13c4
child 127
bfd5be72d711
equal deleted inserted replaced
123:c85e29bc13c4 124:7bd52be093b2
260 rtn = readISAWIG(f.read(), kwargs) 260 rtn = readISAWIG(f.read(), kwargs)
261 f.close() 261 f.close()
262 return rtn 262 return rtn
263 263
264 264
265 def writeSLHAFile(spcfilename, blocks, decays, **kwargs):
266 """
267 Write an SLHA file from the supplied blocks and decays dicts.
268
269 Other keyword parameters are passed to writeSLHA.
270 """
271 f = open(spcfilename, "w")
272 f.write(writeSLHA(blocks, decays, kwargs))
273 f.close()
274
275
276 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays}
277
278
279 def writeSLHA(blocks, decays, ignorenobr=False):
280 """
281 Return an SLHA definition as a string, from the supplied blocks and decays dicts.
282 """
283 sep = " "
284 out = ""
285 def dict_hier_strs(d, s=""):
286 if type(d) is dict:
287 for k, v in sorted(d.iteritems()):
288 for s2 in dict_hier_strs(v, s + sep + _autostr(k)):
289 yield s2
290 else:
291 yield s + sep + _autostr(d)
292 ## Blocks
293 for bname, b in sorted(blocks.iteritems()):
294 namestr = b.name
295 if b.q is not None:
296 namestr += " Q= %e" % b.q
297 out += "BLOCK %s\n" % namestr
298 for s in dict_hier_strs(b.entries):
299 out += sep + s + "\n"
300 out += "\n"
301 ## Decays
302 for pid, particle in sorted(decays.iteritems()):
303 out += "DECAY %d %e\n" % (particle.pid, particle.totalwidth or -1)
304 for d in sorted(particle.decays):
305 if d.br > 0.0 or not ignorenobr:
306 products_str = " ".join(map(str, d.ids))
307 out += sep + "%e" % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n"
308 out += "\n"
309 return out
310
311
312
313 ###############################################################################
314 ## ISAWIG handling
315
316 ## Static array of HERWIG IDHW codes mapped to PDG MC ID codes, based on
317 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
318 ## + the IDPDG array and section 4.13 of the HERWIG manual.
319 _HERWIGID2PDGID = {}
320 _HERWIGID2PDGID[7] = -1
321 _HERWIGID2PDGID[8] = -2
322 _HERWIGID2PDGID[9] = -3
323 _HERWIGID2PDGID[10] = -4
324 _HERWIGID2PDGID[11] = -5
325 _HERWIGID2PDGID[12] = -6
326 _HERWIGID2PDGID[13] = 21
327 _HERWIGID2PDGID[59] = 22
328 _HERWIGID2PDGID[121] = 11
329 _HERWIGID2PDGID[122] = 12
330 _HERWIGID2PDGID[123] = 13
331 _HERWIGID2PDGID[124] = 14
332 _HERWIGID2PDGID[125] = 15
333 _HERWIGID2PDGID[126] = 16
334 _HERWIGID2PDGID[127] = -11
335 _HERWIGID2PDGID[128] = -12
336 _HERWIGID2PDGID[129] = -13
337 _HERWIGID2PDGID[130] = -14
338 _HERWIGID2PDGID[131] = -15
339 _HERWIGID2PDGID[132] = -16
340 _HERWIGID2PDGID[203] = 25 ## HIGGSL0 (== PDG standard in this direction)
341 _HERWIGID2PDGID[204] = 35 ## HIGGSH0
342 _HERWIGID2PDGID[205] = 36 ## HIGGSA0
343 _HERWIGID2PDGID[206] = 37 ## HIGGS+
344 _HERWIGID2PDGID[207] = -37 ## HIGGS-
345 _HERWIGID2PDGID[401] = 1000001 ## SSDLBR
346 _HERWIGID2PDGID[407] = -1000001 ## SSDLBR
347 _HERWIGID2PDGID[402] = 1000002 ## SSULBR
348 _HERWIGID2PDGID[408] = -1000002 ## SSUL
349 _HERWIGID2PDGID[403] = 1000003 ## SSSLBR
350 _HERWIGID2PDGID[409] = -1000003 ## SSSL
351 _HERWIGID2PDGID[404] = 1000004 ## SSCLBR
352 _HERWIGID2PDGID[410] = -1000004 ## SSCL
353 _HERWIGID2PDGID[405] = 1000005 ## SSB1BR
354 _HERWIGID2PDGID[411] = -1000005 ## SSB1
355 _HERWIGID2PDGID[406] = 1000006 ## SST1BR
356 _HERWIGID2PDGID[412] = -1000006 ## SST1
357 _HERWIGID2PDGID[413] = 2000001 ## SSDR
358 _HERWIGID2PDGID[419] = -2000001 ## SSDRBR
359 _HERWIGID2PDGID[414] = 2000002 ## SSUR
360 _HERWIGID2PDGID[420] = -2000002 ## SSURBR
361 _HERWIGID2PDGID[415] = 2000003 ## SSSR
362 _HERWIGID2PDGID[421] = -2000003 ## SSSRBR
363 _HERWIGID2PDGID[416] = 2000004 ## SSCR
364 _HERWIGID2PDGID[422] = -2000004 ## SSCRBR
365 _HERWIGID2PDGID[417] = 2000005 ## SSB2
366 _HERWIGID2PDGID[423] = -2000005 ## SSB2BR
367 _HERWIGID2PDGID[418] = 2000006 ## SST2
368 _HERWIGID2PDGID[424] = -2000006 ## SST2BR
369 _HERWIGID2PDGID[425] = 1000011 ## SSEL-
370 _HERWIGID2PDGID[431] = -1000011 ## SSEL+
371 _HERWIGID2PDGID[426] = 1000012 ## SSNUEL
372 _HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
373 _HERWIGID2PDGID[427] = 1000013 ## SSMUL-
374 _HERWIGID2PDGID[433] = -1000013 ## SSMUL+
375 _HERWIGID2PDGID[428] = 1000014 ## SSNUMUL
376 _HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
377 _HERWIGID2PDGID[429] = 1000015 ## SSTAU1-
378 _HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
379 _HERWIGID2PDGID[430] = 1000016 ## SSNUTL
380 _HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
381 _HERWIGID2PDGID[437] = 2000011 ## SSEL-
382 _HERWIGID2PDGID[443] = -2000011 ## SSEL+
383 _HERWIGID2PDGID[438] = 2000012 ## SSNUEL
384 _HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
385 _HERWIGID2PDGID[439] = 2000013 ## SSMUL-
386 _HERWIGID2PDGID[445] = -2000013 ## SSMUL+
387 _HERWIGID2PDGID[440] = 2000014 ## SSNUMUL
388 _HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
389 _HERWIGID2PDGID[441] = 2000015 ## SSTAU1-
390 _HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
391 _HERWIGID2PDGID[442] = 2000016 ## SSNUTL
392 _HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
393 _HERWIGID2PDGID[449] = 1000021 ## GLUINO
394 _HERWIGID2PDGID[450] = 1000022 ## NTLINO1
395 _HERWIGID2PDGID[451] = 1000023 ## NTLINO2
396 _HERWIGID2PDGID[452] = 1000025 ## NTLINO3
397 _HERWIGID2PDGID[453] = 1000035 ## NTLINO4
398 _HERWIGID2PDGID[454] = 1000024 ## CHGINO1+
399 _HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
400 _HERWIGID2PDGID[455] = 1000037 ## CHGINO2+
401 _HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
402 _HERWIGID2PDGID[458] = 1000039 ## GRAVTINO
403
404 def herwigid2pdgid(hwid):
405 """
406 Convert a particle ID code in the HERWIG internal IDHW format (as used by
407 ISAWIG) into its equivalent in the standard PDG ID code definition.
408 """
409 return _HERWIGID2PDGID.get(hwid, hwid)
410
265 411
266 def readISAWIG(isastr, ignorenobr=False): 412 def readISAWIG(isastr, ignorenobr=False):
267 """ 413 """
268 Read a spectrum definition from a string in the ISAWIG format, returning 414 Read a spectrum definition from a string in the ISAWIG format, returning
269 dictionaries of blocks and decays. While this is not an SLHA format, it is 415 dictionaries of blocks and decays. While this is not an SLHA format, it is
275 421
276 If the ignorenobr parameter is True, do not store decay entries with a 422 If the ignorenobr parameter is True, do not store decay entries with a
277 branching ratio of zero. 423 branching ratio of zero.
278 """ 424 """
279 425
280 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
281 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
282 ## + the IDPDG array and section 4.13 of the HERWIG manual.
283 HERWIGID2PDGID = {}
284 HERWIGID2PDGID[7] = -1
285 HERWIGID2PDGID[8] = -2
286 HERWIGID2PDGID[9] = -3
287 HERWIGID2PDGID[10] = -4
288 HERWIGID2PDGID[11] = -5
289 HERWIGID2PDGID[12] = -6
290 HERWIGID2PDGID[13] = 21
291 HERWIGID2PDGID[59] = 22
292 HERWIGID2PDGID[121] = 11
293 HERWIGID2PDGID[122] = 12
294 HERWIGID2PDGID[123] = 13
295 HERWIGID2PDGID[124] = 14
296 HERWIGID2PDGID[125] = 15
297 HERWIGID2PDGID[126] = 16
298 HERWIGID2PDGID[127] = -11
299 HERWIGID2PDGID[128] = -12
300 HERWIGID2PDGID[129] = -13
301 HERWIGID2PDGID[130] = -14
302 HERWIGID2PDGID[131] = -15
303 HERWIGID2PDGID[132] = -16
304 HERWIGID2PDGID[203] = 25 ## HIGGSL0 (== PDG standard in this direction)
305 HERWIGID2PDGID[204] = 35 ## HIGGSH0
306 HERWIGID2PDGID[205] = 36 ## HIGGSA0
307 HERWIGID2PDGID[206] = 37 ## HIGGS+
308 HERWIGID2PDGID[207] = -37 ## HIGGS-
309 HERWIGID2PDGID[401] = 1000001 ## SSDLBR
310 HERWIGID2PDGID[407] = -1000001 ## SSDLBR
311 HERWIGID2PDGID[402] = 1000002 ## SSULBR
312 HERWIGID2PDGID[408] = -1000002 ## SSUL
313 HERWIGID2PDGID[403] = 1000003 ## SSSLBR
314 HERWIGID2PDGID[409] = -1000003 ## SSSL
315 HERWIGID2PDGID[404] = 1000004 ## SSCLBR
316 HERWIGID2PDGID[410] = -1000004 ## SSCL
317 HERWIGID2PDGID[405] = 1000005 ## SSB1BR
318 HERWIGID2PDGID[411] = -1000005 ## SSB1
319 HERWIGID2PDGID[406] = 1000006 ## SST1BR
320 HERWIGID2PDGID[412] = -1000006 ## SST1
321 HERWIGID2PDGID[413] = 2000001 ## SSDR
322 HERWIGID2PDGID[419] = -2000001 ## SSDRBR
323 HERWIGID2PDGID[414] = 2000002 ## SSUR
324 HERWIGID2PDGID[420] = -2000002 ## SSURBR
325 HERWIGID2PDGID[415] = 2000003 ## SSSR
326 HERWIGID2PDGID[421] = -2000003 ## SSSRBR
327 HERWIGID2PDGID[416] = 2000004 ## SSCR
328 HERWIGID2PDGID[422] = -2000004 ## SSCRBR
329 HERWIGID2PDGID[417] = 2000005 ## SSB2
330 HERWIGID2PDGID[423] = -2000005 ## SSB2BR
331 HERWIGID2PDGID[418] = 2000006 ## SST2
332 HERWIGID2PDGID[424] = -2000006 ## SST2BR
333 HERWIGID2PDGID[425] = 1000011 ## SSEL-
334 HERWIGID2PDGID[431] = -1000011 ## SSEL+
335 HERWIGID2PDGID[426] = 1000012 ## SSNUEL
336 HERWIGID2PDGID[432] = -1000012 ## SSNUELBR
337 HERWIGID2PDGID[427] = 1000013 ## SSMUL-
338 HERWIGID2PDGID[433] = -1000013 ## SSMUL+
339 HERWIGID2PDGID[428] = 1000014 ## SSNUMUL
340 HERWIGID2PDGID[434] = -1000014 ## SSNUMLBR
341 HERWIGID2PDGID[429] = 1000015 ## SSTAU1-
342 HERWIGID2PDGID[435] = -1000015 ## SSTAU1+
343 HERWIGID2PDGID[430] = 1000016 ## SSNUTL
344 HERWIGID2PDGID[436] = -1000016 ## SSNUTLBR
345 HERWIGID2PDGID[437] = 2000011 ## SSEL-
346 HERWIGID2PDGID[443] = -2000011 ## SSEL+
347 HERWIGID2PDGID[438] = 2000012 ## SSNUEL
348 HERWIGID2PDGID[444] = -2000012 ## SSNUELBR
349 HERWIGID2PDGID[439] = 2000013 ## SSMUL-
350 HERWIGID2PDGID[445] = -2000013 ## SSMUL+
351 HERWIGID2PDGID[440] = 2000014 ## SSNUMUL
352 HERWIGID2PDGID[446] = -2000014 ## SSNUMLBR
353 HERWIGID2PDGID[441] = 2000015 ## SSTAU1-
354 HERWIGID2PDGID[447] = -2000015 ## SSTAU1+
355 HERWIGID2PDGID[442] = 2000016 ## SSNUTL
356 HERWIGID2PDGID[448] = -2000016 ## SSNUTLBR
357 HERWIGID2PDGID[449] = 1000021 ## GLUINO
358 HERWIGID2PDGID[450] = 1000022 ## NTLINO1
359 HERWIGID2PDGID[451] = 1000023 ## NTLINO2
360 HERWIGID2PDGID[452] = 1000025 ## NTLINO3
361 HERWIGID2PDGID[453] = 1000035 ## NTLINO4
362 HERWIGID2PDGID[454] = 1000024 ## CHGINO1+
363 HERWIGID2PDGID[456] = -1000024 ## CHGINO1-
364 HERWIGID2PDGID[455] = 1000037 ## CHGINO2+
365 HERWIGID2PDGID[457] = -1000037 ## CHGINO2-
366 HERWIGID2PDGID[458] = 1000039 ## GRAVTINO
367 426
368 blocks = {} 427 blocks = {}
369 decays = {} 428 decays = {}
370 LINES = isastr.splitlines() 429 LINES = isastr.splitlines()
371 430
389 masses = Block("MASS") 448 masses = Block("MASS")
390 numentries = int(getnextvalidline()) 449 numentries = int(getnextvalidline())
391 for i in xrange(numentries): 450 for i in xrange(numentries):
392 hwid, mass, lifetime = getnextvalidlineitems() 451 hwid, mass, lifetime = getnextvalidlineitems()
393 width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds 452 width = 1.0/(lifetime * 1.51926778e24) ## width in GeV == hbar/lifetime in seconds
394 pdgid = HERWIGID2PDGID.get(hwid, hwid) 453 pdgid = herwigid2pdgid(hwid)
395 masses.add_entry((pdgid, mass)) 454 masses.add_entry((pdgid, mass))
396 decays[pdgid] = Particle(pdgid, width, mass) 455 decays[pdgid] = Particle(pdgid, width, mass)
397 #print pdgid, mass, width 456 #print pdgid, mass, width
398 blocks["MASS"] = masses 457 blocks["MASS"] = masses
399 458
402 numdecays = int(getnextvalidline()) 461 numdecays = int(getnextvalidline())
403 for d in xrange(numdecays): 462 for d in xrange(numdecays):
404 #print n, numentries-1, d, numdecays-1 463 #print n, numentries-1, d, numdecays-1
405 decayitems = getnextvalidlineitems() 464 decayitems = getnextvalidlineitems()
406 hwid = decayitems[0] 465 hwid = decayitems[0]
407 pdgid = HERWIGID2PDGID.get(hwid, hwid) 466 pdgid = herwigid2pdgid(hwid)
408 br = decayitems[1] 467 br = decayitems[1]
409 nme = decayitems[2] 468 nme = decayitems[2]
410 daughter_hwids = decayitems[3:] 469 daughter_hwids = decayitems[3:]
411 daughter_pdgids = [] 470 daughter_pdgids = []
412 for hw in daughter_hwids: 471 for hw in daughter_hwids:
413 if hw != 0: 472 if hw != 0:
414 daughter_pdgids.append(HERWIGID2PDGID.get(hw, hw)) 473 daughter_pdgids.append(herwigid2pdgid(hw))
415 if not decays.has_key(pdgid): 474 if not decays.has_key(pdgid):
416 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid) 475 #print "Decay for unlisted particle %d, %d" % (hwid, pdgid)
417 decays[pdgid] = Particle(pdgid) 476 decays[pdgid] = Particle(pdgid)
418 decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids) 477 decays[pdgid].add_decay(br, len(daughter_pdgids), daughter_pdgids)
419 478
476 blocks["MINPAR"].add_entry((4, MUSS)) 535 blocks["MINPAR"].add_entry((4, MUSS))
477 # 536 #
478 return blocks, decays 537 return blocks, decays
479 538
480 539
481 def writeSLHAFile(spcfilename, blocks, decays, **kwargs): 540 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
482 """ 541 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
483 Write an SLHA file from the supplied blocks and decays dicts. 542 ## + the IDPDG array and section 4.13 of the HERWIG manual.
484 543 _PDGID2HERWIGID = {}
485 Other keyword parameters are passed to writeSLHA. 544 _PDGID2HERWIGID[ -1] = 7
486 """ 545 _PDGID2HERWIGID[ -2] = 8
487 f = open(spcfilename, "w") 546 _PDGID2HERWIGID[ -3] = 9
488 f.write(writeSLHA(blocks, decays, kwargs)) 547 _PDGID2HERWIGID[ -4] = 10
489 f.close() 548 _PDGID2HERWIGID[ -5] = 11
490 549 _PDGID2HERWIGID[ -6] = 12
491 550 _PDGID2HERWIGID[ 21] = 13
492 # TODO: Split writeSLHA into writeSLHA{Blocks,Decays} 551 _PDGID2HERWIGID[ 22] = 59
493 552 _PDGID2HERWIGID[ 11] = 121
494 553 _PDGID2HERWIGID[ 12] = 122
495 def writeSLHA(blocks, decays, ignorenobr=False): 554 _PDGID2HERWIGID[ 13] = 123
496 """ 555 _PDGID2HERWIGID[ 14] = 124
497 Return an SLHA definition as a string, from the supplied blocks and decays dicts. 556 _PDGID2HERWIGID[ 15] = 125
498 """ 557 _PDGID2HERWIGID[ 16] = 126
499 sep = " " 558 _PDGID2HERWIGID[ -11] = 127
500 out = "" 559 _PDGID2HERWIGID[ -12] = 128
501 def dict_hier_strs(d, s=""): 560 _PDGID2HERWIGID[ -13] = 129
502 if type(d) is dict: 561 _PDGID2HERWIGID[ -14] = 130
503 for k, v in sorted(d.iteritems()): 562 _PDGID2HERWIGID[ -15] = 131
504 for s2 in dict_hier_strs(v, s + sep + _autostr(k)): 563 _PDGID2HERWIGID[ -16] = 132
505 yield s2 564 _PDGID2HERWIGID[ 25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW)
506 else: 565 _PDGID2HERWIGID[ 26] = 203 ## HIGGSL0
507 yield s + sep + _autostr(d) 566 _PDGID2HERWIGID[ 35] = 204 ## HIGGSH0
508 ## Blocks 567 _PDGID2HERWIGID[ 36] = 205 ## HIGGSA0
509 for bname, b in sorted(blocks.iteritems()): 568 _PDGID2HERWIGID[ 37] = 206 ## HIGGS+
510 namestr = b.name 569 _PDGID2HERWIGID[ -37] = 207 ## HIGGS-
511 if b.q is not None: 570 _PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
512 namestr += " Q= %e" % b.q 571 _PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
513 out += "BLOCK %s\n" % namestr 572 _PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
514 for s in dict_hier_strs(b.entries): 573 _PDGID2HERWIGID[-1000002] = 408 ## SSUL
515 out += sep + s + "\n" 574 _PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
516 out += "\n" 575 _PDGID2HERWIGID[-1000003] = 409 ## SSSL
517 ## Decays 576 _PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
518 for pid, particle in sorted(decays.iteritems()): 577 _PDGID2HERWIGID[-1000004] = 410 ## SSCL
519 out += "DECAY %d %e\n" % (particle.pid, particle.totalwidth or -1) 578 _PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
520 for d in sorted(particle.decays): 579 _PDGID2HERWIGID[-1000005] = 411 ## SSB1
521 if d.br > 0.0 or not ignorenobr: 580 _PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
522 products_str = " ".join(map(str, d.ids)) 581 _PDGID2HERWIGID[-1000006] = 412 ## SST1
523 out += sep + "%e" % d.br + sep + "%d" % len(d.ids) + sep + products_str + "\n" 582 _PDGID2HERWIGID[ 2000001] = 413 ## SSDR
524 out += "\n" 583 _PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
525 return out 584 _PDGID2HERWIGID[ 2000002] = 414 ## SSUR
585 _PDGID2HERWIGID[-2000002] = 420 ## SSURBR
586 _PDGID2HERWIGID[ 2000003] = 415 ## SSSR
587 _PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
588 _PDGID2HERWIGID[ 2000004] = 416 ## SSCR
589 _PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
590 _PDGID2HERWIGID[ 2000005] = 417 ## SSB2
591 _PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
592 _PDGID2HERWIGID[ 2000006] = 418 ## SST2
593 _PDGID2HERWIGID[-2000006] = 424 ## SST2BR
594 _PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
595 _PDGID2HERWIGID[-1000011] = 431 ## SSEL+
596 _PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
597 _PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
598 _PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
599 _PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
600 _PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
601 _PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
602 _PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
603 _PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
604 _PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
605 _PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
606 _PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
607 _PDGID2HERWIGID[-2000011] = 443 ## SSEL+
608 _PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
609 _PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
610 _PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
611 _PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
612 _PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
613 _PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
614 _PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
615 _PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
616 _PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
617 _PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
618 _PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
619 _PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
620 _PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
621 _PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
622 _PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
623 _PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
624 _PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
625 _PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
626 _PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
627 _PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
628
629 def pdgid2herwigid(pdgid):
630 """
631 Convert a particle ID code in the standard PDG ID code definition into
632 its equivalent in the HERWIG internal IDHW format (as used by ISAWIG).
633 """
634 return _PDGID2HERWIGID.get(pdgid, pdgid)
635
526 636
527 637
528 638
529 def writeISAWIGFile(isafilename, blocks, decays, **kwargs): 639 def writeISAWIGFile(isafilename, blocks, decays, **kwargs):
530 """ 640 """
547 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html 657 http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/file.html
548 658
549 If the ignorenobr parameter is True, do not write decay entries with a 659 If the ignorenobr parameter is True, do not write decay entries with a
550 branching ratio of zero. 660 branching ratio of zero.
551 """ 661 """
552 ## PDG MC ID codes mapped to HERWIG IDHW codes, based on
553 ## http://www.hep.phy.cam.ac.uk/~richardn/HERWIG/ISAWIG/susycodes.html
554 ## + the IDPDG array and section 4.13 of the HERWIG manual.
555 PDGID2HERWIGID = {}
556 PDGID2HERWIGID[ -1] = 7
557 PDGID2HERWIGID[ -2] = 8
558 PDGID2HERWIGID[ -3] = 9
559 PDGID2HERWIGID[ -4] = 10
560 PDGID2HERWIGID[ -5] = 11
561 PDGID2HERWIGID[ -6] = 12
562 PDGID2HERWIGID[ 21] = 13
563 PDGID2HERWIGID[ 22] = 59
564 PDGID2HERWIGID[ 11] = 121
565 PDGID2HERWIGID[ 12] = 122
566 PDGID2HERWIGID[ 13] = 123
567 PDGID2HERWIGID[ 14] = 124
568 PDGID2HERWIGID[ 15] = 125
569 PDGID2HERWIGID[ 16] = 126
570 PDGID2HERWIGID[ -11] = 127
571 PDGID2HERWIGID[ -12] = 128
572 PDGID2HERWIGID[ -13] = 129
573 PDGID2HERWIGID[ -14] = 130
574 PDGID2HERWIGID[ -15] = 131
575 PDGID2HERWIGID[ -16] = 132
576 PDGID2HERWIGID[ 25] = 203 ## HIGGSL0 (added for PDG standard -> HERWIG IDHW)
577 PDGID2HERWIGID[ 26] = 203 ## HIGGSL0
578 PDGID2HERWIGID[ 35] = 204 ## HIGGSH0
579 PDGID2HERWIGID[ 36] = 205 ## HIGGSA0
580 PDGID2HERWIGID[ 37] = 206 ## HIGGS+
581 PDGID2HERWIGID[ -37] = 207 ## HIGGS-
582 PDGID2HERWIGID[ 1000001] = 401 ## SSDLBR
583 PDGID2HERWIGID[-1000001] = 407 ## SSDLBR
584 PDGID2HERWIGID[ 1000002] = 402 ## SSULBR
585 PDGID2HERWIGID[-1000002] = 408 ## SSUL
586 PDGID2HERWIGID[ 1000003] = 403 ## SSSLBR
587 PDGID2HERWIGID[-1000003] = 409 ## SSSL
588 PDGID2HERWIGID[ 1000004] = 404 ## SSCLBR
589 PDGID2HERWIGID[-1000004] = 410 ## SSCL
590 PDGID2HERWIGID[ 1000005] = 405 ## SSB1BR
591 PDGID2HERWIGID[-1000005] = 411 ## SSB1
592 PDGID2HERWIGID[ 1000006] = 406 ## SST1BR
593 PDGID2HERWIGID[-1000006] = 412 ## SST1
594 PDGID2HERWIGID[ 2000001] = 413 ## SSDR
595 PDGID2HERWIGID[-2000001] = 419 ## SSDRBR
596 PDGID2HERWIGID[ 2000002] = 414 ## SSUR
597 PDGID2HERWIGID[-2000002] = 420 ## SSURBR
598 PDGID2HERWIGID[ 2000003] = 415 ## SSSR
599 PDGID2HERWIGID[-2000003] = 421 ## SSSRBR
600 PDGID2HERWIGID[ 2000004] = 416 ## SSCR
601 PDGID2HERWIGID[-2000004] = 422 ## SSCRBR
602 PDGID2HERWIGID[ 2000005] = 417 ## SSB2
603 PDGID2HERWIGID[-2000005] = 423 ## SSB2BR
604 PDGID2HERWIGID[ 2000006] = 418 ## SST2
605 PDGID2HERWIGID[-2000006] = 424 ## SST2BR
606 PDGID2HERWIGID[ 1000011] = 425 ## SSEL-
607 PDGID2HERWIGID[-1000011] = 431 ## SSEL+
608 PDGID2HERWIGID[ 1000012] = 426 ## SSNUEL
609 PDGID2HERWIGID[-1000012] = 432 ## SSNUELBR
610 PDGID2HERWIGID[ 1000013] = 427 ## SSMUL-
611 PDGID2HERWIGID[-1000013] = 433 ## SSMUL+
612 PDGID2HERWIGID[ 1000014] = 428 ## SSNUMUL
613 PDGID2HERWIGID[-1000014] = 434 ## SSNUMLBR
614 PDGID2HERWIGID[ 1000015] = 429 ## SSTAU1-
615 PDGID2HERWIGID[-1000015] = 435 ## SSTAU1+
616 PDGID2HERWIGID[ 1000016] = 430 ## SSNUTL
617 PDGID2HERWIGID[-1000016] = 436 ## SSNUTLBR
618 PDGID2HERWIGID[ 2000011] = 437 ## SSEL-
619 PDGID2HERWIGID[-2000011] = 443 ## SSEL+
620 PDGID2HERWIGID[ 2000012] = 438 ## SSNUEL
621 PDGID2HERWIGID[-2000012] = 444 ## SSNUELBR
622 PDGID2HERWIGID[ 2000013] = 439 ## SSMUL-
623 PDGID2HERWIGID[-2000013] = 445 ## SSMUL+
624 PDGID2HERWIGID[ 2000014] = 440 ## SSNUMUL
625 PDGID2HERWIGID[-2000014] = 446 ## SSNUMLBR
626 PDGID2HERWIGID[ 2000015] = 441 ## SSTAU1-
627 PDGID2HERWIGID[-2000015] = 447 ## SSTAU1+
628 PDGID2HERWIGID[ 2000016] = 442 ## SSNUTL
629 PDGID2HERWIGID[-2000016] = 448 ## SSNUTLBR
630 PDGID2HERWIGID[ 1000021] = 449 ## GLUINO
631 PDGID2HERWIGID[ 1000022] = 450 ## NTLINO1
632 PDGID2HERWIGID[ 1000023] = 451 ## NTLINO2
633 PDGID2HERWIGID[ 1000025] = 452 ## NTLINO3
634 PDGID2HERWIGID[ 1000035] = 453 ## NTLINO4
635 PDGID2HERWIGID[ 1000024] = 454 ## CHGINO1+
636 PDGID2HERWIGID[-1000024] = 456 ## CHGINO1-
637 PDGID2HERWIGID[ 1000037] = 455 ## CHGINO2+
638 PDGID2HERWIGID[-1000037] = 457 ## CHGINO2-
639 PDGID2HERWIGID[ 1000039] = 458 ## GRAVTINO
640 662
641 masses = blocks["MASS"].entries 663 masses = blocks["MASS"].entries
642 664
643 ## Init output string 665 ## Init output string
644 out = "" 666 out = ""
656 width = decays[pid].totalwidth 678 width = decays[pid].totalwidth
657 if width and width > 0: 679 if width and width > 0:
658 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV 680 lifetime = 1.0/(width * 1.51926778e24) ## lifetime in seconds == hbar/width in GeV
659 except: 681 except:
660 pass 682 pass
661 massout += "%d %e %e\n" % (PDGID2HERWIGID.get(pid, pid), masses[pid], lifetime) 683 massout += "%d %e %e\n" % (pdgid2herwigid(pid), masses[pid], lifetime)
662 out += "%d\n" % massout.count("\n") 684 out += "%d\n" % massout.count("\n")
663 out += massout 685 out += massout
664 686
665 assert(len(masses) == len(decays)) 687 assert(len(masses) == len(decays))
666 688
673 ## the decay mode. NME is a code for the matrix element to be used, either from the 695 ## the decay mode. NME is a code for the matrix element to be used, either from the
674 ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products. 696 ## SUSY elements or the main HERWIG MEs. IDKPRD are the HERWIG identity codes of the decay products.
675 for i, pid in enumerate(decays.keys()): 697 for i, pid in enumerate(decays.keys()):
676 # if not decays.has_key(pid): 698 # if not decays.has_key(pid):
677 # continue 699 # continue
678 hwid = PDGID2HERWIGID.get(pid, pid) 700 hwid = pdgid2herwigid(pid)
679 decayout = "" 701 decayout = ""
680 #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid) 702 #decayout += "@@@@ %d %d %d\n" % (i, pid, hwid)
681 for i_d, d in enumerate(decays[pid].decays): 703 for i_d, d in enumerate(decays[pid].decays):
682 ## Skip decay if it has no branching ratio 704 ## Skip decay if it has no branching ratio
683 if ignorenobr and d.br == 0: 705 if ignorenobr and d.br == 0:

mercurial