From 1671c4ec07337beae5b9a89fe5a2865a748c92c9 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Fri, 30 Jan 2026 11:00:17 -0600 Subject: [PATCH 1/8] * Warn if spectral normalization results in zero incident photons for any Ei bin when computing effective area correction for response. * Don't allow inf or nan in effective area in the above case -- set the correction to zero for these bins. * Actually honor emin and emax in Linear normalization, as it is honored for the powerlaw case * Consistently strip units when computing effective area correction * Try to correct (currently untestable) validation code for gaussian spectral normalization --- cosipy/response/RspConverter.py | 35 +++++++++++++++------- tests/response/test_response_conversion.py | 33 +++++++++++++++----- 2 files changed, 51 insertions(+), 17 deletions(-) diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index 49f257cb..9d5b56ae 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -436,23 +436,22 @@ def _get_eff_area_correction(self, axes, hdr): """ - ewidth = axes['Ei'].widths - norm = hdr["norm"] # If we have one single bin, treat the Gaussian normalization # like the mono one. Also check that the Gaussian spectrum is # fully contained in that bin - if norm == "Gaussian" and len(ewidth) == 1: + if norm == "Gaussian" and axes['Ei'].nbins == 1: from scipy.special import erf Gauss_mean = hdr["norm_params"][0] + Gauss_sdev = hdr["norm_params"][1] + + edges = axes['Ei'].edges.value # only two edges for 1 bin - edges = axes['Ei'].edges - gauss_int = \ - 0.5 * (1 + erf( (edges[0] - Gauss_mean)/(4*np.sqrt(2)) ) ) + \ - 0.5 * (1 + erf( (edges[1] - Gauss_mean)/(4*np.sqrt(2)) ) ) + z = (edges - Gauss_mean)/(Gauss_sdev * np.sqrt(2)) + gauss_int = np.diff(0.5*(1 + erf(z))) assert gauss_int == 1, "The gaussian spectrum is not fully contained in this single bin!" @@ -470,7 +469,16 @@ def _get_eff_area_correction(self, axes, hdr): if not self.quiet: logger.info(f"normalization: linear with energy range [{emin}-{emax}]") - nperchannel_norm = ewidth / (emax - emin) + e_lo = axes['Ei'].lower_bounds.value + e_hi = axes['Ei'].upper_bounds.value + + e_lo = np.minimum(emax, e_lo) + e_hi = np.minimum(emax, e_hi) + + e_lo = np.maximum(emin, e_lo) + e_hi = np.maximum(emin, e_hi) + + nperchannel_norm = (e_hi - e_lo) / (emax - emin) case "Mono" : if not self.quiet: @@ -484,7 +492,6 @@ def _get_eff_area_correction(self, axes, hdr): if not self.quiet: logger.info(f"normalization: powerlaw with index {alpha} with energy range [{emin}-{emax}]keV") - # From powerlaw e_lo = axes['Ei'].lower_bounds.value e_hi = axes['Ei'].upper_bounds.value @@ -507,8 +514,16 @@ def _get_eff_area_correction(self, axes, hdr): # We assume all FISBEL pixels have the same area. nperchannel = nperchannel_norm * hdr["nevents_sim"] / axes["NuLambda"].nbins + zero_weights = (nperchannel == 0.) + + if np.any(zero_weights): + logger.warning("Spectral normalization gives zero incident photons in some Ei bins; " + "eff_area correction for those bins will be set to zero!") + # Area - eff_area = hdr["area_sim"] / nperchannel + eff_area = np.zeros(len(nperchannel)) + eff_area[~zero_weights] = hdr["area_sim"] / nperchannel[~zero_weights] + return eff_area diff --git a/tests/response/test_response_conversion.py b/tests/response/test_response_conversion.py index 8d0f483c..200c9e22 100644 --- a/tests/response/test_response_conversion.py +++ b/tests/response/test_response_conversion.py @@ -1,5 +1,7 @@ +import numpy as np + from cosipy import test_data -7 + from cosipy.response import FullDetectorResponse, RspConverter rspgz_response_path = test_data.path / "test_full_detector_response.rsp.gz" @@ -54,7 +56,7 @@ def test_default_norms(tmp_path): c = RspConverter(bufsize = 100000, norm = "Linear", - norm_params = [50, 1000]) + norm_params = [50, 10000]) c.convert_to_h5(rspgz_nonorm_response_path, h5_filename = tmp_h5_filename, @@ -62,7 +64,7 @@ def test_default_norms(tmp_path): fdr = FullDetectorResponse.open(tmp_h5_filename) norm_info = fdr.headers["SP"] - assert norm_info == "Linear 50 1000" + assert norm_info == "Linear 50 10000" fdr.close() c = RspConverter(bufsize = 100000, @@ -79,7 +81,20 @@ def test_default_norms(tmp_path): c = RspConverter(bufsize = 100000, norm = "powerlaw", - norm_params = [50, 1000, 0.9]) + norm_params = [50, 10000, 0.9]) + + c.convert_to_h5(rspgz_nonorm_response_path, + h5_filename = tmp_h5_filename, + overwrite = True) + + fdr = FullDetectorResponse.open(tmp_h5_filename) + norm_info = fdr.headers["SP"] + assert norm_info == "powerlaw 50 10000 0.9" + fdr.close() + + c = RspConverter(bufsize = 100000, + norm = "powerlaw", + norm_params = [50, 10000, 1]) c.convert_to_h5(rspgz_nonorm_response_path, h5_filename = tmp_h5_filename, @@ -87,12 +102,12 @@ def test_default_norms(tmp_path): fdr = FullDetectorResponse.open(tmp_h5_filename) norm_info = fdr.headers["SP"] - assert norm_info == "powerlaw 50 1000 0.9" + assert norm_info == "powerlaw 50 10000 1.0" fdr.close() c = RspConverter(bufsize = 100000, norm = "powerlaw", - norm_params = [50, 1000, 1]) + norm_params = [500, 1000, 1]) c.convert_to_h5(rspgz_nonorm_response_path, h5_filename = tmp_h5_filename, @@ -100,7 +115,11 @@ def test_default_norms(tmp_path): fdr = FullDetectorResponse.open(tmp_h5_filename) norm_info = fdr.headers["SP"] - assert norm_info == "powerlaw 50 1000 1.0" + assert norm_info == "powerlaw 500 1000 1.0" + + # bins with no defined spectral normalization should + # have their eff_area set to zero, not -inf or NaN + assert all(np.isfinite(fdr.eff_area_correction)) fdr.close() # cannot test Gaussian norm because it can only be used on From dfe44eeb527e59e83a8ee941f1bd1e45900bd012 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Mon, 2 Feb 2026 11:50:05 -0600 Subject: [PATCH 2/8] * don't try to guess element type of response file until we know that it's a valid file * slightly more robust stripping of extensions from input filename (but we still don't support uncompressed .rsp, though perhaps we should) --- cosipy/response/RspConverter.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index 9d5b56ae..ef7b1a83 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -121,7 +121,7 @@ def convert_to_h5(self, Parameters ---------- rsp_filename: string - name of input file (must end with .rsp.gz) + name of input file (must end with .rsp or .rsp.gz) h5_filename : string (optional) name of output file (should end with .h5); if not specified, use base name of rsp_filename with .h5 extension @@ -142,20 +142,27 @@ def convert_to_h5(self, """ if h5_filename is None: - h5_filename = str(rsp_filename).replace(".rsp.gz", ".h5") + rsp_path = Path(rsp_filename) + + # strip .rsp and .gz from end of path and add .h5 + h5_path = rsp_path.parent / rsp_path.stem + while h5_path.suffix in {".rsp", "gz"}: + h5_path = h5_path.with_suffix("") + + h5_filename = str(rsp_filename) + ".h5" if Path(h5_filename).exists() and not overwrite: raise RuntimeError(f"Not overwriting existing HDF5 file {h5_filename}") - if elt_type is None: - elt_type = self._get_min_elt_type(rsp_filename) - # read all info from the .rsp file with gzip.open(rsp_filename, "rt") as f: axes, hdr = self._read_response_header(f) eff_area = self._get_eff_area_correction(axes, hdr) + if elt_type is None: + elt_type = self._get_min_elt_type(rsp_filename) + nbins = hdr["nbins"] counts = self._read_counts(f, axes, nbins, elt_type) @@ -678,6 +685,7 @@ def _get_min_elt_type(self, rsp_filename): with gzip.open(rsp_filename, "rt") as rsp_file: + nbins = None for line in rsp_file: # consume the file header line = line.split() @@ -685,6 +693,9 @@ def _get_min_elt_type(self, rsp_filename): nbins = int(line[1]) break + if nbins is None: + ValueError("Could not find start of data in .rsp file") + tq = tqdm(total=nbins, desc="Getting type for counts", disable=self.quiet) From f00a02e76b9a1e413afa9db76d3fb384dad29872 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Mon, 2 Feb 2026 12:28:46 -0600 Subject: [PATCH 3/8] * Transparently handle both compressed and uncompressed .rsp files, for both reading and writing. * Add test cases for uncompressed read/write * FDR error message should not mention .rsp, as it only opens .h5 files --- cosipy/response/FullDetectorResponse.py | 2 +- cosipy/response/RspConverter.py | 62 ++++++++++++++++------ tests/response/test_response_conversion.py | 51 ++++++++++++++++-- 3 files changed, 95 insertions(+), 20 deletions(-) diff --git a/cosipy/response/FullDetectorResponse.py b/cosipy/response/FullDetectorResponse.py index 1463acc9..012d1d35 100644 --- a/cosipy/response/FullDetectorResponse.py +++ b/cosipy/response/FullDetectorResponse.py @@ -67,7 +67,7 @@ def open(cls, filename, dtype=None, pa_convention=None, cache_size=None): return cls._open_h5(filename, dtype, pa_convention, cache_size) else: raise ValueError( - "Unsupported file format. Only .h5 and .rsp.gz extensions are supported.") + "Unsupported file format. Only .h5 extension is supported.") @classmethod def _open_h5(cls, filename, dtype=None, pa_convention=None, cache_size=None): diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index ef7b1a83..58e8b20d 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -19,13 +19,16 @@ class RspConverter(): """ - Converter between response files stored in .rsp.gz format and + Converter between response files stored in .rsp format and optimized HDF5 format on disk. - Use method convert_to_h5() to convert a .rsp.gz file to .h5. + Use method convert_to_h5() to convert a .rsp file to .h5. Use method convert_to_rsp() to convert a FullDetectorResponse - (backed by an .h5 file) to .rsp.gz. + (backed by an .h5 file) to .rsp. + + To read and write compressed .rsp files, add ".gz" to the end of + the .rsp filename. """ @@ -107,6 +110,33 @@ def __init__(self, self.norm = norm self.norm_params = norm_params + @staticmethod + def _open_rsp(rsp_filename, mode): + """ + Open an .rsp file with or without .gz extension. If .gz, + use gzip to open; otherwise, open as a regular file. + The resulting file object can be used in a context manager. + + Parameters + ---------- + rsp_filename : Path or string + name of file + mode : str + mode in which to open file + + Returns + ------- + open file or gzip stream object + + """ + + rsp_path = Path(rsp_filename) + + if rsp_path.suffix == ".gz": + return gzip.open(rsp_filename, mode) + else: + return open(rsp_filename, mode) + def convert_to_h5(self, rsp_filename, h5_filename = None, @@ -115,7 +145,7 @@ def convert_to_h5(self, elt_type = None): """ - Given a response file in .rsp.gz format, read it + Given a response file in .rsp format, read it and write it out as an HDF5 file Parameters @@ -144,7 +174,7 @@ def convert_to_h5(self, if h5_filename is None: rsp_path = Path(rsp_filename) - # strip .rsp and .gz from end of path and add .h5 + # strip any .rsp and .gz from end of path and add .h5 h5_path = rsp_path.parent / rsp_path.stem while h5_path.suffix in {".rsp", "gz"}: h5_path = h5_path.with_suffix("") @@ -155,7 +185,7 @@ def convert_to_h5(self, raise RuntimeError(f"Not overwriting existing HDF5 file {h5_filename}") # read all info from the .rsp file - with gzip.open(rsp_filename, "rt") as f: + with self._open_rsp(rsp_filename, "rt") as f: axes, hdr = self._read_response_header(f) eff_area = self._get_eff_area_correction(axes, hdr) @@ -200,7 +230,7 @@ def _read_response_header(self, rsp_file): Parameters ---------- - rsp_file : file handle to open .rsp.gz file + rsp_file : file handle to open .rsp file Returns ------- @@ -683,7 +713,7 @@ def _get_min_elt_type(self, rsp_filename): """ - with gzip.open(rsp_filename, "rt") as rsp_file: + with self._open_rsp(rsp_filename, "rt") as rsp_file: nbins = None for line in rsp_file: @@ -740,7 +770,7 @@ def _read_counts(self, rsp_file, axes, nbins, elt_type): Parameters ---------- rsp_file : file handle - open .rsp.gz file + open .rsp file nbins : int number of bins to read from file axes : Axes object @@ -839,7 +869,7 @@ def convert_to_rsp(self, overwrite = False): """ Convert a FullDetectorResponse object backed by an HDF5 file - into a textual .rsp.gz response. We reuse the header + into a textual .rsp or .rsp.gz response. We reuse the header information stored in the HDF5 file, along with its axes and counts. @@ -848,14 +878,15 @@ def convert_to_rsp(self, fullDetectorResponse : FullDetectorResponse object to be converted rsp_filename : string - path to write .rsp.gz file (should end with .rsp.gz) + path to write .rsp file; if extension is .rsp.gz, the file + will be gzipped. overwrite : bool if true, overwrite existing response if it exists """ if Path(rsp_filename).exists() and not overwrite: - raise RuntimeError(f"Not overwriting existing .rsp.gz file {rsp_filename}") + raise RuntimeError(f"Not overwriting existing file {rsp_filename}") # reorder axes if needed to match the expected order for an .rsp file axes = fullDetectorResponse._axes @@ -875,7 +906,7 @@ def convert_to_rsp(self, def _write_rsp(self, headers, axes, counts, rsp_filename): """ - Write an .rsp.gz file with all necessary info. + Write an .rsp file with all necessary info. Parameters ---------- @@ -886,7 +917,8 @@ def _write_rsp(self, headers, axes, counts, rsp_filename): counts : counts of Histogram rsp_filename : - name of response file to write (should be .rsp.gz). + path to write .rsp file; if extension is .rsp.gz, the file + will be gzipped. """ @@ -895,7 +927,7 @@ def _write_rsp(self, headers, axes, counts, rsp_filename): for desc in RspConverter.axis_name_map: axis_names[RspConverter.axis_name_map[desc]] = desc - with gzip.open(rsp_filename, "wt") as f: + with self._open_rsp(rsp_filename, "wt") as f: f.write("# computed reduced response\n") diff --git a/tests/response/test_response_conversion.py b/tests/response/test_response_conversion.py index 200c9e22..c1108aca 100644 --- a/tests/response/test_response_conversion.py +++ b/tests/response/test_response_conversion.py @@ -13,15 +13,17 @@ def test_convert_rsp_to_h5(tmp_path): + import gzip + tmp_h5_filename = tmp_path / "fdr.h5" c = RspConverter(bufsize = 100000) + # test opening compressed .rsp and writing compressed .h5 c.convert_to_h5(rspgz_response_path, h5_filename = tmp_h5_filename, overwrite = True) - fdr = FullDetectorResponse.open(h5_response_path) fdr2 = FullDetectorResponse.open(tmp_h5_filename) @@ -33,7 +35,7 @@ def test_convert_rsp_to_h5(tmp_path): assert h1 == h2 - + # test opening compressed .rsp and writing uncompressed .h5 c.convert_to_h5(rspgz_response_path, h5_filename = tmp_h5_filename, overwrite = True, @@ -50,6 +52,27 @@ def test_convert_rsp_to_h5(tmp_path): assert h1 == h2 + # test opening uncompressed .rsp + with gzip.open(rspgz_response_path) as gzfile: + content = gzfile.read() + with open(tmp_path / "response.rsp", "wb") as rspfile: + rspfile.write(content) + + c.convert_to_h5(rspgz_response_path, + h5_filename = tmp_h5_filename, + overwrite = True) + + fdr = FullDetectorResponse.open(h5_response_path) + fdr2 = FullDetectorResponse.open(tmp_h5_filename) + + h1 = fdr.to_dr() + h2 = fdr2.to_dr() + + fdr.close() + fdr2.close() + + assert h1 == h2 + def test_default_norms(tmp_path): tmp_h5_filename = tmp_path / "fdr.h5" @@ -127,12 +150,32 @@ def test_default_norms(tmp_path): def test_convert_h5_to_rsp(tmp_path): - tmp_rsp_filename = tmp_path / "fdr.rsp.gz" + tmp_rsp_filename = tmp_path / "fdr.rsp" + tmp_rspgz_filename = tmp_path / "fdr.rsp.gz" tmp_h5_filename = tmp_path / "fdr.h5" + c = RspConverter(bufsize = 100000) + + # test writing compressed .rsp.gz fdr = FullDetectorResponse.open(h5_response_path) - c = RspConverter(bufsize = 100000) + c.convert_to_rsp(fdr, tmp_rspgz_filename, overwrite=True) + + tmp_h5_filename = c.convert_to_h5(tmp_rspgz_filename, overwrite=True) + + fdr2 = FullDetectorResponse.open(tmp_h5_filename) + + + h1 = fdr.to_dr() + h2 = fdr2.to_dr() + + fdr.close() + fdr2.close() + + assert h1 == h2 + + # test writing uncompressed .rsp + fdr = FullDetectorResponse.open(h5_response_path) c.convert_to_rsp(fdr, tmp_rsp_filename, overwrite=True) From bbfba4dfdba44e4bf13417975e7968a29a27af77 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Tue, 3 Feb 2026 08:45:48 -0600 Subject: [PATCH 4/8] * Mono normalization takes an *optional* parameter, which is the mono-energetic energy. If not specified, normalization is ill-defined if Ei axis has multiple bins. If specified, normalization is zero for bins not containing the energy (which is weird and throws a warning, but at least is well-defined and should not crash). * test case for Mono in fact used an Ei axis with multiple bins but did not specify an energy! Fix to provide an energy. --- cosipy/response/RspConverter.py | 43 +++++++++++++--------- tests/response/test_response_conversion.py | 5 ++- 2 files changed, 29 insertions(+), 19 deletions(-) diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index 58e8b20d..b75eefd4 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -400,9 +400,10 @@ def _validate_norm_params(self, norm, params): match norm: case 'Mono': - if len(params) > 0: - raise ValueError(f"Mono normalization takes zero params; {len(params)} given") - params = () + if len(params) > 1: + raise ValueError(f"Mono normalization takes a most one param; {len(params)} given") + # emono if given, else dummy value + params = (None,) if len(params) == 0 else ( int(params[0]), ) case 'Linear': if len(params) != 2: @@ -474,18 +475,22 @@ def _get_eff_area_correction(self, axes, hdr): """ norm = hdr["norm"] + params = hdr["norm_params"] + + ei_axis = axes['Ei'] + e_lo = ei_axis.lower_bounds.value + e_hi = ei_axis.upper_bounds.value # If we have one single bin, treat the Gaussian normalization # like the mono one. Also check that the Gaussian spectrum is # fully contained in that bin - if norm == "Gaussian" and axes['Ei'].nbins == 1: + if norm == "Gaussian" and ei_axis.nbins == 1: from scipy.special import erf - Gauss_mean = hdr["norm_params"][0] - Gauss_sdev = hdr["norm_params"][1] + Gauss_mean, Gauss_sdev = params - edges = axes['Ei'].edges.value # only two edges for 1 bin + edges = ei_axis.edges.value # only two edges for 1 bin z = (edges - Gauss_mean)/(Gauss_sdev * np.sqrt(2)) gauss_int = np.diff(0.5*(1 + erf(z))) @@ -496,19 +501,16 @@ def _get_eff_area_correction(self, axes, hdr): logger.info("Only one bin so we will use the Mono normalisation") norm = "Mono" + params = (None,) match norm: case "Linear": - - emin, emax = hdr["norm_params"] + emin, emax = params if not self.quiet: logger.info(f"normalization: linear with energy range [{emin}-{emax}]") - e_lo = axes['Ei'].lower_bounds.value - e_hi = axes['Ei'].upper_bounds.value - e_lo = np.minimum(emax, e_lo) e_hi = np.minimum(emax, e_hi) @@ -521,17 +523,24 @@ def _get_eff_area_correction(self, axes, hdr): if not self.quiet: logger.info("normalization: mono") - nperchannel_norm = np.array([1.]) + emono = params[0] + if emono is None: + if ei_axis.nbins > 1: + raise ValueError("Cannot specify Mono norm without energy for Ei axis with multiple bins") + else: + # all energy is in the single bin + nperchannel_norm = np.array([1.]) + else: + # set just the Ei bin containing the mono energy to 1 + nperchannel_norm = np.zeros(ei_axis.nbins) + nperchannel_norm[(e_lo <= emono) & (e_hi >= emono)] = 1. case "powerlaw": - emin, emax, alpha = hdr["norm_params"] + emin, emax, alpha = params if not self.quiet: logger.info(f"normalization: powerlaw with index {alpha} with energy range [{emin}-{emax}]keV") - e_lo = axes['Ei'].lower_bounds.value - e_hi = axes['Ei'].upper_bounds.value - e_lo = np.minimum(emax, e_lo) e_hi = np.minimum(emax, e_hi) diff --git a/tests/response/test_response_conversion.py b/tests/response/test_response_conversion.py index c1108aca..45221e3a 100644 --- a/tests/response/test_response_conversion.py +++ b/tests/response/test_response_conversion.py @@ -91,7 +91,8 @@ def test_default_norms(tmp_path): fdr.close() c = RspConverter(bufsize = 100000, - norm = "Mono") + norm = "Mono", + norm_params = [511]) c.convert_to_h5(rspgz_nonorm_response_path, h5_filename = tmp_h5_filename, @@ -99,7 +100,7 @@ def test_default_norms(tmp_path): fdr = FullDetectorResponse.open(tmp_h5_filename) norm_info = fdr.headers["SP"] - assert norm_info == "Mono" + assert norm_info == "Mono 511" fdr.close() c = RspConverter(bufsize = 100000, From 47a8720b4f722d07d77dd0a9e0f1c64e1fc55420 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Tue, 3 Feb 2026 09:04:29 -0600 Subject: [PATCH 5/8] test more cases of Mono norm; we can't test one-Ei-bin case without a suitable response file --- tests/response/test_response_conversion.py | 61 ++++++++++++++-------- 1 file changed, 38 insertions(+), 23 deletions(-) diff --git a/tests/response/test_response_conversion.py b/tests/response/test_response_conversion.py index 45221e3a..00fb8e0a 100644 --- a/tests/response/test_response_conversion.py +++ b/tests/response/test_response_conversion.py @@ -4,6 +4,8 @@ from cosipy.response import FullDetectorResponse, RspConverter +from pytest import raises + rspgz_response_path = test_data.path / "test_full_detector_response.rsp.gz" rspgz_nonorm_response_path = test_data.path / "test_full_detector_response_no_norm.rsp.gz" @@ -85,10 +87,9 @@ def test_default_norms(tmp_path): h5_filename = tmp_h5_filename, overwrite = True) - fdr = FullDetectorResponse.open(tmp_h5_filename) - norm_info = fdr.headers["SP"] - assert norm_info == "Linear 50 10000" - fdr.close() + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "Linear 50 10000" c = RspConverter(bufsize = 100000, norm = "Mono", @@ -98,10 +99,26 @@ def test_default_norms(tmp_path): h5_filename = tmp_h5_filename, overwrite = True) - fdr = FullDetectorResponse.open(tmp_h5_filename) - norm_info = fdr.headers["SP"] - assert norm_info == "Mono 511" - fdr.close() + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "Mono 511" + + # bins with no defined spectral normalization should + # have their eff_area set to zero, not -inf or NaN + assert all(np.isfinite(fdr.eff_area_correction)) + + # if no monoenergetic energy is given, Mono is + # undefined with multiple Ei bins + with raises(ValueError): + c = RspConverter(bufsize = 100000, + norm = "Mono") + + c.convert_to_h5(rspgz_nonorm_response_path, + h5_filename = tmp_h5_filename, + overwrite = True) + + # cannot test Mono with no param without an .rsp file + # containing only a single Ei axis bin c = RspConverter(bufsize = 100000, norm = "powerlaw", @@ -111,10 +128,9 @@ def test_default_norms(tmp_path): h5_filename = tmp_h5_filename, overwrite = True) - fdr = FullDetectorResponse.open(tmp_h5_filename) - norm_info = fdr.headers["SP"] - assert norm_info == "powerlaw 50 10000 0.9" - fdr.close() + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "powerlaw 50 10000 0.9" c = RspConverter(bufsize = 100000, norm = "powerlaw", @@ -124,11 +140,11 @@ def test_default_norms(tmp_path): h5_filename = tmp_h5_filename, overwrite = True) - fdr = FullDetectorResponse.open(tmp_h5_filename) - norm_info = fdr.headers["SP"] - assert norm_info == "powerlaw 50 10000 1.0" - fdr.close() + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "powerlaw 50 10000 1.0" + # powerlaw norm that does not span all of Ei c = RspConverter(bufsize = 100000, norm = "powerlaw", norm_params = [500, 1000, 1]) @@ -137,14 +153,13 @@ def test_default_norms(tmp_path): h5_filename = tmp_h5_filename, overwrite = True) - fdr = FullDetectorResponse.open(tmp_h5_filename) - norm_info = fdr.headers["SP"] - assert norm_info == "powerlaw 500 1000 1.0" + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "powerlaw 500 1000 1.0" - # bins with no defined spectral normalization should - # have their eff_area set to zero, not -inf or NaN - assert all(np.isfinite(fdr.eff_area_correction)) - fdr.close() + # bins with no defined spectral normalization should + # have their eff_area set to zero, not -inf or NaN + assert all(np.isfinite(fdr.eff_area_correction)) # cannot test Gaussian norm because it can only be used on # a response with a single Ei bin From 64ebc3ccaa15d47a97a4f0ff82f54a6e02029335 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Wed, 4 Feb 2026 08:56:21 -0600 Subject: [PATCH 6/8] * fix output renaming for .h5 after compressed/uncompressed .rsp support * Gaussian normalization takes 3 arguments --- cosipy/response/RspConverter.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index b75eefd4..810b4033 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -179,7 +179,7 @@ def convert_to_h5(self, while h5_path.suffix in {".rsp", "gz"}: h5_path = h5_path.with_suffix("") - h5_filename = str(rsp_filename) + ".h5" + h5_filename = h5_path.parent / (h5_path.stem + ".h5") if Path(h5_filename).exists() and not overwrite: raise RuntimeError(f"Not overwriting existing HDF5 file {h5_filename}") @@ -488,11 +488,11 @@ def _get_eff_area_correction(self, axes, hdr): from scipy.special import erf - Gauss_mean, Gauss_sdev = params + mean, sdev, cutoff = params edges = ei_axis.edges.value # only two edges for 1 bin - z = (edges - Gauss_mean)/(Gauss_sdev * np.sqrt(2)) + z = (edges - mean)/(sdev * np.sqrt(2)) gauss_int = np.diff(0.5*(1 + erf(z))) assert gauss_int == 1, "The gaussian spectrum is not fully contained in this single bin!" From d341fc607569b835bfe4497c95b3a57c93115c4c Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Wed, 4 Feb 2026 13:53:04 -0600 Subject: [PATCH 7/8] support Gaussian normalization --- cosipy/response/RspConverter.py | 39 +++++++++------------- tests/response/test_response_conversion.py | 19 +++++++++-- 2 files changed, 33 insertions(+), 25 deletions(-) diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index 810b4033..7775e4e7 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -474,6 +474,11 @@ def _get_eff_area_correction(self, axes, hdr): """ + def gauss_int(x, mu, sigma): + from scipy.special import erf + z = (x - mu)/(sigma * np.sqrt(2)) + return 0.5*(1 + erf(z)) + norm = hdr["norm"] params = hdr["norm_params"] @@ -481,28 +486,6 @@ def _get_eff_area_correction(self, axes, hdr): e_lo = ei_axis.lower_bounds.value e_hi = ei_axis.upper_bounds.value - # If we have one single bin, treat the Gaussian normalization - # like the mono one. Also check that the Gaussian spectrum is - # fully contained in that bin - if norm == "Gaussian" and ei_axis.nbins == 1: - - from scipy.special import erf - - mean, sdev, cutoff = params - - edges = ei_axis.edges.value # only two edges for 1 bin - - z = (edges - mean)/(sdev * np.sqrt(2)) - gauss_int = np.diff(0.5*(1 + erf(z))) - - assert gauss_int == 1, "The gaussian spectrum is not fully contained in this single bin!" - - if not self.quiet: - logger.info("Only one bin so we will use the Mono normalisation") - - norm = "Mono" - params = (None,) - match norm: case "Linear": @@ -554,7 +537,17 @@ def _get_eff_area_correction(self, axes, hdr): nperchannel_norm = (e_hi**a - e_lo**a) / (emax**a - emin**a) case "Gaussian" : - raise NotImplementedError("Gaussian normalization for multiple bins not yet implemented") + mean, sdev, cutoff = params + emin = mean - cutoff * sdev + emax = mean + cutoff * sdev + + e_lo = np.minimum(emax, e_lo) + e_hi = np.minimum(emax, e_hi) + + e_lo = np.maximum(emin, e_lo) + e_hi = np.maximum(emin, e_hi) + + nperchannel_norm = (gauss_int(e_hi, mean, sdev) - gauss_int(e_lo, mean, sdev)) / (emax - emin) # If Nulambda is full-sky, its nbins will be 1, so division is a no-op. # We assume all FISBEL pixels have the same area. diff --git a/tests/response/test_response_conversion.py b/tests/response/test_response_conversion.py index 00fb8e0a..bed15533 100644 --- a/tests/response/test_response_conversion.py +++ b/tests/response/test_response_conversion.py @@ -161,8 +161,23 @@ def test_default_norms(tmp_path): # have their eff_area set to zero, not -inf or NaN assert all(np.isfinite(fdr.eff_area_correction)) - # cannot test Gaussian norm because it can only be used on - # a response with a single Ei bin + # Gaussian norm + c = RspConverter(bufsize = 100000, + norm = "Gaussian", + norm_params = [100, 100, 3]) + + c.convert_to_h5(rspgz_nonorm_response_path, + h5_filename = tmp_h5_filename, + overwrite = True) + + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "Gaussian 100.0 100.0 3.0" + + # bins with no defined spectral normalization should + # have their eff_area set to zero, not -inf or NaN + assert all(np.isfinite(fdr.eff_area_correction)) + def test_convert_h5_to_rsp(tmp_path): From dc747fbe4cf575495d8fcd73b6d8ecef5c7d1f7c Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Wed, 4 Feb 2026 16:46:19 -0600 Subject: [PATCH 8/8] * add single-Ei-bin test response to test Mono normalization with no energy * make sure we use an empty norm parameter list, not [None], to describe a Mono norm without energy internally, so that we don't emit 'None' in the saved response headers --- cosipy/response/RspConverter.py | 13 ++++++++--- ...full_detector_response_mono_no_norm.rsp.gz | Bin 0 -> 6515 bytes tests/response/test_response_conversion.py | 22 +++++++++++++++--- 3 files changed, 29 insertions(+), 6 deletions(-) create mode 100644 cosipy/test_data/test_full_detector_response_mono_no_norm.rsp.gz diff --git a/cosipy/response/RspConverter.py b/cosipy/response/RspConverter.py index 7775e4e7..1c86e82d 100644 --- a/cosipy/response/RspConverter.py +++ b/cosipy/response/RspConverter.py @@ -403,7 +403,7 @@ def _validate_norm_params(self, norm, params): if len(params) > 1: raise ValueError(f"Mono normalization takes a most one param; {len(params)} given") # emono if given, else dummy value - params = (None,) if len(params) == 0 else ( int(params[0]), ) + params = () if len(params) == 0 else ( int(params[0]), ) case 'Linear': if len(params) != 2: @@ -506,14 +506,18 @@ def gauss_int(x, mu, sigma): if not self.quiet: logger.info("normalization: mono") - emono = params[0] - if emono is None: + if params == (): if ei_axis.nbins > 1: raise ValueError("Cannot specify Mono norm without energy for Ei axis with multiple bins") else: # all energy is in the single bin nperchannel_norm = np.array([1.]) + + # fill in a header value so we don't write "None" + hdr["norm_params"] = (ei_axis.centers[0],) else: + emono = params[0] + # set just the Ei bin containing the mono energy to 1 nperchannel_norm = np.zeros(ei_axis.nbins) nperchannel_norm[(e_lo <= emono) & (e_hi >= emono)] = 1. @@ -538,9 +542,12 @@ def gauss_int(x, mu, sigma): case "Gaussian" : mean, sdev, cutoff = params + emin = mean - cutoff * sdev emax = mean + cutoff * sdev + logger.info(f"normalization: Gaussian with energy range [{emin}-{emax}]keV") + e_lo = np.minimum(emax, e_lo) e_hi = np.minimum(emax, e_hi) diff --git a/cosipy/test_data/test_full_detector_response_mono_no_norm.rsp.gz b/cosipy/test_data/test_full_detector_response_mono_no_norm.rsp.gz new file mode 100644 index 0000000000000000000000000000000000000000..11093d1ebfb83bd15d1d2cbfa00f63923acb887e GIT binary patch literal 6515 zcmV-(8I0y1iwFp@T7zi-19D|^aBpsNWnVF3X>MO`Z*FgLZ7yaCJTWq0Bt{{k7bFU>|9$HFYF|n=#6dvv+{@E_cJD{ks#UA% z{QCO*+poX-@!gv*uJ7M`@#E+I{NcNAzxm{&#QReD%fefB5#@_uqc^ z_VZW&^2znzufK|q{`v3Uy}AB$2KvM6ci+GL+aEvu?A`ZoUVr`Sqff72{ra09ufKlr z=F9*5Rs8#zUoi`VD>uRi+oRX@J|?URrG{XgFR?bSygU-$K=-+%hAe(d+(ynXlf z^;g%MZ{B?W<=?Nr{_mT=;Hd# z?APD@&)c6(>Ywhte{)^ejF(b*+3`|)yfiK}E_5z*E(|VAE=(>gE-Wr=E<9X#;zGGA zFL9xixKQfN3l|y}8W%bj1{VewCKo0b78myYGW8$c{V$$(` zd1BCY&3KJ_V#8*P*5lo_d#sXMYWy3Ix0DfgwQI!gwiW}#6UU0TpA{FJLbu3yNhuuJ5_8xknsNi7{jR z!9!!BF@qW-|8h$_AtCLmeD!Au-JHF%(OAlcmfk+9XG2TZcE(D z|KhQ6%q0ft+_?R<#9M1&%b2a!9<#48Ozd)4ISY%~?A){IEzQlzMT7~*VaHHm$hkmg z93PKavW)v0UM);K&MS_hGUy7!c$`!mEswo(|2wH-B)I|qbO-=|xpTPlpi2DzVz&jDU;fi$LGuNGuAlC#vF_49goU-oP;p=}hj|_X zFyPBGwkW+43@RAdE<#*k9&ulf<%BEZ7qN>C8;M0#2QGkGY@^fM@tEE?*6l7fWtKWiFY)!m29`CUY-~oqfX(Bb8poWtQXctSxPOqoM`tbg6 zn#q7$S1%mv<7~sPV{^7v8k-nayuh3A7wYTr=^{9==q<3fwAxtd;u8b3eSA(~>nDc9`6RS+7;G?6 z0LzMvFv`Yz8HUT!)xu}u(?HrInh<81%j_i)j33nz;!`ofI9(uUvN5qETG0tlTc@jD z68uEzxeW@j)40PFYuMKZIE&YqRNtaq2S>U)j&ZfbQ#PgRBDeq1y$b-uj)Hv9i|2$S zW`1e%j|+UX+{6{40&jyHpezrGJ$(|4hNLCY0;pLZgdwo(d4IjWe;;p}hEyO4sB7Zu z?XsXJ<>6FutqGjhp18O66yOB42l)pvW&2FEW?lk8GundIv9|)F^(FJA27)|;J-{G-e!|P^MIP=v12wH zwgw1`1O-5eql8MmkfRFiSm|0@-HBA@HSWXWfSmxAm=)kh7$b-TD} zt?^`e4o?~~D5nXXs5PVh;&#YZY>6MK0=C@fk@f<#$G$+)BFw^4b&eu8B8UJJtXMRH z1%(k2w_^KT=6xPv>R>pPmd0MRBsqCGOf3eF&kJCfeB3`hrcuPgsS^fTjXWdC)03>p zT`droj+riF7Wd&cLOOw+({uYDOtTq*I}oCI%i!@&J)&iHL2^@gu8_oYblW1uD^9h; zW`rlRm9x=E7FM zD$}Zf6hJi(N4L=cMH$R>%IYg7H8`q)m3Q5%L()tm8!l|B2KA4Jg+4j89|D+>|$wdxOQ@2?F45h zZZSF6UIeTBE15q@>eOM#Vkh1b>+9>qEBbSi*i6|vPkLIi5`YJui^)PqGNJMDkAN6( z_CWVwhooELum=7x7?Y4wJ!npanEzzz zM<~EIVg&elNiG2YtRBH-I?WUa3Df*?3RNtt9u!J9S{}a{6ui6{1`4)(OJNKA`W^}V zXH_N&B+J?~3>bQ8WDT4qyWI?2Ai^xHoli=*GmGe~MBM<^SXS(&c zsUC3!2t|5vaT(%9hKz0i5-d7VCjLLb8)kNAzPQb#HkTgB5dwuK4YR5>^T^BGe7cv_ z#HzDxrk~T*!h3pVkWd%$^&_0{atVn;RH+Al$_pv_utXHwas6}VR;JTq&zKta5bw^H z6pzKnlA=2g;v4M|~WnI2-ZU zEUYrTb3ccGC0ee)fewtw_eI~WL4}>4+`nLMWNib7?;|J4$a@jJCB8O*M7b(s#)T5jB=6NILlb!!Id zImLXqr5IFUox~2nX~CyAEbd4nSOcLkQ?)P!V(OhdXvq7)dUCvLW#Z%}Fy_9KSCCH= zI4YkBgGy*7hrUERFd2`U>7LhNYh~#81ZLHneyhh;GxaC4 z9?}GY-42*TEgB=EeU)nPH27B)`mxud|Db2aQ-?fxr&1NuSK~Z-<4n2EJs=DLQ$6W^ z^}d#%GI9i2%0{b|uo0uoyfX)hPFRqizHmjd95Xr`6TY)_0cLV5()np4GX`*@ z>g5|O!$cQV5<;qjw&TJ~h{c0UvKlBJs`xaV68K{fI5O|;ycbCpXP(3lb);I}4hw#>qk7z-f%gP@dh_s{=-{te3sN;M$umYA+od<4i< zO|_290%z;0z}Nkv^8x`6hu5ue=%z!}@k;z&K#k^u-PXf}OUg7UFmtFgH9ONtsxAd0 zWQjzGX%X2~Ixovd(ncZ@h#1@{rr8Wd)^IR1t5?R#C$k`eWL0`)YVnDc^Qszf7F6nb zFR8#f9;{r+3t1d39CyiB=rU% zK+4q9?4gJys|*0x!}NRNQx>&EaH!&)dvYt}DKoGEHx#X9QkHCsxd-~rv(S5Q%UWCQ zxu}hyNiHFsvzzxzH7wp=863&kz!-CXO^7Im>)d&Iz8|)NBmldS^OcM{MY|Gbs_8RN ziKv`8CMjZB+SSS>r3$@YDg_dpp$>j3Yh5;1#Fq?lschy?8Mb0DM2(q`#AnaAFT;5l z$}_U`$_70oMLkNyFvJbgO0oHPzWx zl_VNz;_?1d&=$wCJc0{6SulRG5Ysb-+F3oTX#iT;W6YhZIOBx^Rkfv-Ld?M_0O!N% z(hG!@e^8fnOTy9@61GSL}MdM79o2#?X`+%`QYiMD1cOvopw4g%=+eu z(oddozGVGJjlVCW0YQ8wTeZWWn*>IrLHjaE>PKu-+`!Ms$COry^`;7UYLp7*na|a% z-R|^EgqC~Zfsee~DiC+xRyOCcLuypWc&vCEH96>$&fThI06}{ZZ&Rq{QL$yiaWYZ0 z7-WHKgA*;*)Ch%PI=Ch>bT>;&4xj|$p3HU_Xuh)L9En8}o_b#+vYx9hHeIrA@%@(+=pX&@jD(<&Z_8Vq$c`>wF1crh(y zp}7CH&q8)=re7p&pU)P=jK?~1anB*LGR3z02Lg4N8ISzjgL5ifr3_bLW5jQ}fM zNz5Q9ycF~89 znBdjroWqT$NHh<#r5B}~>af8{64A8NB)-{4Bs~(xArHjJerH{;9|<9Mj^CgYCV*#% zYDW_(+`O~+eGEYwH2r4+856&uRkzjec{eS*W&&kM$0mp36;dh(3RG~a5Pa})8JG>^ zmKGHSJ)>csBNIe+HtTP8c~GF6r?}l2R3^71F2s+R9QCK zY&u2P-A*r`M?kKZ?LbqNY!Z?u4S)0alL~&VCR(XJ6TX-sGAJ~?r7X1*PiXXctXE%T z9+WOu%E#uENp^HBqC-oZKdRFfTvkFCNIb-omp5P99zTRi?Mk+HTU`5zC{bO=<{e9t z#%ebUF3PX;w3wc%a1B9S>u==D>jjyrrAOS%LMJ=?(BF`u;>iGXTIIIFpKUtQE~@_Q zzTQ@qYajW1Z7C0w%-UP#FtDh{>Y5o0&%BErF5_sjJ?pgW(#{@_v9dy5i&Y>90QP?? zU8q#R5KOg1)@GR;h+q-fvLSbGx#q-CW^Pf3jQ`1omE5;wM_kQ_b|1VUR#BmBXd!l2 zSrM+1Me88;bShBVX=?kQT?Wg-ZO@Evh>g=!1x=}e{Eu)o?PMJH33g2n1zjQbT{5VM zG2S;jWb2s8v_W9G%NW<49gPGO1C0trtZkewLxv51;5XoUXNf3bX=eT242KMN9L-7r z1P-Q$#v~zu+Y_s}Z{kusuVvjcHvLgLKn#K`fJ_CuMZ8IwvAgrPIhSnkQMn9d<5BP(cm$QiO;J-mTXUy zaU{xWteCW=8-pK@K^u{zJn3ilP7J$yyN!X^MBYt4xDS)+-n010i8cBvvU%)NT}`UVUdf*-nLw zUem{-PgZ?X09EpN36Y9wgMr&Vo07t4cEoK(f@+hr%V%OR$0HUXLxMy*>umH8`n$S@ z5JLCiI3j13=Ha=6EDO=BnmU-aW(OW_WBWTy#5vOkF+VT<}>QI6_e%03hxvR0Dv+_t>Z zMi%5~Crn8qJD6I^5+<5d{cxqxRmzf33Ft~QHs3>G88+VQe2>#`!}Ivdies zhCeD$qL*ox^)~j#_>i(k#`&6c=0^@@;soZ5agw%W35Ev`!O1M1x$b2bYoUr$*7?T5 zJWWqq$+!%qSQHgVsy%BnyD@9#9Lq5Yd)rE5qd^a@nB2$N_QO4TjtLZP6#}p+Egxta zNB6xYIfdr6OhV%{yWbKp=GT(#%VuaycBHw|Q%)~#DmmliUP_`}8G*jmEcN#Hf)Kjn zqwkLLB5~4sYhIU_W9MyC3dp9#e{4Yd5F};CQ~lZ3NL0eJtNGTN3|6O23OKDcW6;hm zWnac>+0C$1gmN6jZwhHW1-+9I%L0ZS*0$_m0!7wU{F21*iP;EyXFo|&`f6LsofR!C z|C|QY!YbWM^O%*Kl2x^lwJ6$fhzhf$Qn5_S)*)(C$B9Se^$a7wx9|!y8%S4(Jxv%_w0`qA(M|&y9>7TRL*t*T_IZabC!X#g=!M}QDZysCK*itVH zx-@gkk6$HDv#)mXPR|qp)u3F+8#|s&Vpx_X_L)*|&G_0fA#OISb$XZY`y@NWTpPPn zWcE8Z$;s$Fg+@$!I)d4Z7q3q~+pQvEe?9##VCzwUGwUlQHJALBL4Gwu81w5CGIFIrmA%FSy2rxq$G4d-Y+3>gy Z3*zFd&)$9e-T7YP{{ib1`{XS-002?Evy%V- literal 0 HcmV?d00001 diff --git a/tests/response/test_response_conversion.py b/tests/response/test_response_conversion.py index bed15533..a95780f4 100644 --- a/tests/response/test_response_conversion.py +++ b/tests/response/test_response_conversion.py @@ -10,6 +10,8 @@ rspgz_nonorm_response_path = test_data.path / "test_full_detector_response_no_norm.rsp.gz" +rspgz_mono_nonorm_response_path = test_data.path / "test_full_detector_response_mono_no_norm.rsp.gz" + h5_response_path = test_data.path / "test_full_detector_response.h5" @@ -107,6 +109,23 @@ def test_default_norms(tmp_path): # have their eff_area set to zero, not -inf or NaN assert all(np.isfinite(fdr.eff_area_correction)) + # mono norm with no energy should work for single-bin response + c = RspConverter(bufsize = 100000, + norm = "Mono", + norm_params = []) + + c.convert_to_h5(rspgz_mono_nonorm_response_path, + h5_filename = tmp_h5_filename, + overwrite = True) + + with FullDetectorResponse.open(tmp_h5_filename) as fdr: + norm_info = fdr.headers["SP"] + assert norm_info == "Mono" + + # bins with no defined spectral normalization should + # have their eff_area set to zero, not -inf or NaN + assert all(np.isfinite(fdr.eff_area_correction)) + # if no monoenergetic energy is given, Mono is # undefined with multiple Ei bins with raises(ValueError): @@ -117,9 +136,6 @@ def test_default_norms(tmp_path): h5_filename = tmp_h5_filename, overwrite = True) - # cannot test Mono with no param without an .rsp file - # containing only a single Ei axis bin - c = RspConverter(bufsize = 100000, norm = "powerlaw", norm_params = [50, 10000, 0.9])