254 lines
8.5 KiB
Python
254 lines
8.5 KiB
Python
#!/usr/bin/python3
|
|
|
|
import numpy as np
|
|
|
|
def read_DWBA(fileName):
|
|
dataX = [] # List for the first column
|
|
data = [] # 2D list for other columns
|
|
headers = [] # List to store headers
|
|
|
|
with open(fileName, 'r') as file:
|
|
header_found = False # Flag to indicate if the header has been found
|
|
for line in file:
|
|
# Skip lines that start with '#' and empty lines
|
|
if line.startswith('#') or not line.strip():
|
|
continue
|
|
|
|
if not header_found:
|
|
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
|
|
|
|
# Split the line by whitespace
|
|
parts = line.split()
|
|
if len(parts) > 0: # Make sure there is at least one column
|
|
dataX.append(float(parts[0])) # First column
|
|
# Append the rest of the columns to data
|
|
if len(data) == 0:
|
|
# Initialize the data array with the right number of sublists
|
|
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
|
|
|
|
# print(self.headers)
|
|
return headers, dataX, data
|
|
|
|
|
|
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
|
|
|
|
print("============ Python Extraction Code ==============")
|
|
|
|
# --- 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
|
|
reaction_flag = 0
|
|
pre_find_for_elastic = False
|
|
continue
|
|
|
|
# ----- start extracting Xsec
|
|
if start_extract:
|
|
# print(line)
|
|
# print(len(line))
|
|
|
|
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
|
|
|
|
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}")
|
|
|
|
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}")
|
|
|
|
if num_angle != num_xsec and num_xsec > 0.:
|
|
if not angle_filled:
|
|
angle.append(num_angle)
|
|
data_xsec.append(num_xsec)
|
|
|
|
|
|
# ================================== 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 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(f"{'Angle':>8}")
|
|
for t in title:
|
|
file_out.write(f"{t:>40}")
|
|
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):
|
|
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")
|
|
|
|
except Exception as e:
|
|
print(f"An error occurred: {e}")
|
|
|
|
def is_float(value):
|
|
try:
|
|
float(value)
|
|
return True
|
|
except ValueError:
|
|
return False
|