Fro version 1.3:
COMPLEX FUNCTION XXX\*16(xxx) to COMPLEX\*16 FUNCTION XXX(xxx)REWIND(ninlha)#!/usr/bin/python
import re
from subprocess import call
from copy import deepcopy
import os.path
import shutil
import matplotlib.pyplot as mpl
class SusyHitPlotter:
def __init__(self,options,debug=False):
self.runCommand = ["./run"]
self.options = options
self.debug = debug
def changeAllUnlistedOptions(self,outfile,mapToChange):
for option in mapToChange:
value = mapToChange[option]
if value != "written":
outfile.write(" {0} {1} #Added Option by susyHitPlotter.py\n".format(option,value))
def setupInputs(self,mapToChange):
mapToChange = deepcopy(mapToChange)
print("Run Options: "+str(mapToChange))
if not os.path.isfile("suspect2_lha.in.bak"):
shutil.copy("suspect2_lha.in", "suspect2_lha.in.bak")
infile = open("suspect2_lha.in.bak",'r')
outfile = open("suspect2_lha.in",'w')
currentBlock = ""
for line in infile:
blockMatch = re.search(r"Block[\s]+([\w]+)",line)
optionMatch = re.search(r"[\s]*([\d]+)[\s]+([.\ddeDE+-]+)",line)
#print(line[0])
if(line[0] == '#'):
outfile.write(line)
elif(blockMatch):
#print("block: "+blockMatch.group(1))
tmpBlock = blockMatch.group(1)
if(mapToChange.has_key(currentBlock)):
#print("changeingUnlistedOptions for block: {0}".format(currentBlock))
self.changeAllUnlistedOptions(outfile,mapToChange[currentBlock])
currentBlock = tmpBlock
outfile.write(line)
elif(optionMatch):
#print("option: " + optionMatch.group(1)+" "+optionMatch.group(2))
if(mapToChange.has_key(currentBlock)):
option = optionMatch.group(1)
currValue = optionMatch.group(2)
if(mapToChange[currentBlock].has_key(option)):
newValue = mapToChange[currentBlock][option]
outfile.write('#'+line)
outfile.write(" {0} {1} #Changed Option by susyHitPlotter.py\n".format(option,newValue))
outfile.flush()
mapToChange[currentBlock][option] = "written"
else:
outfile.write(line)
else:
outfile.write(line)
else:
print("Warning: Line not identified: "+line)
outfile.write(line)
if(mapToChange.has_key(currentBlock)): # add final unwritten lines
self.changeAllUnlistedOptions(outfile,mapToChange[currentBlock])
infile.close()
outfile.close()
def getOutput(self):
infile = open("susyhit_slha.out",'r')
currentBlock = ""
currentDecay = ""
masses = {}
decays = {}
for line in infile:
blockMatch = re.search(r"BLOCK[\s]+([\w]+)",line)
decayMatch = re.search(r"DECAY[\s]+([\d]+)[\s]+([.\ddeDE+-]+)",line)
massMatch = re.search(r"[\s]*([\d]+)[\s]+([.\ddeDE+-]+)",line)
brMatch = re.search(r"[\s]*([.\ddeDE+-]+)[\s]+([\d]+)[\s]+([\d]+)[\s]+([\d]+)[\s]+",line)
#print(line[0])
if(line[0] == '#'):
pass
if(re.match(r"[\s]*#.*",line)):
pass
elif(blockMatch):
#print("block: "+blockMatch.group(1))
tmpBlock = blockMatch.group(1)
currentBlock = tmpBlock
currentDecay = ""
elif(currentBlock == "MASS" and massMatch):
#print("mass: " + massMatch.group(1)+" "+massMatch.group(2))
pdgid = massMatch.group(1)
massValue = massMatch.group(2)
masses[pdgid] = massValue
elif(decayMatch):
#print("decay: " + decayMatch.group(1)+" "+decayMatch.group(2))
pdgid = massMatch.group(1)
width = massMatch.group(2)
decays[pdgid] = {}
decays[pdgid]["width"] = width
currentDecay = pdgid
elif(currentDecay != "" and brMatch):
#print("br: " + brMatch.group(1)+" "+brMatch.group(2)+" "+brMatch.group(3)+" "+brMatch.group(4))
if(brMatch.group(2) != '2'):
continue
br = brMatch.group(1)
id1 = brMatch.group(3)
id2 = brMatch.group(4)
decays[currentDecay][id1+","+id2] = br
infile.close()
#print decays
return masses, decays
def run(self,optionToChange,massesToWatch,brsToWatch,parameterLabel,logx=False,logy=False):
block = ""
optionToChangeName = ""
for b in optionToChange:
block = b
for o in optionToChange[b]:
optionToChangeName = o
if not self.options.has_key(block):
self.options[block] = {}
massesToWatchData = {}
brsToWatchData = {}
for key in massesToWatch:
massesToWatchData[str(key)] = []
for key in brsToWatch:
brsToWatchData[str(key)] = []
debugOutputFile = ""
if self.debug:
debugOutputFile = open("logRun","w")
else:
debugOutputFile = open("/dev/null","w")
for i in optionToChange[block][optionToChangeName]:
self.options[block][optionToChangeName] = str(i)
self.setupInputs(self.options)
call(self.runCommand,stdout=debugOutputFile,stderr=debugOutputFile)
masses, decays = self.getOutput()
#print("\n\n\n")
#print(decays)
#print("\n\n\n")
for key in massesToWatchData:
if masses.has_key(key):
massesToWatchData[key].append(masses[key])
else:
massesToWatchData[key].append(-10.0)
for key in brsToWatchData:
if decays.has_key(key):
brsToWatchData[key].append(decays[key])
else:
brsToWatchData[key].append(-10.0)
#print("\n\n\n")
#print(brsToWatchData)
#print("\n\n\n")
debugOutputFile.close()
func = getattr(self,"makePlots")
func(optionToChange[block][optionToChangeName],massesToWatchData,brsToWatchData,parameterLabel,logx,logy)
print("Done.")
def getParticleLatex(self,number):
sModel = ["??","d","u","s","c","b","t","b'","t'","b''","t''", #0-10
"e","#nu_{e}","#mu","#nu_{#mu}","#tau","#nu_{#tau}", #11-16
"#tau'","#nu_{#tau}'","#tau''","#nu_{#tau}''", #17-20
"g","#gamma","Z","W","h^{0}" #21-25
]
if type(number) == str:
number = int(number)
number = abs(number)
result = ""
if number < 26:
result = sModel[number]
return result
if number == 39:
return "G"
if number == 41:
return "R^{0}"
if number == 42:
return "LQ"
if number > 2000000:
number = number - 2000000
if number < 5:
return "#tilde{"+sModel[number]+"}_{R}"
elif number == 5 or number == 6 or number == 17:
return "#tilde{"+sModel[number]+"}_{2}"
elif number == 11 or number == 13:
return "#tilde{"+sModel[number]+"}_{R}"
else:
return str(2000000+number)
elif number > 1000000:
number = number - 1000000
if number < 5:
return "#tilde{"+sModel[number]+"}_{L}"
elif number == 5 or number == 6 or number == 17:
return "#tilde{"+sModel[number]+"}_{1}"
elif number == 11 or number == 13:
return "#tilde{"+sModel[number]+"}_{L}"
elif number == 21:
return "#tilde{"+sModel[number]+"}"
elif number == 22:
return "#tilde{#chi_{1}^{0}}"
elif number == 23:
return "#tilde{#chi_{2}^{0}}"
elif number == 24:
return "#tilde{#chi_{1}^{+}}"
elif number == 25:
return "#tilde{#chi_{3}^{0}}"
elif number == 35:
return "#tilde{#chi_{4}^{0}}"
elif number == 37:
return "#tilde{#chi_{2}^{+}}"
elif number == 39:
return "#tilde{G}"
else:
return str(1000000+number)
elif number > 400000:
return sModel[number]+"^{*}"
elif number == 2212:
return "p"
elif number == 2112:
return "n"
return sModel[number]+"^{*}"
def convertRealLatex(self,latex):
latex = re.sub(r"#",r"\\",latex)
return latex
def makePlots(self,xData,massData,brData,xLabel,logx,logy):
x = []
for i in xData:
x.append(float(i))
fig = mpl.figure()
ax1 = fig.add_subplot(111)
ax1.set_xlabel(xLabel)
ax1.set_ylabel("Particle Mass [GeV]")
plotfunc = ax1.plot
if logx and logy:
plotfunc=ax1.loglog
elif logx:
plotfunc=ax1.semilogx
elif logy:
plotfunc=ax1.semilogy
legendEntries = []
for key in massData:
particleLatex = self.getParticleLatex(key)
particleLatex = self.convertRealLatex(particleLatex)
particleLatex = "$"+ particleLatex+"$"
legendEntries.append(particleLatex)
data = []
for i in massData[key]:
data.append(float(i))
if len(data) != len(x):
print("Error massData doesn't have the same length as varying parameter list")
print("{0} != {1}".format(len(data),len(x)))
plot = plotfunc(x,data,label=particleLatex)
ax1.legend()
fig.savefig("masses.png")
fig.savefig("masses.pdf")
for key in brData:
fig.clf()
ax1 = fig.add_subplot(111)
ax1.set_xlabel(xLabel)
ax1.set_ylabel("Particle Mass [GeV]")
plotfunc = ax1.plot
if logx and logy:
plotfunc=ax1.loglog
elif logx:
plotfunc=ax1.semilogx
elif logy:
plotfunc=ax1.semilogy
particleLatex = self.getParticleLatex(key)
particleLatex = self.convertRealLatex(particleLatex)
motherLatex = "$"+ particleLatex+"$"
dictOfDecays = brData[key][0]
#print "mother: {0}".format(key)
#print "dictofDecays: {0}".format(dictOfDecays)
for decay in dictOfDecays:
if decay == "width":
continue
matchObj = re.match(r"([\d]+),([\d]+)",decay)
if not matchObj:
continue
decay1 = matchObj.group(1)
decay2 = matchObj.group(2)
daughter1Latex = self.getParticleLatex(decay1)
daughter1Latex = self.convertRealLatex(daughter1Latex)
daughter1Latex = "$"+daughter1Latex+"$"
daughter2Latex = self.getParticleLatex(decay2)
daughter2Latex = self.convertRealLatex(daughter2Latex)
daughter2Latex = "$"+daughter2Latex+"$"
decayString = motherLatex+r" $\rightarrow$ "+daughter1Latex+""+daughter2Latex
data = []
for i in brData[key]:
data.append(float(i[decay]))
plot = plotfunc(x,data,label=decayString)
ax1.legend()
fig.savefig("decays"+str(key)+".png")
fig.savefig("decays"+str(key)+".pdf")
if __name__ == "__main__":
options ={
"MINPAR": {
"1":"500.0", #m0
"2":"500.0", #m1/2
"3":"10", #tan beta
"5":"0.0" #A0
}
}
#options = {
#"MINPAR": {
# "1":"10000.0", #m0
# "2":"10000.0", #m1/2
# "3":"0", #tan beta
# "5":"0.0" #A0
# },
#"EXTPAR": {
# "1":"100.0",
# "2":"100.0",
# "3":"100.0",
# "31":"5000.0",
# "32":"5000.0",
# "33":"5000.0",
# "34":"5000.0",
# "35":"5000.0",
# "36":"5000.0",
# "41":"100.0",
# "42":"100.0",
# "43":"100.0", #M_q3L
# "44":"5000.0",
# "45":"5000.0",
# "46":"300.0", #M_tR
# "47":"5000.0",
# "48":"5000.0",
# "49":"300.0" #M_bR
# }
#}
#change = {"EXTPAR":{"1":[50.,100.,200.0,500.0,1000.0,2000.0,5000.0]}}
#label = r"$m_{q_3,L}$ [GeV]"
change = {"MINPAR":{"2":[50.,100.,200.0,500.0,1000.0,2000.0,5000.0]}}
label = r"$m_{1/2}$ [GeV]"
#change = {"MINPAR":{"1":[50.,100.,200.0,500.0,1000.0,2000.0,5000.0]}}
#label = r"$m_{0}$ [GeV]"
#change = {"MINPAR":{"3":[0.1,1,10,100]}}
#label = r"$\tan \beta$ [GeV]"
#change = {"MINPAR":{"5":[-1000.,-100.,0.0,100,1000]}}
#label = r"$A_0$ [GeV]"
#change = {"MINPAR":{"5":[0,-100]}}
#label = r"$A_0$ [GeV]"
massesToWatch = [1000022,1000006,2000006,1000005,2000005,1000002,1000013]
susy = SusyHitPlotter(options,debug=False)
susy.run(change,massesToWatch,massesToWatch,label,logy=True)