diff --git a/src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv b/src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv new file mode 100644 index 0000000..e71ec6d --- /dev/null +++ b/src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv @@ -0,0 +1,37 @@ +angle,662,1173,1332 +0,0.015,0.030,0.033 +10,0.011,0.021,0.024 +20,0.086,0.127,0.146 +30,0.294,0.356,0.397 +40,0.661,0.700,0.734 +50,1.054,1.057,1.057 +60,1.154,1.140,1.137 +70,1.186,1.152,1.138 +80,1.151,1.114,1.097 +90,1.000,1.000,1.000 +100,1.020,1.040,1.047 +110,1.074,1.093,1.103 +120,1.113,1.092,1.102 +130,1.139,1.122,1.113 +140,1.146,1.152,1.140 +150,1.113,1.118,1.104 +160,1.113,1.096,1.099 +170,1.091,1.076,1.083 +180,1.076,1.066,1.078 +-170,1.102,1.091,1.093 +-160,1.122,1.100,1.102 +-150,1.128,1.105,1.093 +-140,1.144,1.112,1.123 +-130,1.140,1.117,1.095 +-120,1.146,1.127,1.098 +-110,1.068,1.068,1.045 +-100,1.013,1.025,1.016 +-90,1.004,1.018,1.021 +-80,1.150,1.137,1.132 +-70,1.184,1.167,1.164 +-60,1.158,1.140,1.138 +-50,1.090,1.068,1.064 +-40,0.595,0.620,0.631 +-30,0.332,0.430,0.430 +-20,0.055,0.081,0.096 +-10,0.009,0.018,0.019 \ No newline at end of file diff --git a/src/pg_rad/data/detectors.csv b/src/pg_rad/data/detectors.csv new file mode 100644 index 0000000..f5ae726 --- /dev/null +++ b/src/pg_rad/data/detectors.csv @@ -0,0 +1,4 @@ +name,type,is_isotropic +dummy,NaI,true +LU_NaI_3inch,NaI,true +LU_HPGe_90,HPGe,false \ No newline at end of file diff --git a/src/pg_rad/data/field_efficiencies/LU_HPGe_90.csv b/src/pg_rad/data/field_efficiencies/LU_HPGe_90.csv new file mode 100644 index 0000000..c94f571 --- /dev/null +++ b/src/pg_rad/data/field_efficiencies/LU_HPGe_90.csv @@ -0,0 +1,17 @@ +energy_keV,field_efficiency_m2,uncertainty +59.5,0.00140,0.00005 +81.0,0.00310,0.00010 +122.1,0.00420,0.00013 +136.5,0.00428,0.00017 +160.6,0.00426,0.00030 +223.2,0.00418,0.00024 +276.4,0.00383,0.00012 +302.9,0.00370,0.00012 +356.0,0.00338,0.00010 +383.8,0.00323,0.00010 +511.0,0.00276,0.00008 +661.7,0.00241,0.00007 +834.8,0.00214,0.00007 +1173.2,0.00179,0.00005 +1274.5,0.00168,0.00005 +1332.5,0.00166,0.00005 \ No newline at end of file diff --git a/src/pg_rad/data/field_efficiencies/LU_NaI_3inch.csv b/src/pg_rad/data/field_efficiencies/LU_NaI_3inch.csv new file mode 100644 index 0000000..e3f4b1f --- /dev/null +++ b/src/pg_rad/data/field_efficiencies/LU_NaI_3inch.csv @@ -0,0 +1,64 @@ +energy_keV,field_efficiency_m2 +10,5.50129E-09 +11.22018454,2.88553E-07 +12.58925412,4.81878E-06 +14.12537545,3.55112E-05 +15.84893192,0.000146367 +17.7827941,0.000397029 +19.95262315,0.000803336 +22.38721139,0.00131657 +25.11886432,0.001862377 +28.18382931,0.002373449 +31.6227766,0.002811046 +33.164,0.00269554 +33.164,0.002698792 +35.48133892,0.002509993 +39.81071706,0.002801304 +44.66835922,0.003015877 +50.11872336,0.003227431 +56.23413252,0.00341077 +63.09573445,0.003562051 +70.79457844,0.00368852 +79.43282347,0.003788875 +89,0.003867423 +100,0.003925025 +112.2018454,0.003967222 +125.8925412,0.003991551 +141.2537545,0.004000729 +158.4893192,0.003993145 +177.827941,0.003969163 +199.5262315,0.003925289 +223.8721139,0.003856247 +251.1886432,0.00375596 +281.8382931,0.003619634 +316.227766,0.003446087 +354.8133892,0.003242691 +398.1071706,0.003021761 +446.6835922,0.002791816 +501.1872336,0.002568349 +562.3413252,0.002350052 +630.9573445,0.002147662 +707.9457844,0.001957893 +794.3282347,0.001785694 +891,0.001626634 +1000,0.001482571 +1122.018454,0.00135047 +1258.925412,0.001231358 +1412.537545,0.001116695 +1584.893192,0.001011833 +1778.27941,0.000917017 +1995.262315,0.000828435 +2238.721139,0.000746854 +2511.886432,0.000672573 +2818.382931,0.00060493 +3162.27766,0.000544458 +3548.133892,0.000488446 +3981.071706,0.000438438 +4466.835922,0.000392416 +5011.872336,0.00035092 +5623.413252,0.000313959 +6309.573445,0.000279409 +7079.457844,0.000247794 +7943.282347,0.000218768 +8913,0.000190209 +10000,0.000164309 \ No newline at end of file diff --git a/src/pg_rad/data/field_efficiencies/dummy.csv b/src/pg_rad/data/field_efficiencies/dummy.csv new file mode 100644 index 0000000..bd40665 --- /dev/null +++ b/src/pg_rad/data/field_efficiencies/dummy.csv @@ -0,0 +1,3 @@ +energy_keV,field_efficiency_m2 +0,1.0 +10000,1.0 \ No newline at end of file diff --git a/src/pg_rad/detector/builder.py b/src/pg_rad/detector/builder.py deleted file mode 100644 index 940a2f9..0000000 --- a/src/pg_rad/detector/builder.py +++ /dev/null @@ -1,20 +0,0 @@ -from pg_rad.inputparser.specs import DetectorSpec - -from .detectors import IsotropicDetector, AngularDetector - - -class DetectorBuilder: - def __init__( - self, - detector_spec: DetectorSpec, - ): - self.detector_spec = detector_spec - - def build(self) -> IsotropicDetector | AngularDetector: - if self.detector_spec.is_isotropic: - return IsotropicDetector( - self.detector_spec.name, - self.detector_spec.efficiency - ) - else: - raise NotImplementedError("Angular detector not supported yet.") diff --git a/src/pg_rad/detector/detector.py b/src/pg_rad/detector/detector.py new file mode 100644 index 0000000..89a6101 --- /dev/null +++ b/src/pg_rad/detector/detector.py @@ -0,0 +1,43 @@ +from importlib.resources import files + +from pandas import read_csv + +from pg_rad.utils.interpolators import ( + get_field_efficiency, get_angular_efficiency +) + + +class Detector: + def __init__( + self, + name: str, + type: str, + is_isotropic: bool + ): + self.name = name + self.type = type + self.is_isotropic = is_isotropic + + def get_efficiency(self, energy_keV, angle=None): + f_eff = get_field_efficiency(self.name, energy_keV) + + if self.is_isotropic or angle is None: + return f_eff + else: + f_eff = get_field_efficiency(self.name, energy_keV) + a_eff = get_angular_efficiency(self.name, energy_keV, *angle) + return f_eff * a_eff + + +def load_detector(name) -> Detector: + csv = files('pg_rad.data').joinpath('detectors.csv') + data = read_csv(csv) + dets = data['name'].values + if name in dets: + row = data[data['name'] == name].iloc[0] + return Detector(row['name'], row['type'], row['is_isotropic']) + else: + raise NotImplementedError( + f"Detector with name '{name}' not in detector library. Available:" + f"{', '.join(dets)}" + ) diff --git a/src/pg_rad/detector/detectors.py b/src/pg_rad/detector/detectors.py deleted file mode 100644 index 439ee56..0000000 --- a/src/pg_rad/detector/detectors.py +++ /dev/null @@ -1,38 +0,0 @@ -from abc import ABC - - -class BaseDetector(ABC): - def __init__( - self, - name: str, - efficiency: float - ): - self.name = name - self.efficiency = efficiency - - def get_efficiency(self): - pass - - -class IsotropicDetector(BaseDetector): - def __init__( - self, - name: str, - efficiency: float, - ): - super().__init__(name, efficiency) - - def get_efficiency(self, energy): - return self.efficiency - - -class AngularDetector(BaseDetector): - def __init__( - self, - name: str, - efficiency: float - ): - super().__init__(name, efficiency) - - def get_efficiency(self, angle, energy): - pass diff --git a/src/pg_rad/inputparser/parser.py b/src/pg_rad/inputparser/parser.py index c6dc8b2..9f4d981 100644 --- a/src/pg_rad/inputparser/parser.py +++ b/src/pg_rad/inputparser/parser.py @@ -353,41 +353,11 @@ class ConfigParser: return specs def _parse_detector(self) -> DetectorSpec: - det_dict = self.config.get("detector") - required = {"name", "is_isotropic"} - if not isinstance(det_dict, dict): - raise InvalidConfigValueError( - "detector is not specified correctly. Must contain at least" - f"the subkeys {required}" - ) + det_name = self.config.get("detector") + if not det_name: + raise MissingConfigKeyError("detector") - missing = required - det_dict.keys() - if missing: - raise MissingConfigKeyError("detector", missing) - - name = det_dict.get("name") - is_isotropic = det_dict.get("is_isotropic") - eff = det_dict.get("efficiency") - - default_detectors = defaults.DETECTOR_EFFICIENCIES - - if name in default_detectors.keys() and not eff: - eff = default_detectors[name] - elif eff: - pass - else: - raise InvalidConfigValueError( - f"The detector {name} not found in library. Either " - f"specify {name}.efficiency or " - "choose a detector from the following list: " - f"{default_detectors.keys()}." - ) - - return DetectorSpec( - name=name, - efficiency=eff, - is_isotropic=is_isotropic - ) + return DetectorSpec(name=det_name) def _warn_unknown_keys(self, section: str, provided: set, allowed: set): unknown = provided - allowed diff --git a/src/pg_rad/inputparser/specs.py b/src/pg_rad/inputparser/specs.py index db00bae..be3cd00 100644 --- a/src/pg_rad/inputparser/specs.py +++ b/src/pg_rad/inputparser/specs.py @@ -67,8 +67,6 @@ class RelativePointSourceSpec(PointSourceSpec): @dataclass class DetectorSpec: name: str - efficiency: float - is_isotropic: bool @dataclass diff --git a/src/pg_rad/utils/interpolators.py b/src/pg_rad/utils/interpolators.py new file mode 100644 index 0000000..d07bd9c --- /dev/null +++ b/src/pg_rad/utils/interpolators.py @@ -0,0 +1,56 @@ +from importlib.resources import files + +import numpy as np +from pandas import read_csv +from scipy.interpolate import interp1d, CubicSpline + +from pg_rad.configs.filepaths import ATTENUATION_TABLE + + +def get_mass_attenuation_coeff(*args) -> float: + csv = files('pg_rad.data').joinpath(ATTENUATION_TABLE) + data = read_csv(csv) + x = data["energy_mev"].to_numpy() + y = data["mu"].to_numpy() + f = interp1d(x, y) + return f(*args) + + +def get_field_efficiency(name: str, energy_keV: float) -> float: + csv = files('pg_rad.data.field_efficiencies').joinpath(name+'.csv') + data = read_csv(csv) + data = data.groupby("energy_keV", as_index=False).mean() + x = data["energy_keV"].to_numpy() + y = data["field_efficiency_m2"].to_numpy() + f = CubicSpline(x, y) + return f(energy_keV) + + +def get_angular_efficiency(name: str, energy_keV: float, *angle: float): + csv = files('pg_rad.data.angular_efficiencies').joinpath(name+'.csv') + data = read_csv(csv) + + # check all energies at which angular eff. is available for this detector. + # this is done within 1% tolerance + energy_cols = [col for col in data.columns if col != "angle"] + energies = np.array([float(col) for col in energy_cols]) + rel_diff = np.abs(energies - energy_keV) / energies + match_idx = np.where(rel_diff <= 0.01)[0] + + if len(match_idx) == 0: + raise NotImplementedError( + f"No angular efficiency defined for {energy_keV} keV " + f"in detector '{name}'. Available: {energies}" + ) + + best_idx = match_idx[np.argmin(rel_diff[match_idx])] + selected_energy_col = energy_cols[best_idx] + + x = data["angle"].to_numpy() + y = data[selected_energy_col].to_numpy() + idx = np.argsort(x) + x = x[idx] + y = y[idx] + f = interp1d(x, y) + + return f(angle)