SOLARIS_Web_Simulation/simpleInput.py
2024-01-04 16:46:19 -05:00

345 lines
11 KiB
Python
Executable File

#!/usr/bin/env /usr/bin/python3
import os
import cgi, cgitb
import subprocess
def FindName(A, Z):
result = subprocess.run(["../Cleopatra/IsotopeShort", str(A), str(Z)], stdout=subprocess.PIPE).stdout.decode('utf-8')
ppp1 = result.find('Name:')
ppp2 = result.find('|')
return result[ppp1+5:ppp2]
def FindAZ(AZ):
result = subprocess.run(["../Cleopatra/IsotopeShort", AZ], stdout=subprocess.PIPE).stdout.decode('utf-8')
pA = result.find('A:')
pZ = result.find('Z:')
return [int(result[pA+2:pA+5]), int(result[pZ+2:pZ+5])]
form = cgi.FieldStorage()
beam_AZ = form.getvalue('beam_AZ')
beam_Jpi = form.getvalue('beam_Jpi')
beam_Ex = form.getvalue('beam_Ex')
target_AZ = form.getvalue('target_AZ')
lRecoil_AZ = form.getvalue('lRecoil_AZ')
numEvent = form.getvalue('numEvent')
KEA = form.getvalue('KEA')
BField = form.getvalue('BField')
posArray = form.getvalue('posArray')
posRecoil = form.getvalue('posRecoil')
isDWBA = form.getvalue('DWBA')
if isDWBA == "On" :
isDWBA = True
else:
isDWBA = False
onlyDWBA = form.getvalue('onlyDWBA')
minAng = form.getvalue('minAng')
maxAng = form.getvalue('maxAng')
if onlyDWBA == "On" :
onlyDWBA = True
isDWBA = True
else :
onlyDWBA = False
SSType = form.getvalue('SSType')
ExList = form.getvalue('Ex')
JpiList = form.getvalue('Jpi')
OrbList = form.getvalue('Orb')
op1 = form.getvalue('op1')
op2 = form.getvalue('op2')
plot = form.getvalue('plot')
gate = form.getvalue('gate')
#====== ch directory
os.chdir(r"files")
#================== write reactionConfig.txt
beam = FindAZ(beam_AZ)
target = FindAZ(target_AZ)
lRecoil = FindAZ(lRecoil_AZ)
fn1 = "sim_reactionConfig.txt"
f = open(fn1, "w")
f.write("%s //beam_A\n" % beam[0])
f.write("%s //beam_Z\n" % beam[1])
f.write("%s //target_A\n" % target[0])
f.write("%s //target_Z\n" % target[1])
f.write("%s //recoil_light_A\n" % lRecoil[0])
f.write("%s //recoil_light_Z\n" % lRecoil[1])
f.write("%.3f //beam_energy_in_MeV/u\n" % float(KEA))
f.write("0.000 //beam_energy_sigma_in_MeV/u\n")
f.write("0.000 //beam_angle_in_mrad\n")
f.write("0.000 //beam_emittance_in_mrad\n")
f.write("0.000 //beam_x-offset_in_mm\n")
f.write("0.000 //beam_y-offset_in_mm\n")
f.write("%s //number_of_event_generated\n" % numEvent)
f.write("false //isTargetScattering\n")
f.write("0.913 //Target_density_in_g/cm3\n")
f.write("2.2e-4 //Target_thickness_in_cm\n")
f.write("../SRIM/20F_in_CD2.txt //stopping_power_for_beam\n")
f.write("../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil\n")
f.write("../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil\n")
f.write("false //isDacay\n")
f.write("32 //decayNucleus_A\n")
f.write("14 //decayNucleus_Z\n")
f.write("false //isReDo\n")
f.write("%s //excitation_energy_of_A[MeV]\n" % beam_Ex)
f.write("#===== end of file")
f.close()
#================== write detectorGeo.txt
bore="462.5"
prepDistant="11.5"
detWidth="10.0"
detLength="50.0"
nSide="4"
pos=['58.6', '117.9', '176.8', '235.8', '294.0']
if SSType == "SOLARIS" :
nSide="6"
pos=['58.6', '117.9', '176.8', '235.8', '294.0', '353.2', '412.1', '471.0', '530.0']
if SSType == "ISS":
nSide="6"
prepDistant="28.75"
detWidth="27"
detLength="125"
pos=['125.5', '251', '376.5']
fn2 = "sim_detectorGeo.txt"
f = open(fn2, "w")
f.write("%s //Bfield_[T]\n" % BField)
f.write("0.00 //Bfield_direction_to_z-axis_[deg]_should_not_use\n")
f.write("%s //bore_[mm]\n" % bore)
f.write("%s //distance_from_axis_[mm]\n" % prepDistant)
f.write("%s //width_of_detector_[mm]\n" % detWidth)
f.write("%s //length_of_detector_[mm]\n" % detLength)
f.write("%s //recoil_position_+_for_downstream_[mm]\n" % posRecoil)
f.write("10.0 //inner_radius_of_recoil_detector_[mm]\n")
f.write("40.2 //outter_radius_of_recoil_detector_[mm]\n")
f.write("false //is_coincident_with_recoil\n")
f.write("0 //Recoil_1_position_[mm]_when_0_disable_tree_branch\n")
f.write("0 //Recoil_2_position_[mm]\n")
f.write("0.00 //Elum_1_position_[mm]_(just_another_recoil_detector_but_for_light_recoil)\n")
f.write("0.00 //Elum_2_position_[mm]_when_Elum=0_disable_tree_branch\n")
f.write("0 //support_length_[mm]\n")
f.write("%s //first_position_-_for_upstream_[mm]\n" % posArray)
f.write("0.03 //energy_resolution_of_PSD_array_[MeV]\n")
f.write("1.00 //position_resolution_of_PSD_array_[mm]\n")
f.write("Out //detector_facing_Out_or_In\n")
f.write("%s //number_of_detector_as_same_side\n" % nSide)
f.write("0.00 //1st_detector_near_position_in_reference_to_det6_[mm]\n")
for a in pos:
f.write("%s\n" % a)
f.write("#============= end of file")
f.close()
#================== if DWBA, write DWBA.txt, else write Ex.txt
fn4 = "sim_DWBA.txt"
hRecoil = [ beam[0] + target[0] - lRecoil[0], beam[1] + target[1] - lRecoil[1] ]
hRecoil_AZ = FindName(hRecoil[0], hRecoil[1]);
if isDWBA :
#===== find the symbol of Beam
reactionName = beam_AZ + "(" + target_AZ +","+lRecoil_AZ+")" + hRecoil_AZ;
f = open(fn4, "w")
if isinstance(ExList, list):
nState = len(ExList)
for i in range(0, nState):
f.write("%s %s %s %s %.3f %sMeV/u %s%s\n" % (reactionName, beam_Jpi, OrbList[i], JpiList[i], float(ExList[i]), KEA, op1, op2))
else :
f.write("%s %s %s %s %.3f %sMeV/u %s%s\n" % (reactionName, beam_Jpi, OrbList, JpiList, float(ExList), KEA, op1, op2))
f.close()
fnEx="sim_Ex.txt"
if isDWBA == False:
f = open(fnEx, "w")
f.write("//Ex rel_xsec SF sigma_in_MeV\n")
if isinstance(ExList, list):
for i in range(0, len(ExList)):
f.write("%.3f 1.0 1.0 0.000\n" % float(ExList[i]))
else:
f.write("%.3f 1.0 1.0 0.000\n" % float(ExList))
f.write("#==============_End_of_file")
f.close()
#================== write PlotConfig.txt (to be done)
fn5="sim_PlotConfig.txt"
plotStr="{"
plotSize=len(plot)
for i in range(plotSize):
plotStr += plot[i] + ", "
if( i == (int)(plotSize/2) and plotSize > 3 ) :
plotStr += "break, "
plotStr +="pInfo}\n"
gateStr =""
gateSize=len(gate)
for i in range(gateSize):
gateStr += gate[i]
if( i < gateSize - 1 ):
gateStr += " && "
gateStr +="\n"
f = open(fn5, "w")
f.write(plotStr)
f.write(gateStr)
f.write("60 //elum range\n")
f.write("{0, 50} //thetaCM range\n")
f.write("false //showKElines\n")
f.write("false //isOverRideEx\n")
f.write("{-0.5, 4.0} //over-ride Ex range\n")
f.close()
#==================== Run Simulation
dwba_1=""
outFile=""
dwba_3=""
dwba_4=""
isPtolemy=False
dwbaPNG=""
dwbaRootFile=""
dwbaExFile=""
pngName=""
result=""
haha=""
if isDWBA :
if onlyDWBA :
dwba_1 = subprocess.run(["../Cleopatra/InFileCreator", fn4, minAng, maxAng, "0.1"] , stdout=subprocess.PIPE).stdout.decode('utf-8')
else :
dwba_1 = subprocess.run(["../Cleopatra/InFileCreator", fn4, "0", "180", "0.1"] , stdout=subprocess.PIPE).stdout.decode('utf-8')
cmd = "../Cleopatra/ptolemy <%s.in> %s.out" % (fn4, fn4)
outFile = fn4+ ".out"
if os.path.exists(outFile):
os.remove(outFile)
os.system(cmd)
if os.path.exists(outFile) :
if os.path.exists(fn4 + ".root"):
os.remove(fn4 + ".root")
if os.path.exists(fn4 + ".Ex.txt"):
os.remove(fn4 + ".Ex.txt")
if os.path.exists(fn4 + ".Xsec.txt"):
os.remove(fn4 + ".Xsec.txt")
dwba_3 = subprocess.run(["../Cleopatra/ExtractXSec", outFile, "1"], stdout=subprocess.PIPE).stdout.decode('utf-8')
dwbaRootFile = fn4 + ".root"
dwbaExFile = fn4 + ".Ex.txt"
dwba_4 = subprocess.run(["../Cleopatra/PlotTGraphTObjArray", dwbaRootFile, "1"], stdout=subprocess.PIPE).stdout.decode('utf-8')
dwbaPNG = dwba_4.splitlines().pop()
#if os.path.exists("transfer.root"):
# os.remove("transfer.root");
if onlyDWBA == False:
if isDWBA :
result = subprocess.run(["../Cleopatra/Transfer", fn1, fn2, dwbaExFile, dwbaRootFile, "transfer.root", "reaction.dat", "0"], stdout=subprocess.PIPE).stdout.decode('utf-8')
else:
result = subprocess.run(["../Cleopatra/Transfer", fn1, fn2, fnEx, "", "transfer.root", "reaction.dat", "0"], stdout=subprocess.PIPE).stdout.decode('utf-8')
if fn5 =="":
haha = subprocess.run(["../Cleopatra/PlotSimulation", "transfer.root"], stdout=subprocess.PIPE).stdout.decode('utf-8')
else:
haha = subprocess.run(["../Cleopatra/PlotSimulation", "transfer.root", fn5], stdout=subprocess.PIPE).stdout.decode('utf-8')
pngName=haha.splitlines().pop()
print ("Content-Type: text/html\r\n\r\n")
print ("<html>")
print ("<head>")
print ("<title>Result of DWBA and Monte Carlo Simulation</title>")
print ("<link rel=\"icon\" type=\"image/x-icon\" href=\"SOLARIS_favicon.png\">")
print ("</head>")
print ("<style> img {max-width: 100%} </style>")
print ("<style> body { font-family: courier, courier new, serif; background: #6DB33E} </style>")
#print ("<style> a { color: #F7CF3C;} </style>")
print ("<style> a { color: #1364A7;} </style>")
print ("<body>")
'''
print("<br> %s" % beam_A)
print("<br> %s" % beam_Z)
print("<br> %s" % beam_Jpi)
print("<br> %s" % beam_Ex)
print("<br> %s" % target_A)
print("<br> %s" % target_Z)
print("<br> %s" % lRecoil_A)
print("<br> %s" % lRecoil_Z)
print("<br> %s" % KEA)
print("<br> %s" % BField)
print("<br> %s" % posArray)
print("<br> %s" % posRecoil)
print("<br> %s" % SSType)
print("<br> do DWBA : |%s|" % isDWBA)
print("<br> %s" % ExList)
print("<br> %s" % JpiList)
print("<br> %s" % OrbList)
print("<br> %s" % op1)
print("<br> %s" % op2)
print("<br> %s" % plot)
print("<br> %s" % gate)
'''
if dwbaPNG != "" :
print ("<img src=\"files/%s\">" % dwbaPNG)
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">PNG</a>" % (dwbaPNG, dwbaPNG))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">%s</a>" % (fn4, fn4, fn4))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">%s</a>" % (fn4+".in", fn4+".in", fn4+".in"))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">%s</a>" % (fn4+".out", fn4+".out", fn4+".out"))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">%s</a>" % (fn4+".Xsec.txt", fn4+".Xsec.txt", fn4+".Xec.txt"))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">%s</a>" % (fn4+".root", fn4+".root", fn4+".root"))
print ("<br>")
print ("<p style=\"color:#1363A7;font-size:20px;\"><b>If the DWBA result different from the setting or not show, probably the DWBA setting is incorrect that causes the calculation fail.</b></p>")
if pngName != "" :
print ("<img src=\"files/%s\">" % pngName)
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">PNG</a>" % (pngName, pngName))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">reactionConfig file</a>" % (fn1, fn1))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">detectorGeo file</a>" % (fn2, fn2))
print ("<br> Download the <a href=\"files/%s\" download=\"%s\">Ex file</a>" % (fnEx, fnEx))
print ("<br> Download the simulation <a href=\"files/transfer.root\" download=\"transfer.root\">root file</a>")
if fn5 != "" :
print ("<br> Downlaod the <a href=\"files/%s\" download=\"%s\">Plot Config file</a>" % (fn5, fn5))
print ("<br><br><br>#################################### programs text output")
if isDWBA:
print ("<br><span style=\"white-space: pre-line\"> %s </span>" % dwba_1)
print ("<br><span style=\"white-space: pre-line\"> %s </span>" % dwba_3)
print ("<br><span style=\"white-space: pre-line\"> %s </span>" % dwba_4)
print ("<br><span style=\"white-space: pre-line\"> %s </span>" % result)
print ("<br><span style=\"white-space: pre-line\"> %s </span>" % haha)
print ("</body>")
print ("</html>")