From 8a9a77e8575e74a3f9b9de706eacd866a9cf63a0 Mon Sep 17 00:00:00 2001 From: "Ryan@SOLARIS_testStation" Date: Mon, 4 Nov 2024 14:14:32 -0500 Subject: [PATCH] use python ExtractXSecPy for extracting xsec. now is cern root independent --- PyGUIQt6/DWBA | 4 +- PyGUIQt6/ExtractXsecPy.py | 386 ++++++++++++++++++++------------------ PyGUIQt6/PtolemyGUIPy.py | 53 +++--- 3 files changed, 232 insertions(+), 211 deletions(-) diff --git a/PyGUIQt6/DWBA b/PyGUIQt6/DWBA index 16da375..5d0c236 100644 --- a/PyGUIQt6/DWBA +++ b/PyGUIQt6/DWBA @@ -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 \ No newline at end of file diff --git a/PyGUIQt6/ExtractXsecPy.py b/PyGUIQt6/ExtractXsecPy.py index dee5e2c..ae0a05a 100644 --- a/PyGUIQt6/ExtractXsecPy.py +++ b/PyGUIQt6/ExtractXsecPy.py @@ -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 diff --git a/PyGUIQt6/PtolemyGUIPy.py b/PyGUIQt6/PtolemyGUIPy.py index 062295f..a6aae6a 100755 --- a/PyGUIQt6/PtolemyGUIPy.py +++ b/PyGUIQt6/PtolemyGUIPy.py @@ -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") :