working on ExtracXsecPy.py
This commit is contained in:
Ryan@Home
parent
67a5e1e459
commit
6c35437ffe
207
PyGUIQt6/ExtractXsecPy.py
Normal file
207
PyGUIQt6/ExtractXsecPy.py
Normal file
|
|
@ -0,0 +1,207 @@
|
|||
#!/usr/bin/python3
|
||||
|
||||
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
|
||||
|
||||
# --- 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:])
|
||||
pos = title[-1].find("(")
|
||||
ex = title[-1][3:pos]
|
||||
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
|
||||
|
||||
# ----- start extracting Xsec
|
||||
if start_extract:
|
||||
if len(line) < 20:
|
||||
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: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
|
||||
|
||||
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
|
||||
|
||||
if num_angle != num_xsec and num_xsec > 0.:
|
||||
if not angle_filled:
|
||||
angle.append(num_angle)
|
||||
data_xsec.append(num_xsec)
|
||||
|
||||
# ------ 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:
|
||||
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
|
||||
|
||||
# ================================== 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}")
|
||||
|
||||
def is_float(value):
|
||||
try:
|
||||
float(value)
|
||||
return True
|
||||
except ValueError:
|
||||
return False
|
||||
Loading…
Reference in New Issue
Block a user