From c64346092e27ce8df7af1ca9c175b65fadea2e31 Mon Sep 17 00:00:00 2001 From: fzeiser Date: Wed, 17 Jun 2020 10:47:23 +0200 Subject: [PATCH] Fix response interpolation above Compt.edge This closes #131. Changed response function interpolation between Compton edge and he chosen max. energy. Before, there was a misunderstanding of the *bin-by-bin interpolation* in Guttormsen1996. It now used a fan-like interpolation, too. Most noticable for response functions with small fwhm (like 1/10 Magne recommends for unfolding). --- ompy/response.py | 143 +++++++++++++++++++++++++++++++++++++---------- release_note.md | 51 +++++++++++++++++ setup.py | 4 +- 3 files changed, 165 insertions(+), 33 deletions(-) create mode 100644 release_note.md diff --git a/ompy/response.py b/ompy/response.py index 683517df..10c554bf 100644 --- a/ompy/response.py +++ b/ompy/response.py @@ -10,7 +10,7 @@ import numpy as np import pandas as pd from pathlib import Path -from typing import Union, Optional, Tuple +from typing import Union, Optional, Tuple, Dict, Any from scipy.interpolate import interp1d import logging @@ -144,7 +144,8 @@ def LoadZip(self, compton_matrix[i, 0:len(cmp_current)] = cmp_current i += 1 - return resp, compton_matrix, np.linspace(a0_cmp, a1_cmp * (N_cmp - 1), N_cmp) + return resp, compton_matrix, np.linspace(a0_cmp, a1_cmp * (N_cmp - 1), + N_cmp) def LoadDir(self, path: Union[str, Path], @@ -201,7 +202,6 @@ def LoadDir(self, # "Eg_sim" means "gamma, simulated", and refers to the gamma energies # where we have simulated Compton spectra. - # Get calibration and array length from highest-energy spectrum, # because the spectra may have differing length, # but the last is bound to be the longest. @@ -270,7 +270,7 @@ def interpolate(y, fill_value="extrapolate"): fwhm_rel_1330 = (self.fwhm_abs / 1330 * 100) self.f_fwhm_rel_perCent = interpolate(self.resp['FWHM_rel_norm'] * fwhm_rel_1330) - def f_fwhm_abs(E): + def f_fwhm_abs(E): # noqa return E * self.f_fwhm_rel_perCent(E)/100 self.f_fwhm_abs = f_fwhm_abs @@ -278,7 +278,8 @@ def f_fwhm_abs(E): def interpolate(self, Eout: np.ndarray = None, fwhm_abs: float = None, - return_table: bool = False) -> Union[Matrix, Tuple[Matrix, pd.DataFrame]]: + return_table: bool = False + ) -> Union[Matrix, Tuple[Matrix, pd.DataFrame]]: """ Interpolated the response matrix Perform the interpolation for the energy range specified in Eout with @@ -287,7 +288,13 @@ def interpolate(self, The interpolation is split into energy regions. Below the back-scattering energy Ebsc we interpolate linearly, then we apply the "fan method" (Guttormsen 1996) in the region from Ebsc up to the - Compton edge, then linear extrapolation again the rest of the way. + Compton edge, with a Compton scattering angle dependent interpolation. + From the Compton edge to Egmax we also use a fan, but with a linear + interpolation. + + Note: + Below the ~350 keV we only use a linear interpolation, as the + fan method does not work. This is not described in Guttormsen 1996. Args: folderpath: The path to the folder containing Compton spectra and @@ -319,25 +326,26 @@ def interpolate(self, # Loop over rows of the response matrix # TODO for speedup: Change this to a cython for j, E in enumerate(Eout): - - # Find maximal energy for current response (+n*sigma) function, - # -> Better if the lowest energies of the simulated spectra are - # above the gamma energy to be extrapolated - oneSigma = fwhm_abs * self.f_fwhm_rel_perCent_norm(E) / 2.35 + oneSigma = fwhm_abs_array[j] / 2.35 Egmax = E + 6 * oneSigma i_Egmax = min(index(Eout, Egmax), N_out-1) - # LOG.debug("Response for E: {E:.0f} calc. up to {Egmax:.0f}") + LOG.debug(f"Response for E {E:.0f} calc up to {Eout[i_Egmax]:.0f}") compton = self.get_closest_compton(E) R_low, i_bsc = self.linear_backscatter(E, compton) - R_fan, i_last = self.fan_method(E, compton, + R_fan, i_last = \ + self.fan_method_compton(E, compton, + i_start=i_bsc+1, i_stop=i_Egmax) + if i_last <= i_bsc+2: # fan method didn't do anything + R_high = self.linear_to_end(E, compton, i_start=i_bsc+1, i_stop=i_Egmax) - if i_last == i_bsc+1: # fan method didn't do anything - i_last -= 1 - R_high = self.linear_to_end(E, compton, - i_start=i_last+1, i_stop=i_Egmax) - R[j, :] += R_low + R_fan + R_high + R[j, :] += R_low + R_high + else: + R_high = self.fan_to_end(E, compton, + i_start=i_last+1, i_stop=i_Egmax, + fwhm_abs_array=fwhm_abs_array) + R[j, :] += R_low + R_fan + R_high # coorecton below E_sim[0] if E < self.resp['Eg'][0]: @@ -364,7 +372,8 @@ def interpolate(self, response = Matrix(values=R, Eg=Eout, Ex=Eout) if return_table: - # Return the response matrix, as well as the other structures, FWHM and efficiency, interpolated to the Eout_array + # Return the response matrix, as well as the other structures, + # FWHM and efficiency, interpolated to the Eout_array response_table = {'E': Eout, 'fwhm_abs': fwhm_abs_array, 'fwhm_rel_%': self.f_fwhm_rel_perCent(Eout), @@ -400,7 +409,7 @@ def iterpolate_checks(self): LOG.info(f"Note: Spectra outside of {self.resp['Eg'].min()} and " f"{self.resp['Eg'].max()} are extrapolation only.") - def get_closest_compton(self, E: float) -> Tuple[int, int]: + def get_closest_compton(self, E: float) -> Dict[str, Any]: """Find and rebin closest energies from available response functions If `E < self.resp['Eg'].min()` the compton matrix will be replaced @@ -410,7 +419,9 @@ def get_closest_compton(self, E: float) -> Tuple[int, int]: E (float): Description Returns: - Tuple[int, int]: `ilow` and `ihigh` are indices of closest energies + Dict with entries `Elow` and `Ehigh`, and `ilow` and `ihigh`, the + (indices) of closest energies. The arrays `counts_low` and + `counts_high` are the corresponding arrays of `compton` spectra. """ N = len(self.resp['Eg']) # ilow = 0 @@ -429,8 +440,10 @@ def get_closest_compton(self, E: float) -> Tuple[int, int]: Elow = self.resp['Eg'][ilow] cmp_low = self.cmp_matrix[ilow, :] - cmp_low = rebin_1D(cmp_low, self.Ecmp_array, self.Eout) - cmp_high = rebin_1D(cmp_high, self.Ecmp_array, self.Eout) + Enew = np.arange(self.Eout[0], self.Ecmp_array[-1], + self.Eout[1] - self.Eout[0]) + cmp_low = rebin_1D(cmp_low, self.Ecmp_array, Enew) + cmp_high = rebin_1D(cmp_high, self.Ecmp_array, Enew) compton = {"ilow": ilow, "ihigh": ihigh, @@ -442,7 +455,8 @@ def get_closest_compton(self, E: float) -> Tuple[int, int]: return compton def linear_cmp_interpolation(self, E: float, compton: dict, - fill_value: str = "extrapolate")-> np.ndarray: + fill_value: str = "extrapolate" + ) -> np.ndarray: """Linear interpolation between the compton spectra Args: @@ -505,8 +519,77 @@ def linear_to_end(self, E: float, compton: dict, R[R < 0] = 0 return R - def fan_method(self, E: float, compton: dict, - i_start: int, i_stop: int) -> Tuple[np.ndarray, int]: + def fan_to_end(self, E: float, compton: dict, + i_start: int, i_stop: int, + fwhm_abs_array: np.ndarray) -> np.ndarray: + """Linear(!) fan interpolation from Compton edge to Emax + + The fan-part is "scaled" by the distance between the Compton edge and + max(E). To get a reasonable scaling, we have to use ~6 sigma. + + Note: + We extrapolate up to self.N_out, and not i_stop, as a workaround + connected to Magne's 1/10th FWHM unfolding [which results + in a very small i_stop.] + + Args: + E (float): Incident energy + compton (dict): Dict. with information about the compton spectra + to interpolate between + i_start (int): Index where to start (usually end of fan method) + i_stop (int): Index where to stop (usually E+n*resolution) + fwhm_abs_array (np.ndarray): FHWM array, absolute values + + Returns: + np.ndarray: Response for `E` + """ + + R = np.zeros(self.N_out) + + Esim_low = compton["Elow"] + Esim_high = compton["Ehigh"] + Ecmp1 = self.E_compton(Esim_low, np.pi) + Ecmp2 = self.E_compton(Esim_high, np.pi) + i_low_edge = index(self.Eout, Ecmp1) + i_high_edge = index(self.Eout, Ecmp2) + + oneSigma = fwhm_abs_array[index(self.Eout, Esim_low)] / 2.35 + Egmax1 = Esim_low + 6 * oneSigma + scale1 = (Egmax1 - Ecmp1) / (self.Eout[i_stop] - self.Eout[i_start]) + + oneSigma = fwhm_abs_array[index(self.Eout, Esim_high)] / 2.35 + Egmax2 = Esim_high + 6 * oneSigma + scale2 = (Egmax2 - Ecmp2) / (self.Eout[i_stop] - self.Eout[i_start]) + + def lin_interpolation(x, x0, y0, x1, y1): + return y0 + (y1-y0)*(x-x0)/(x1-x0) + + i_edge = i_start-1 + # for i in range(i_edge+1, i_stop): + for i in range(i_edge+1, self.N_out): + i1 = int(i_low_edge + scale1 * (i - i_edge)) + i2 = int(i_high_edge + scale2 * (i - i_edge)) + + if i1 >= len(compton["counts_low"]): + i1 = len(compton["counts_low"])-1 + if i2 >= len(compton["counts_high"]): + i2 = len(compton["counts_high"])-1 + + c1 = compton["counts_low"][i1] + c2 = compton["counts_high"][i2] + y = lin_interpolation(E, Esim_low, c1, Esim_high, c2) + R[i] = y + + if len(R[R < 0]) != 0: + LOG.debug(f"In linear fan method, {len(R[R < 0])} entries in R" + "are negative and now set to 0") + R[R < 0] = 0 + + return R + + def fan_method_compton(self, E: float, compton: dict, + i_start: int, + i_stop: int) -> Tuple[np.ndarray, int]: """Fan method Args: @@ -628,21 +711,19 @@ def E_compton(Eg, theta): """ Calculates the energy of an electron that is scattered an angle theta by a gamma-ray of energy Eg. - Adapted from MAMA, file "folding.f", which references - Canberra catalog ed.7, p.2. Note: For `Eg <= 0.1` it returns `Eg`. (workaround) Args: - Eg: Energy of gamma-ray in keV + Eg: Energy of incident gamma-ray in keV theta: Angle of scatter in radians Returns: Energy Ee of scattered electron """ - scattered = Eg / (1 + Eg / 511 * (1 - np.cos(theta))) - electron = scattered * Eg / 511 * (1 - np.cos(theta)) + Eg_scattered = Eg / (1 + Eg / 511 * (1 - np.cos(theta))) + electron = Eg - Eg_scattered return np.where(Eg > 0.1, electron, Eg) @staticmethod diff --git a/release_note.md b/release_note.md new file mode 100644 index 00000000..4867f0cf --- /dev/null +++ b/release_note.md @@ -0,0 +1,51 @@ +# Release notes for OMpy + +## v.1.1.0 +Most important changes: +- Changed response function interpolation between Compton edge and he chosen max. energy. Before, there was a + misunderstanding of the *bin-by-bin interpolation* in Guttormsen1996. It now used a fan-like interpolation, + too. Most noticable for response functions with small fwhm (like 1/10 Magne recommends for unfolding). + +## v1.0.0 +Several changes, most important: + +**theoretical framework**: +- Corrected likelihood: Forgot non-constant K term when K=K(theta): c8c0046153b1eb269d00280b572f742b1a3cf4d7 + +**parameters and choices**: +- unfolding parameters: 0fcafe2ff7770be8c2bb107256201af79739cdb3 +- unfolder and fg method use remove negatives only, no fill: 9edb48537cca1f88c3120a73fa8eb92f6ebb5177 +- Randomize p0 for decomposition 77dec9db9a3a34d5fd6195752c84cfbca0c26c39 + +**implementation and convenience**: +- different save/load for vectors e5f7e52ce13cff04e8b23f50a00902be1d098bfc and parent commits +- Enable pickling of normalizer instances via dill: 896b352686594a8c7dbe52904645cc5b900ba800 + + +## v0.9.1 +Changes: + +- corrected version number +(v 0.9.0 has still v.0.8.0 as the version number) + +## v0.9 +Many changes to the API have occured since v.0.2. Here a (incomplete) summary of the main changes: + +- `Vector` and `Matrix` are now in mid-bin calibration. Most or all other functions have been adopted. +- Many changes (bugfix & readability) to the ensemble, decomposition and normalization classes. +- Normalization of nld and gsf ensembles working +- Parallelization, even though it could be more efficient for multinest (see #94 ) +- Renamed response functions submodule; run `git submodule update --init --recursive` after `git pull` to get the new files +- remember to run `pip install -e .` such that changes to the cython files will be recompiled +- Documentation now available via https://ompy.readthedocs.io +- Installation requirements are hopefully all specified; docker file is provided with integration at https://hub.docker.com/r/oslocyclotronlab/ompy and [mybinder](https://mybinder.org/v2/gh/oslocyclotronlab/ompy/master?filepath=ompy%2Fnotebooks%2Fgetting_started.ipynb) can be used to rund the examples. +- We have clean-up the history of the repo to downsize it. + Here the warning message: *NB: Read this (only) if you have cloned the repo before October 2019: We cleaned the repository from old comits clogging the repo (big data files that should never have been there). Unfortunetely, this has the sideeffect that the history had to be rewritten: Previous commits now have a different SHA1 (git version keys). If you need anything from the previous repo, see ompy_Archive_Sept2019. This will unfortunately also destroy references in issues. The simplest way to get the new repo is to rerun the installation instructions below.* + +## v0.2-beta +This is the first public beta version of the OMpy library, the Oslo Method in Python. + +**NB: Read this (only) if you have cloned the repo before October 2019 (which affects this release, v0.2-beta)**: +We cleaned the repository from old comits clogging the repo (big data files that should never have been there). Unfortunetely, this has the sideeffect that the history had to be rewritten: Previous commits now have a different SHA1 (git version keys). If you need anything from the previous repo, see [ompy_Archive_Sept2019](https://github.com/oslocyclotronlab/ompy_Archive_Sept2019). This will unfortunately also destroy references in issues. The simplest way to get the new repo is to rerun the installation instructions below. + +**In essence**: This tag does not work any longer; you have to download the version from https://zenodo.org/record/2654604 diff --git a/setup.py b/setup.py index a60720cb..99a108c7 100755 --- a/setup.py +++ b/setup.py @@ -20,8 +20,8 @@ # Version rutine taken from numpy MAJOR = 1 -MINOR = 0 -MICRO = 1 +MINOR = 1 +MICRO = 0 VERSION = '%d.%d.%d' % (MAJOR, MINOR, MICRO)