Pr 90

libra_toolbox/neutron_detection/activation_foils/calibration.py

@@ -1,5 +1,5 @@
-from dataclasses import dataclass
-from typing import List
+from dataclasses import dataclass, field
+from typing import List, Dict, Optional
 import datetime
 import numpy as np
 
@@ -26,24 +26,44 @@ class Nuclide:
     """
 
     name: str
-    energy: List[float] = None
-    intensity: List[float] = None
-    half_life: float = None
-    atomic_mass: float = None
+    energy: Optional[List[float]] = None
+    intensity: Optional[List[float]] = None
+    half_life: Optional[float] = None
+    atomic_mass: Optional[float] = None
     abundance: float = 1.00
+    _uncalibrated_measured_energies: Optional[Dict] = field(default=None, repr=False)
 
     @property
     def decay_constant(self):
         """
         Returns the decay constant of the nuclide in 1/s.
         """
         return np.log(2) / self.half_life
+
+    def calibrated_measured_energies(self, channel_nb, calibration_coeffs):
+        """
+        Returns the calibrated measured energies of the nuclide in keV.
+        """
+        if self._uncalibrated_measured_energies is None:
+            return None
+        else:
+            uncalibrated = np.array(
+                self._uncalibrated_measured_energies.get(channel_nb, []),
+                dtype=float
+            )
+            return np.polyval(calibration_coeffs, uncalibrated)
 
 
+# ba133 = Nuclide(
+#     name="Ba133",
+#     energy=[80.9979, 276.3989, 302.8508, 356.0129, 383.8485],
+#     intensity=[0.329, 0.0716, 0.1834, 0.6205, 0.0894],
+#     half_life=10.551 * 365.25 * 24 * 3600,
+# )
 ba133 = Nuclide(
     name="Ba133",
-    energy=[80.9979, 276.3989, 302.8508, 356.0129, 383.8485],
-    intensity=[0.329, 0.0716, 0.1834, 0.6205, 0.0894],
+    energy=[276.3989, 302.8508, 356.0129, 383.8485],
+    intensity=[0.0716, 0.1834, 0.6205, 0.0894],
     half_life=10.551 * 365.25 * 24 * 3600,
 )
 co60 = Nuclide(

libra_toolbox/neutron_detection/activation_foils/compass.py

@@ -169,6 +169,7 @@ class Measurement:
     stop_time: Union[datetime.datetime, None]
     name: str
     detectors: List[Detector]
+    detector_type: str = "NaI"  # Default detector type, can be 'NaI' or 'HPGe'
 
     def __init__(self, name: str) -> None:
         """
@@ -415,6 +416,8 @@ def compute_detection_efficiency(
         calibration_coeffs: np.ndarray,
         channel_nb: int,
         search_width: float = 800,
+        threshold_overlap: float = 200,
+        summing_method: str = 'sum_gaussian',
     ) -> Union[np.ndarray, float]:
         """
         Computes the detection efficiency of a check source given the
@@ -433,6 +436,12 @@ def compute_detection_efficiency(
             calibration_coeffs: the calibration polynomial coefficients for the detector
             channel_nb: the channel number of the detector
             search_width: the search width for the peak fitting
+            threshold_overlap: the threshold width for considering two peaks as overlapping
+            summing_method: method to sum counts under the peak, either 'sum_gaussian' or 'sum_histogram'
+                            with 'sum_gaussian' fitting a Gaussian to the peak and integrating it, 
+                            and with 'sum_histogram' summing the histogram counts under the peak.
+                            'sum_histogram' SHOULD NOT BE USED for OVERLAPPING peaks as it will
+                            overestimate the counts.
 
         Returns:
             the detection efficiency
@@ -448,11 +457,18 @@ def compute_detection_efficiency(
 
         calibrated_bin_edges = np.polyval(calibration_coeffs, bin_edges)
 
+        peak_energies = self.check_source.nuclide.calibrated_measured_energies(channel_nb, calibration_coeffs)
+        if peak_energies is None:
+            print("TOOLBOX: No calibrated measured energies found for the check source. Cannot compute detection efficiency.")
+            peak_energies = self.check_source.nuclide.energy
+
         nb_counts_measured = get_multipeak_area(
             hist,
             calibrated_bin_edges,
-            self.check_source.nuclide.energy,
+            peak_energies,
             search_width=search_width,
+            threshold_overlap=threshold_overlap,
+            summing_method=summing_method,
         )
 
         nb_counts_measured = np.array(nb_counts_measured)
@@ -497,26 +513,58 @@ def get_peaks(self, hist: np.ndarray, **kwargs) -> np.ndarray:
             the peak indices in ``hist``
         """
 
-        # peak finding parameters
-        start_index = 100
-        prominence = 0.10 * np.max(hist[start_index:])
-        height = 0.10 * np.max(hist[start_index:])
-        width = [10, 150]
-        distance = 30
-        if self.check_source.nuclide == na22:
+        if self.detector_type.lower() == 'nai':
+            # peak finding parameters
             start_index = 100
-            height = 0.1 * np.max(hist[start_index:])
-            prominence = 0.1 * np.max(hist[start_index:])
+            prominence = 0.10 * np.max(hist[start_index:])
+            height = 0.10 * np.max(hist[start_index:])
             width = [10, 150]
             distance = 30
-        elif self.check_source.nuclide == co60:
-            start_index = 400
-            height = 0.60 * np.max(hist[start_index:])
-            prominence = None
-        elif self.check_source.nuclide == ba133:
-            width = [10, 200]
-        elif self.check_source.nuclide == mn54:
-            height = 0.6 * np.max(hist[start_index:])
+            if self.check_source.nuclide == na22:
+                start_index = 100
+                height = 0.4 * np.max(hist[start_index:])
+                prominence = 0.4 * np.max(hist[start_index:])
+                width = [10, 150]
+                distance = 30
+            elif self.check_source.nuclide == co60:
+                start_index = 400
+                height = 0.60 * np.max(hist[start_index:])
+                prominence = None
+            elif self.check_source.nuclide == ba133:
+                start_index = 150
+                height = 0.10 * np.max(hist[start_index:])
+                prominence = 0.10 * np.max(hist[start_index:])
+            elif self.check_source.nuclide == mn54:
+                height = 0.6 * np.max(hist[start_index:])
+        elif self.detector_type.lower() == 'hpge':
+            # peak finding parameters for HPGe detectors
+            start_index = 10
+            prominence = 0.50 * np.max(hist[start_index:])
+            height = 0.50 * np.max(hist[start_index:])
+            width = [2, 50]
+            distance = 100
+            if self.check_source.nuclide == na22:
+                start_index = 100
+                height = 0.4 * np.max(hist[start_index:])
+                prominence = 0.4 * np.max(hist[start_index:])
+                distance = 100
+            elif self.check_source.nuclide == co60:
+                height = 0.5 * np.max(hist[start_index:])
+                prominence = 0.5 * np.max(hist[start_index:])
+            elif self.check_source.nuclide == ba133:
+                start_index = 150
+                height = 0.10 * np.max(hist[start_index:])
+                prominence = 0.10 * np.max(hist[start_index:])
+                distance = 10
+            elif self.check_source.nuclide == mn54:
+                start_index = 400
+                height = 0.7 * np.max(hist[start_index:])
+                prominence = 0.7 * np.max(hist[start_index:])
+                distance = 100
+        else:
+            raise ValueError(
+                f"Unknown detector type: {self.detector_type}. Supported types are 'NaI' and 'HPGe'."
+            )
 
         # update the parameters if kwargs are provided
         if kwargs:
@@ -537,6 +585,10 @@ def get_peaks(self, hist: np.ndarray, **kwargs) -> np.ndarray:
         )
         peaks = np.array(peaks) + start_index
 
+        # special case for Mn-54, only keep the first high count energy peak
+        if self.check_source.nuclide == mn54 and len(peaks) > 1:
+            peaks = np.array([peaks[0]])
+
         return peaks
 
 
@@ -550,6 +602,7 @@ def get_gamma_emitted(
         calibration_coeffs,
         channel_nb: int,
         search_width: float = 800,
+        summing_method: str = 'sum_gaussian',
     ):
         # find right background detector
 
@@ -569,6 +622,7 @@ def get_gamma_emitted(
             calibrated_bin_edges,
             energy,
             search_width=search_width,
+            summing_method=summing_method,
         )
 
         nb_counts_measured = np.array(nb_counts_measured)
@@ -713,6 +767,7 @@ def get_calibration_data(
     check_source_measurements: List[CheckSourceMeasurement],
     background_measurement: Measurement,
     channel_nb: int,
+    peak_kwargs: dict = None,
 ):
     background_detector = [
         detector
@@ -731,15 +786,28 @@ def get_calibration_data(
             hist, bin_edges = detector.get_energy_hist_background_substract(
                 background_detector, bins=None
             )
-            peaks_ind = measurement.get_peaks(hist)
+            if peak_kwargs is not None:
+                if measurement.check_source.nuclide in peak_kwargs.keys():
+                    kwargs = peak_kwargs[measurement.check_source.nuclide]
+                else:
+                    kwargs = {}
+            else:
+                kwargs = {}
+
+            peaks_ind = measurement.get_peaks(hist, **kwargs)
             peaks = bin_edges[peaks_ind]
 
             if len(peaks) != len(measurement.check_source.nuclide.energy):
                 raise ValueError(
-                    f"SciPy find_peaks() found {len(peaks)} photon peaks, while {len(measurement.check_source.nuclide.energy)} were expected"
+                    f"SciPy find_peaks() found {len(peaks)} photon peaks, while {len(measurement.check_source.nuclide.energy)} were expected",
+                    f" peaks found: {peaks} for {measurement.check_source.nuclide.name}",
                 )
             calibration_channels += list(peaks)
             calibration_energies += measurement.check_source.nuclide.energy
+            # Store the uncalibrated measured energies in the measurement object for later use
+            if measurement.check_source.nuclide._uncalibrated_measured_energies is None:
+                measurement.check_source.nuclide._uncalibrated_measured_energies = {}
+            measurement.check_source.nuclide._uncalibrated_measured_energies[channel_nb] = list(peaks)
 
     inds = np.argsort(calibration_channels)
     calibration_channels = np.array(calibration_channels)[inds]
@@ -763,7 +831,9 @@ def gauss(x, b, m, *args):
     return out
 
 
-def fit_peak_gauss(hist, xvals, peak_ergs, search_width=600, threshold_overlap=200):
+def fit_peak_gauss(hist, xvals, peak_ergs, 
+                   search_width=600, 
+                   threshold_overlap=200):
 
     if len(peak_ergs) > 1:
         if np.max(peak_ergs) - np.min(peak_ergs) > threshold_overlap:
@@ -772,17 +842,18 @@ def fit_peak_gauss(hist, xvals, peak_ergs, search_width=600, threshold_overlap=2
             )
 
     search_start = np.argmin(
-        np.abs((peak_ergs[0] - search_width / (2 * len(peak_ergs))) - xvals)
+        np.abs((peak_ergs[0] - search_width / ( len(peak_ergs))) - xvals)
     )
     search_end = np.argmin(
-        np.abs((peak_ergs[-1] + search_width / (2 * len(peak_ergs))) - xvals)
+        np.abs((peak_ergs[-1] + search_width / (len(peak_ergs))) - xvals)
     )
 
     slope_guess = (hist[search_end] - hist[search_start]) / (
         xvals[search_end] - xvals[search_start]
     )
 
-    guess_parameters = [0, slope_guess]
+    # guess_parameters = [0, slope_guess]
+    guess_parameters = [0, 0]
 
     for i in range(len(peak_ergs)):
         peak_ind = np.argmin(np.abs((peak_ergs[i]) - xvals))
@@ -792,37 +863,51 @@ def fit_peak_gauss(hist, xvals, peak_ergs, search_width=600, threshold_overlap=2
             search_width / (3 * len(peak_ergs)),
         ]
 
+
     parameters, covariance = curve_fit(
         gauss,
         xvals[search_start:search_end],
         hist[search_start:search_end],
         p0=guess_parameters,
     )
 
+    print("Fitted parameters:", parameters)
+
     return parameters, covariance
 
 
 def get_multipeak_area(
-    hist, bins, peak_ergs, search_width=600, threshold_overlap=200
+    hist, 
+    bins, 
+    peak_ergs, 
+    search_width=600, 
+    threshold_overlap=200,
+    summing_method='sum_gaussian',
+    ax=None,
 ) -> List[float]:
+
+    print(peak_ergs)
 
     if len(peak_ergs) > 1:
         if np.max(peak_ergs) - np.min(peak_ergs) > threshold_overlap:
             areas = []
-            for peak in peak_ergs:
+            for p, peak in enumerate(peak_ergs):
                 area = get_multipeak_area(
                     hist,
                     bins,
                     [peak],
                     search_width=search_width,
                     threshold_overlap=threshold_overlap,
+                    summing_method=summing_method
                 )
                 areas += area
             return areas
+
 
-    # get midpoints of every bin
+    # get midpoints of every binß
     xvals = np.diff(bins) / 2 + bins[:-1]
 
+
     parameters, covariance = fit_peak_gauss(
         hist, xvals, peak_ergs, search_width=search_width
     )
@@ -845,17 +930,24 @@ def get_multipeak_area(
         # Use unimodal gaussian to estimate counts from just one peak
         peak_params = [parameters[0], parameters[1], parameters[2 + 3 * i], mean, sigma]
         all_peak_params += [peak_params]
-        gross_area = np.trapz(
-            gauss(xvals[peak_start:peak_end], *peak_params),
-            x=xvals[peak_start:peak_end],
-        )
 
+        if summing_method == 'sum_gaussian':
+            gross_area = np.trapz(
+                gauss(xvals[peak_start:peak_end], *peak_params),
+                x=xvals[peak_start:peak_end],
+            )
+        elif summing_method == 'sum_histogram':
+            gross_area = np.trapz(
+                hist[peak_start:peak_end],
+                x=xvals[peak_start:peak_end],
+            )
         # Cut off trapezoidal area due to compton scattering and noise
         trap_cutoff_area = np.trapz(
             parameters[0] + parameters[1] * xvals[peak_start:peak_end],
             x=xvals[peak_start:peak_end],
         )
         area = gross_area - trap_cutoff_area
+
         areas += [area]
 
     return areas
@@ -939,7 +1031,7 @@ def get_events(directory: str) -> Tuple[Dict[int, np.ndarray], Dict[int, np.ndar
         time_values[ch] = np.empty(0)
         energy_values[ch] = np.empty(0)
         for i, filename in enumerate(data_filenames[ch]):
-            print(f'Processing File {i}')
+            # print(f'Processing File {i}')
 
             csv_file_path = os.path.join(directory, filename)
 

test/neutron_detection/compass_test_data/background_measurement/calibration_coefficients.json

@@ -0,0 +1,18 @@
+{
+  "HPGe": {
+    "0": [
+      1.0198403099123785,
+      2.689961639523998
+    ]
+  },
+  "NaI": {
+    "4": [
+      0.7959901704589225,
+      -38.917607406086006
+    ],
+    "5": [
+      1.0622115094394577,
+      -30.576666718985596
+    ]
+  }
+}