SOLARIS_Web_Simulation/simpleInput.py

#!/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

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 :

  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 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>")