Development

requirements.txt

@@ -1,5 +1,8 @@
 sphinx
 sphinx-rtd-theme
-numpy
-pandas
-seaborn
+numpy~=1.26.4
+pandas~=2.2.1
+seaborn~=0.13.2
+
+matplotlib~=3.8.4
+scipy~=1.13.1

src/Utilities/fix_no_dose_error_composite.py

@@ -0,0 +1,22 @@
+import sys
+import os
+
+# Add the src directory to the sys.path
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+from src import io_snc
+
+def main():
+    """
+    main function to fix measured files
+
+    Returns
+    -------
+
+    """
+    # ask user for file path of measured file
+    file_path = input("Please enter the file path: ")
+    return file_path
+
+
+if __name__ == "__main__":
+    main()

src/correction_coefficients/correction_coefficients.py

@@ -4,7 +4,7 @@
 The module includes the following functions:
 
 - `get_pr_data_from_acm_files`: Extracts the average counts/50ms from ACM files in a specified directory.
-- `exponential_fit`: Defines an exponential fitting function for curve fitting.
+- `saturation_fit`: Defines an exponential fitting function for curve fitting.
 - `get_correction_coefficients`: Calculates the correction coefficients based on the average counts/50ms and relative signal values.
 - `plot_correction_curve`: Plots the correction curve with the data points and the fitted exponential function.
 - `plot_counts_per_50ms`: Plots the counts per 50ms at different nominal dose rates.
@@ -17,10 +17,14 @@
 import os
 import sys
 
-import numpy as np
 import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
 import seaborn as sns
 from scipy.optimize import curve_fit
+import datetime
+
+
 
 # Add the src directory to the sys.path
 sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
@@ -42,36 +46,49 @@ def get_pr_data_from_acm_files(pr_measurement_folder):
         List of tuples containing the file name and the average counts/50ms for each nominal dose rate.
     """
     diode_numbers = np.array([758, 759, 760, 761, 692, 693, 694, 695, 626, 627, 628, 629]).astype(str)
+
+    # initialize a dictionary to store the average counts/50ms  for each file over all
+    # diodes in the middle 200 frames.
+    counts_per_50ms_dict = {}
+
+    # initialize a list to store the average counts per 50ms for each file
     counts_per_50ms_list = []
 
     for acm_file in os.listdir(pr_measurement_folder):
         acm_file_path = os.path.join(pr_measurement_folder, acm_file)
         frame_data_df, diode_data_df, bkrnd_and_calibration_df = io_snc.parse_acm_file(acm_file_path)
 
-        # Ensure the frame data is sorted by time or frame number if necessary
-        diode_data_df.sort_index(inplace=True)
-
-        # Calculate the count rate by differencing the accumulated counts
-        diode_data_diff = diode_data_df.diff().iloc[1:]
-
         # Grab the count rate from the middle 100 frames of the measurement
-        num_frames = len(diode_data_diff)
+        num_frames = len(diode_data_df)
         start_frame = num_frames // 2 - 50
         end_frame = start_frame + 100
-        middle_frames = diode_data_diff.iloc[start_frame:end_frame]
+        middle_frames = diode_data_df.iloc[start_frame:end_frame]
+
+        # Calculate the count rate by differencing the accumulated counts
+        diode_data_diff_df = middle_frames.diff().iloc[1:]
 
         # Extract the count rate for the specified diodes in the middle frames
-        count_rate_df = middle_frames.loc[:, diode_numbers]
+        count_rate_df = diode_data_diff_df.loc[:, diode_numbers]
+
+        # Calculate the average count rate for each frame across the selected diodes
+        avg_count_rate_per_frame = count_rate_df.mean(axis=1)
+        avg_count_rate_per_frame = avg_count_rate_per_frame.values
+        counts_per_50ms_dict[acm_file] = avg_count_rate_per_frame
 
         # Calculate the average count rate for the selected diodes
         avg_count_rate = count_rate_df.mean().mean()
         counts_per_50ms_list.append((acm_file, avg_count_rate))
 
+    # Convert the dictionary to a DataFrame
+    counts_per_50ms_df = pd.DataFrame.from_dict(counts_per_50ms_dict, orient='index')
+    # Save the DataFrame to a CSV file
+    counts_per_50ms_df.to_csv('counts_per_50ms.csv')
+
     counts_per_50ms_list = sorted(counts_per_50ms_list, key=lambda x: x[1])
     return counts_per_50ms_list
 
 
-def exponential_fit(x, a, b, c):
+def saturation_func(x, a, b, c):
     """
     Exponential fitting function.
 
@@ -91,7 +108,7 @@ def exponential_fit(x, a, b, c):
     array_like
         The calculated dependent variable values.
     """
-    return a * np.exp(b * x) + c
+    return c - a * np.exp(-b * x)
 
 
 def get_correction_coefficients(counts_per_50ms, relative_signal):
@@ -111,10 +128,49 @@ def get_correction_coefficients(counts_per_50ms, relative_signal):
         The fitted coefficients (a, b, c) for the exponential correction function.
     """
     counts_per_50ms = np.array([data[1] for data in counts_per_50ms])
+    initial_guess = [0.035, 5.21e-5, 1.0]
+    # noinspection PyTupleAssignmentBalance
+    params, covariance = curve_fit(saturation_func, counts_per_50ms,
+                                   relative_signal, p0=initial_guess)
+    a_fit, b_fit, c_fit = params
 
-    # Perform the curve fitting
-    popt, _ = curve_fit(exponential_fit, counts_per_50ms, relative_signal)
-    return tuple(popt)
+    return a_fit, b_fit, c_fit, covariance
+
+def get_coefficient_file_path(date=None):
+    """
+    Get the path to the file where the correction coefficients will be written.
+
+    Returns
+    -------
+    str
+        The file path for the correction coefficients.
+    """
+    # get user input on type of whether it DPP or PR
+    coefficient_file_directory = input("Enter the directory to the file where the correction coefficients will be written: ")
+
+    # Get the DRD type, ArcCheck SN from the user
+    DRD_type = input("Enter the DRD type: ")
+    ArcCheck_SN = input("Enter the ArcCheck SN: ")
+
+    # Get today's date
+    today = datetime.date.today()
+    date = today.strftime("%Y-%m-%d")
+
+    # construct coefficient file name from DRD type, ArcCheck SN and Date
+    coefficient_file_name = f"{DRD_type}_{ArcCheck_SN}_{date}.txt"
+
+    # construct the full path to the coefficient file
+    coefficient_file_path = os.path.join(coefficient_file_directory, coefficient_file_name)
+
+    return coefficient_file_path
+
+def write_correction_coefficients_to_file(coefficient_file_path, a_fit, b_fit, c_fit, covariance):
+    with open(coefficient_file_path, 'w') as f:
+        # Write the variable names and their values to the file
+        f.write(f"a_fit = {a_fit}\n")
+        f.write(f"b_fit = {b_fit}\n")
+        f.write(f"c_fit = {c_fit}\n")
+        f.write(f"covariance = {covariance}\n")
 
 
 def plot_correction_curve(counts_per_50ms, relative_signal, a_fit, b_fit, c_fit):
@@ -138,7 +194,7 @@ def plot_correction_curve(counts_per_50ms, relative_signal, a_fit, b_fit, c_fit)
 
     # Generate fitted curve data points
     x_fit = np.linspace(counts_per_50ms.min(), counts_per_50ms.max(), 500)
-    y_fit = exponential_fit(x_fit, a_fit, b_fit, c_fit)
+    y_fit = saturation_func(x_fit, a_fit, b_fit, c_fit)
 
     # Set Seaborn style
     sns.set(style="whitegrid")
@@ -148,7 +204,7 @@ def plot_correction_curve(counts_per_50ms, relative_signal, a_fit, b_fit, c_fit)
     sns.lineplot(x=x_fit, y=y_fit, label=f'Exponential fit, a={a_fit:.4f}, b={b_fit:.2e}, c={c_fit:.3f}')
     sns.scatterplot(x=counts_per_50ms, y=relative_signal, label='Data points')
 
-    plt.ioff()  # Turn off interactive mode
+    #plt.ioff()  # Turn off interactive mode
 
     # Plot the data points and the fitted curve
     plt.figure(figsize=(8, 6))
@@ -162,8 +218,8 @@ def plot_correction_curve(counts_per_50ms, relative_signal, a_fit, b_fit, c_fit)
     plt.title('Pulse rate (MU/min)')
     plt.legend()
     plt.grid(True)
-
     plt.show()
+    plt.savefig('counts_per_50ms.png')
 
 
 def plot_counts_per_50ms(counts_per_50ms):
@@ -179,7 +235,7 @@ def plot_counts_per_50ms(counts_per_50ms):
 
     # Set Seaborn style
     sns.set(style="whitegrid")
-
+    #plt.ioff()  # Turn off interactive mode
     # Create a plot with Seaborn
     plt.figure(figsize=(10, 6))
     sns.scatterplot(x=range(len(counts_per_50ms)), y=counts_per_50ms, label='Counts/50ms')
@@ -190,20 +246,23 @@ def plot_counts_per_50ms(counts_per_50ms):
     plt.title('Counts per 50ms at different nominal dose rates')
     plt.legend()
     plt.grid(True)
-
     plt.show()
+    plt.savefig('counts_per_50ms_nominal.png')
 
 
 def main():
     """
     Main function to execute the correction coefficient calculations and plotting.
     """
+
     pr_measurement_folder = r"P:\02_QA Equipment\02_ArcCheck\05_Commissoning\03_NROAC\Dose Rate Dependence Fix\NRO PR Coefficient"
     counts_per_50ms = get_pr_data_from_acm_files(pr_measurement_folder)
-    plot_counts_per_50ms(counts_per_50ms)
+    # plot_counts_per_50ms(counts_per_50ms)
     relative_signal = np.array([0.964, 0.971, 0.980, 0.988, 0.994, 0.995, 1.000, 1.000])
-    a_fit, b_fit, c_fit = get_correction_coefficients(counts_per_50ms, relative_signal)
-    plot_correction_curve(counts_per_50ms, relative_signal, a_fit, b_fit, c_fit)
+    a_fit, b_fit, c_fit, covariance = get_correction_coefficients(counts_per_50ms, relative_signal)
+    coefficient_file_path = get_coefficient_file_path()
+    write_correction_coefficients_to_file(coefficient_file_path, a_fit, b_fit, c_fit, covariance)
+    # plot_correction_curve(counts_per_50ms, relative_signal, a_fit, b_fit, c_fit)
 
 
 if __name__ == "__main__":

src/corrections.py

@@ -73,7 +73,9 @@ def get_intrinsic_corrections(array_data):
 
 def apply_jager_corrections(counts_accumulated_df, bkrnd_and_calibration_df, intrinsic_corrections=None):
     """Apply Jager pulse rate and dose per pulse corrections."""
-    a_pr, b_pr, c_pr = 0.035, 5.21 * 10 ** -5, 1
+
+    # a_pr, b_pr, c_pr = 0.035, 5.21 * 10 ** -5, 1 # Jager coefficients from paper
+    a_pr, b_pr, c_pr = 0.03921672747199777, 6.499040838838414e-055, 0.999142799734291 # Coefficients from NRO Measurements
     jager_pr_coefficients = np.array([a_pr, b_pr, c_pr])
 
     a_dpp, b_dpp, c_dpp = 0.0978, 3.33 * 10 ** -5, 1.011

src/main.py

@@ -199,3 +199,4 @@ def main():
 if __name__ == "__main__":
     main()
 
+