rtang/PtolemyGUI
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

use python ExtractXSecPy for extracting xsec. now is cern root independent

Browse Source
This commit is contained in:
Ryan Tang 2024-11-04 14:14:32 -05:00
parent 3525cd5b0e
commit 8a9a77e857
3 changed files with 232 additions and 211 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
PyGUIQt6/DWBA
View File

@ -47,5 +47,5 @@
#82Kr(d,p)83Kr 0 1d3/2 3/2+ 1.220 10MeV/u AK
#82Kr(d,p)83Kr 0 2s1/2 1/2+ 1.220 10MeV/u AK
10Be(d,d)10Be 0 none 0+ 0.000 9MeV/u ZZ
10Be(d,d)10Be 0 none 0+ 0.000 9MeV/u AA
11C(d,p)12C 3/2- 0p1/2 2+ 4.4 10MeV/u AK
11C(d,p)12C 3/2- 0p1/2 1+ 15.1 10MeV/u AK

386
PyGUIQt6/ExtractXsecPy.py
View File

@ -2,206 +2,228 @@
import numpy as np
def extract_xsec(read_file, index_for_elastic=1):
# index_for_elastic = 1 ; for Ratio
# index_for_elastic = 2 ; for Total
# index_for_elastic = 3 ; for Rutherford
def extract_xsec(read_file, index_for_elastic):
# index_for_elastic = 1 ; for Ratio
# index_for_elastic = 2 ; for Total
# index_for_elastic = 3 ; for Rutherford
# --- open file
try:
with open(read_file, 'r') as file_in:
# ---- variable for Xsec
title = []
data_matrix = []
angle = []
Ex = []
total_xsec = []
reaction = []
angle_min = 0
angle_max = 0
angle_step = -1
print("============ Python Extraction Code ==============")
# ================================== read file
line_num = 0
data_xsec = []
# --- open file
try:
with open(read_file, 'r') as file_in:
# ---- variable for Xsec
title = []
data_matrix = []
angle = []
Ex = []
total_xsec = []
reaction = []
angle_min = 0
angle_max = 0
angle_step = -1
# ================================== read file
line_num = 0
data_xsec = []
start_extract = False
angle_filled = False
num_cal = 0
reaction_flag = 0
pre_find_for_elastic = False
print("======================================================")
for line in file_in:
line_num += 1
line = line.strip()
# ---- find Title
find_str = "$============================================"
pos = line.find(find_str)
if pos != -1:
title.append(line[pos + len(find_str) + 1:])
pos1 = title[-1].find("=")
pos2 = title[-1].find("(")
ex = title[-1][pos1+1:pos2]
Ex.append(float(ex))
data_xsec.clear()
num_cal += 1
continue
# ---- find Reaction
find_str = "0INPUT... CHANNEL"
pos = line.find(find_str)
if pos != -1:
reaction.append(line[pos + len(find_str) + 1:])
reaction_flag = 1 # 1 for (d,d) or (p,p)
continue
find_str = "REACTION:"
pos = line.find(find_str)
if pos != -1:
reaction.append(line[pos + len(find_str) + 1:])
reaction_flag = 2 # 2 for (d,p) or (p,d)
continue
# ----- pre find for Elastic scattering
find_str = "0 OPTICAL MODEL SCATTERING FOR THE OUTGOING CHANNEL"
pos = line.find(find_str)
if pos != -1:
pre_find_for_elastic = True
continue
# ----- find angle setting when not known
if angle_step == -1:
find_str = "anglemin="
pos = line.find(find_str)
if pos != -1:
pos1 = line.find("anglemax=")
angle_min = float(line[pos + len(find_str):pos1].strip())
pos = pos1
pos1 = line.find("anglestep=")
angle_max = float(line[pos + len(find_str):pos1].strip())
angle_step = float(line[pos1 + len(find_str) + 1:].strip())
continue
# ----- check if start extracting Xsec or not
find_str = "dumpdumpdump"
if reaction_flag == 1 and pre_find_for_elastic:
find_str = "C.M. LAB RUTHERFORD"
if reaction_flag == 2:
find_str = "0 C.M. REACTION REACTION LOW L HIGH L % FROM"
pos = line.find(find_str)
if pos != -1:
start_extract = True
continue
# ------ find total Xsec, if found stop extraction
find_str = "dumpdumpdump"
if reaction_flag == 1 and pre_find_for_elastic:
find_str = "0TOTAL REACTION CROSS SECTION ="
if reaction_flag == 2:
find_str = "0TOTAL:"
pos = line.find(find_str)
if pos != -1:
# print(line)
total_xsec.append(float(line[pos + len(find_str):-3].strip()))
print(f"{num_cal:2d} | {title[-1]} | total Xsec(4pi): {total_xsec[-1]} mb")
data_matrix.append(data_xsec[:])
angle_filled = True
start_extract = False
angle_filled = False
num_cal = 0
reaction_flag = 0
pre_find_for_elastic = False
print("======================================================")
continue
# ----- start extracting Xsec
if start_extract:
# print(line)
# print(len(line))
for line in file_in:
line_num += 1
line = line.strip()
if len(line) < 80:
continue
if line.find("C.M.") > 0 or line.find("ANGLE") > 0:
continue
if line.find("INCOMING ELASTIC") > 0 or line.find("LX") > 0 :
continue
# ---- find Title
find_str = "$============================================"
pos = line.find(find_str)
if pos != -1:
title.append(line[pos + len(find_str) + 1:])
pos = title[-1].find("(")
ex = title[-1][3:pos]
Ex.append(float(ex))
data_xsec.clear()
num_cal += 1
continue
num_angle = num_xsec = None
if reaction_flag == 1: # Elastics (d,d) or (p,p)
num_angle = float(line[0:7].strip())
if index_for_elastic == 1:
num_xsec = float(line[15:28].strip())
elif index_for_elastic == 2:
if len(line) > 60:
num_xsec = float(line[29:40].strip())
else:
num_xsec = -404
else:
if len(line) > 60:
num_xsec = float(line[57:71].strip())
else:
num_xsec = -404
# print(f"{num_angle}, {num_xsec}")
# ---- find Reaction
find_str = "0INPUT... CHANNEL"
pos = line.find(find_str)
if pos != -1:
reaction.append(line[pos + len(find_str) + 1:])
reaction_flag = 1 # 1 for (d,d) or (p,p)
continue
if reaction_flag == 2: # InElastics (d,p) or (p,d)
# print(line)
# print(line[0:6])
# print(line[6:19])
if is_float(line[0:6].strip()):
num_angle = float(line[0:6].strip())
num_xsec = float(line[6:19].strip())
else:
num_angle = -1.0
num_xsec = -1.0
# print(f"{num_angle}, {num_xsec}")
find_str = "REACTION:"
pos = line.find(find_str)
if pos != -1:
reaction.append(line[pos + len(find_str) + 1:])
reaction_flag = 2 # 2 for (d,p) or (p,d)
continue
if num_angle != num_xsec and num_xsec > 0.:
if not angle_filled:
angle.append(num_angle)
data_xsec.append(num_xsec)
# ----- pre find for Elastic scattering
find_str = "0 OPTICAL MODEL SCATTERING FOR THE OUTGOING CHANNEL"
pos = line.find(find_str)
if pos != -1:
pre_find_for_elastic = True
continue
# ----- find angle setting when not known
if angle_step == -1:
find_str = "anglemin="
pos = line.find(find_str)
if pos != -1:
pos1 = line.find("anglemax=")
angle_min = float(line[pos + len(find_str):pos1].strip())
pos = pos1
pos1 = line.find("anglestep=")
angle_max = float(line[pos + len(find_str):pos1].strip())
angle_step = float(line[pos1 + len(find_str) + 1:].strip())
continue
# ================================== summary
print("====================== Total number of lines read : ", line_num)
print(f"Angle : {angle_min:5.2f}, {angle_max:5.2f} | step : {angle_step:5.2f}")
print("Number of Angles read : ", len(angle))
print("Number of data read : ", len(data_xsec))
print("Number of Reactions : ", len(reaction))
# ----- check if start extracting Xsec or not
find_str = "dumpdumpdump"
if reaction_flag == 1 and pre_find_for_elastic:
find_str = "C.M. LAB RUTHERFORD"
if reaction_flag == 2:
find_str = "0 C.M. REACTION REACTION LOW L HIGH L % FROM"
pos = line.find(find_str)
if pos != -1:
start_extract = True
continue
print("----------------------------- list of Calculation")
reaction_str_len = max(len(r) for r in reaction)
# ----- start extracting Xsec
if start_extract:
if len(line) < 20:
continue
partial_xsec = []
for i in range(num_cal):
partial_sum_xsec = 0.0
for j in range(len(data_matrix[i])):
theta = np.radians(angle[j])
d_theta = np.radians(angle_step)
phi = 2 * np.pi
partial_sum_xsec += data_matrix[i][j] * np.sin(theta) * d_theta * phi
partial_xsec.append(partial_sum_xsec)
num_angle = num_xsec = None
if reaction_flag == 1: # Elastics (d,d) or (p,p)
num_angle = float(line[0:7].strip())
if index_for_elastic == 1:
num_xsec = float(line[15:30].strip())
elif index_for_elastic == 2:
if len(line) > 60:
num_xsec = float(line[30:43].strip())
else:
num_xsec = -404
else:
if len(line) > 60:
num_xsec = float(line[57:71].strip())
else:
num_xsec = -404
pos = title[i].find(")")
print(f"{reaction[i]:>{reaction_str_len + 3}} | {title[i][pos + 1:]} | Xsec({angle_min:3.0f}-{angle_max:3.0f} deg) : {partial_sum_xsec:.6f} mb")
print("---------------------------------------------------")
if reaction_flag == 2: # InElastics (d,p) or (p,d)
if is_float(line[0:7].strip()):
num_angle = float(line[0:7].strip())
num_xsec = float(line[7:26].strip())
else:
num_angle = -1.0
num_xsec = -1.0
# ================================== save *.Ex.txt
save_ex_name = f"{read_file[:-4]}.Ex.txt"
print("Output : ", save_ex_name)
with open(save_ex_name, 'w') as file_ex:
file_ex.write("//generated_by_ExtractXSec.h____Ex____Xsec(4pi)____SF____sigma\n")
for i in range(num_cal):
file_ex.write(f"{Ex[i]:9.5f} {partial_xsec[i]:9.5f} 1.0 0.000\n")
file_ex.write("#=====End_of_File\n")
if num_angle != num_xsec and num_xsec > 0.:
if not angle_filled:
angle.append(num_angle)
data_xsec.append(num_xsec)
# ================================== save file.Xsec.txt
save_file_name = f"{read_file[:-4]}.Xsec.txt"
print("Output : ", save_file_name)
with open(save_file_name, 'w') as file_out:
for r in reaction:
file_out.write(f"#{r}\n")
# ------ find total Xsec, if found stop extraction
find_str = "dumpdumpdump"
if reaction_flag == 1 and pre_find_for_elastic:
find_str = "0TOTAL REACTION CROSS SECTION ="
if reaction_flag == 2:
find_str = "0TOTAL:"
pos = line.find(find_str)
file_out.write(f"{"Angle":>8}")
for t in title:
file_out.write(f"{t:>40}")
file_out.write("\n")
if pos != -1:
total_xsec.append(float(line[pos + len(find_str):].strip()))
print(f"{num_cal:2d} | {title[-1]} | total Xsec(4pi): {total_xsec[-1]} mb")
data_matrix.append(data_xsec[:])
angle_filled = True
start_extract = False
reaction_flag = 0
pre_find_for_elastic = False
continue
for i in range(len(angle)):
file_out.write(f"{angle[i]:8.3f}\t")
for j in range(num_cal):
if data_matrix[j][i] > 1e8 :
file_out.write(f"{data_matrix[j][i]:>40.6e}")
else:
file_out.write(f"{data_matrix[j][i]:>40.6f}")
file_out.write("\n")
# ================================== summary
print("====================== Total number of lines read : ", line_num)
print(f"Angle : {angle_min:5.2f}, {angle_max:5.2f} | step : {angle_step:5.2f}")
print("Number of Angles read : ", len(angle))
print("Number of data read : ", len(data_xsec))
print("Number of Reactions : ", len(reaction))
print("----------------------------- list of Calculation")
reaction_str_len = max(len(r) for r in reaction)
partial_xsec = []
for i in range(num_cal):
partial_sum_xsec = 0.0
for j in range(len(data_matrix[i])):
theta = np.radians(angle[j])
d_theta = np.radians(angle_step)
phi = 2 * np.pi
partial_sum_xsec += data_matrix[i][j] * np.sin(theta) * d_theta * phi
partial_xsec.append(partial_sum_xsec)
pos = title[i].find(")")
print(f"{reaction[i]:>{reaction_str_len + 3}} | {title[i][pos + 1:]} | Xsec({angle_min:3.0f}-{angle_max:3.0f} deg) : {partial_sum_xsec:.6f} mb")
print("---------------------------------------------------")
# ================================== save *.Ex.txt
save_ex_name = f"{read_file[:-4]}.Ex.txt"
print("Output : ", save_ex_name)
with open(save_ex_name, 'w') as file_ex:
file_ex.write("//generated_by_ExtractXSec.h____Ex____Xsec(4pi)____SF____sigma\n")
for i in range(num_cal):
file_ex.write(f"{Ex[i]:9.5f} {partial_xsec[i]:9.5f} 1.0 0.000\n")
file_ex.write("#=====End_of_File\n")
# ================================== save file.Xsec.txt
save_file_name = f"{read_file[:-4]}.Xsec.txt"
print("Output : ", save_file_name)
with open(save_file_name, 'w') as file_out:
for r in reaction:
file_out.write(f"#{r}\n")
file_out.write("Angel\t")
for t in title:
file_out.write(f"{t:19}")
file_out.write("\n")
for i in range(len(angle)):
file_out.write(f"{angle[i]:8.3f}\t")
for j in range(num_cal):
file_out.write(f"{data_matrix[j][i]:19.6f}")
file_out.write("\n")
except Exception as e:
print(f"An error occurred: {e}")
except Exception as e:
print(f"An error occurred: {e}")
def is_float(value):
try:
float(value)
return True
except ValueError:
return False
try:
float(value)
return True
except ValueError:
return False

53
PyGUIQt6/PtolemyGUIPy.py
View File

@ -17,6 +17,8 @@ from PyQt6.QtWebEngineWidgets import QWebEngineView
import plotly.graph_objects as go
import tempfile
from ExtractXsecPy import extract_xsec
class PlotWindow(QWidget):
def __init__(self, XsecFile):
super().__init__()
@ -27,7 +29,7 @@ class PlotWindow(QWidget):
self.x = []
self.data = []
self.headers = []
self.x, self.data, self.headers = self.read_data(XsecFile)
self.read_data(XsecFile)
self.log_scale_checkbox = QCheckBox("Use Log Scale for Y-Axis")
self.log_scale_checkbox.setChecked(True)
@ -50,9 +52,9 @@ class PlotWindow(QWidget):
self.plot_plotly_graph()
def read_data(self,file_path):
x = [] # List for the first column
data = [] # 2D list for other columns
headers = [] # List to store headers
self.x = [] # List for the first column
self.data = [] # 2D list for other columns
self.headers = [] # List to store headers
with open(file_path, 'r') as file:
header_found = False # Flag to indicate if the header has been found
@ -62,12 +64,7 @@ class PlotWindow(QWidget):
continue
if not header_found:
parts = line.split('ELab')
elab_parts = [parts[0]] # Start with the first part
for part in parts[1:]:
elab_parts.append('ELab' + part) # Prepend 'ELab' to each subsequent part
headers = elab_parts # Use the split parts as headers
self.headers = line.split() # Use the split parts as headers
header_found = True # Set the flag to True to skip this block in future iterations
# print(f"ELab parts found: {elab_parts}") # Print or process this as needed
continue
@ -75,15 +72,15 @@ class PlotWindow(QWidget):
# Split the line by whitespace
parts = line.split()
if len(parts) > 0: # Make sure there is at least one column
x.append(float(parts[0])) # First column
self.x.append(float(parts[0])) # First column
# Append the rest of the columns to data
if len(data) == 0:
if len(self.data) == 0:
# Initialize the data array with the right number of sublists
data = [[] for _ in range(len(parts) - 1)]
self.data = [[] for _ in range(len(parts) - 1)]
for i in range(len(parts) - 1):
data[i].append(float(parts[i + 1])) # Rest of the columns
return x, data, headers
self.data[i].append(float(parts[i + 1])) # Rest of the columns
# print(self.headers)
def plot_plotly_graph(self):
# Create a Plotly figure
@ -142,15 +139,15 @@ class PlotWindow(QWidget):
x1=1.0,
y1=1.0,
line=dict(
color="black",
width=1,
color="black",
width=1,
)
)
# Save the plot as an HTML file in a temporary location
with tempfile.NamedTemporaryFile(delete=False, suffix=".html") as tmp_file:
fig.write_html(tmp_file.name)
html_file = tmp_file.name
fig.write_html(tmp_file.name)
html_file = tmp_file.name
# Load the HTML file in QWebEngineView
self.web_view.setUrl(QUrl.fromLocalFile(html_file))
@ -255,7 +252,7 @@ class MyWindow(QMainWindow):
self.text_edit = QTextEdit()
self.text_edit.setLineWrapMode(QTextEdit.LineWrapMode.NoWrap)
font = QFont("Courier New", 8) # You can adjust the size as needed
font = QFont("Courier New", 10) # You can adjust the size as needed
self.text_edit.setFont(font)
self.leStatus = QLineEdit("")
@ -334,6 +331,8 @@ class MyWindow(QMainWindow):
def CalDWBA(self):
self.SaveFile()
self.BashCommand("cd ../Cleopatra; make;cd ../PyGUIQt6")
print(" Is Create InFile : " + str(self.chkCreateInFile.isChecked() ))
@ -348,7 +347,7 @@ class MyWindow(QMainWindow):
self.BashCommand("../Cleopatra/InFileCreator " + self.DWBAFileName + aMin + aMax + aSize)
isRunOK = False
isRunOK = True
if self.chkRunPtolemy.isChecked() :
os_name = platform.system()
@ -358,14 +357,14 @@ class MyWindow(QMainWindow):
if os_name == "Darwin":
self.BashCommand("../Cleopatra/ptolemy_mac <" + self.DWBAFileName + ".in>" + " " + self.DWBAFileName + ".out")
if self.bashResult.returncode == 0 :
isRunOK = True
else:
if self.bashResult.returncode != 0 :
isRunOK = False
self.leStatus.setText("Ptolemy Run Error. Should check the out File.")
if isRunOK and self.chkExtracrXsec.isChecked() and self.file_exists(self.DWBAFileName + ".out") :
option = str(self.cbXsec.currentIndex())
self.BashCommand("../Cleopatra/ExtractXSec " + self.DWBAFileName + ".out " + option)
extract_xsec(self.DWBAFileName + ".out", self.cbXsec.currentIndex())
# option = str(self.cbXsec.currentIndex())
# self.BashCommand("../Cleopatra/ExtractXSec " + self.DWBAFileName + ".out " + option)
if self.chkPlot.isChecked() and self.file_exists(self.DWBAFileName + ".Xsec.txt") :
if self.A_file_changed_after_B_file(self.DWBAFileName + ".Xsec.txt", self.DWBAFileName + ".out") :