From 2b838bb452e8cd5be50990346b29b4796e458fd9 Mon Sep 17 00:00:00 2001 From: Julian Saffer <31139467+jsaffer@users.noreply.github.com> Date: Thu, 12 Dec 2024 19:43:21 +0100 Subject: [PATCH 01/10] Adding flux weighting option to the effective_area function of Weighter objects MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Instead of weighting to a 1 GeV⁻¹ m⁻² sr⁻¹ flux, one can now chose to use any of the cosmic ray flux models as well in order take into account the spectral shape and varying mass composition inside energy bins. The scalar flux is still possible and the default (1e-4). --- src/simweights/_weighter.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index 2896109..5112039 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -139,6 +139,7 @@ def effective_area( self: Weighter, energy_bins: ArrayLike, cos_zenith_bins: ArrayLike, + flux: Any = 1e-4, # default is 1 GeV^-1 m^-2 sr^-1 flux mask: ArrayLike | None = None, ) -> NDArray[np.float64]: r"""Calculate The effective area for the given energy and zenith bins. @@ -182,7 +183,7 @@ def effective_area( energy = self.get_weight_column("energy") cos_zen = self.get_weight_column("cos_zen") - weights = self.get_weights(1e-4) + weights = self.get_weights(flux) maska = np.full(weights.size, 1, dtype=bool) if mask is None else np.asarray(mask, dtype=bool) assert maska.shape == weights.shape @@ -198,8 +199,18 @@ def effective_area( assert np.array_equal(enbin, energy_bins) assert np.array_equal(czbin, cos_zenith_bins) - e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin)) - return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) + if np.isscalar(flux): + e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin)) + return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) + elif callable(flux): + flux_pdgids = [pdgid.value for pdgid in flux.pdgids] + flux_func = lambda E: sum([flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](E) for i_species in range(nspecies)]) + from scipy.integrate import quad + flux_integrals = np.asarray([quad(flux_func, energy_bins[bin_index], energy_bins[bin_index+1])[0] for bin_index in range(len(energy_bins)-1)]) + e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) + return np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) + else: + raise ValueError("only scalar flux or cosmic ray flux models are supported right now") def __add__(self: Weighter, other: Weighter | int) -> Weighter: if other == 0: From d122cca3ed585b4845567c8e2a345e722b7d846d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 12 Dec 2024 18:50:59 +0000 Subject: [PATCH 02/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/simweights/_weighter.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index 5112039..afd7996 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -139,7 +139,7 @@ def effective_area( self: Weighter, energy_bins: ArrayLike, cos_zenith_bins: ArrayLike, - flux: Any = 1e-4, # default is 1 GeV^-1 m^-2 sr^-1 flux + flux: Any = 1e-4, # default is 1 GeV^-1 m^-2 sr^-1 flux mask: ArrayLike | None = None, ) -> NDArray[np.float64]: r"""Calculate The effective area for the given energy and zenith bins. @@ -204,9 +204,20 @@ def effective_area( return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) elif callable(flux): flux_pdgids = [pdgid.value for pdgid in flux.pdgids] - flux_func = lambda E: sum([flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](E) for i_species in range(nspecies)]) + flux_func = lambda E: sum( + [ + flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](E) + for i_species in range(nspecies) + ] + ) from scipy.integrate import quad - flux_integrals = np.asarray([quad(flux_func, energy_bins[bin_index], energy_bins[bin_index+1])[0] for bin_index in range(len(energy_bins)-1)]) + + flux_integrals = np.asarray( + [ + quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0] + for bin_index in range(len(energy_bins) - 1) + ] + ) e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) return np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) else: From 8c41dcdf225eff871cfab5dfa8c7090603790beb Mon Sep 17 00:00:00 2001 From: Julian Saffer <31139467+jsaffer@users.noreply.github.com> Date: Fri, 13 Dec 2024 10:20:46 +0100 Subject: [PATCH 03/10] Changed order of arguments of the effective_area function The new optional flux argument is now the last in the list of arguments. This way, confusion with the mask argument, which remains in 3rd place, is avoided. --- src/simweights/_weighter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index afd7996..6ee014e 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -139,8 +139,8 @@ def effective_area( self: Weighter, energy_bins: ArrayLike, cos_zenith_bins: ArrayLike, - flux: Any = 1e-4, # default is 1 GeV^-1 m^-2 sr^-1 flux mask: ArrayLike | None = None, + flux: Any = 1e-4, # default is 1 GeV^-1 m^-2 sr^-1 flux ) -> NDArray[np.float64]: r"""Calculate The effective area for the given energy and zenith bins. From 4121c11b74eec507bd399219f59454429d2781bf Mon Sep 17 00:00:00 2001 From: Julian Saffer <31139467+jsaffer@users.noreply.github.com> Date: Fri, 13 Dec 2024 11:53:37 +0100 Subject: [PATCH 04/10] Addrerssing some of the warnings and suggestions from pylint and ruff This should fix C3001 (lambda expression replaced by def function definition), R1728 (using generator for sum of flux components), D417 (added description for flux argument), E731 (lambda -> def), N803 (E -> energy), TRY003 and EM101 (assigned error message to variable), also I changed the ValueError to a TypeError which should be more appropriate --- src/simweights/_weighter.py | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index 6ee014e..d0af65f 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -158,13 +158,16 @@ def effective_area( Args: - energy_bins : array_like + energy_bins: array_like A length N+1 array of energy bin edges - cos_zenith_bins : array_like + cos_zenith_bins: array_like A length M+1 array of cos(zenith) bin edges mask: array_like boolean array where 1 indicates to use the event in the calculation. Must have the same length as the simulation sample. + flux: Any + A model describing the flux of the primaries to weight against. For + possible types, see get_weights() Returns: array_like @@ -204,24 +207,15 @@ def effective_area( return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) elif callable(flux): flux_pdgids = [pdgid.value for pdgid in flux.pdgids] - flux_func = lambda E: sum( - [ - flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](E) - for i_species in range(nspecies) - ] - ) + def flux_func(energy): + return sum(flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) for i_species in range(nspecies)) from scipy.integrate import quad - - flux_integrals = np.asarray( - [ - quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0] - for bin_index in range(len(energy_bins) - 1) - ] - ) + flux_integrals = np.asarray([quad(flux_func, energy_bins[bin_index], energy_bins[bin_index+1])[0] for bin_index in range(len(energy_bins)-1)]) e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) return np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) else: - raise ValueError("only scalar flux or cosmic ray flux models are supported right now") + mesg = f"flux of type {type(flux)} is supplied but only scalar flux or cosmic ray flux models are supported so far" + raise TypeError(mesg) def __add__(self: Weighter, other: Weighter | int) -> Weighter: if other == 0: From e05dee0839a94d18b0a33c2aeceb2b6a025b2605 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 13 Dec 2024 10:57:51 +0000 Subject: [PATCH 05/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/simweights/_weighter.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index d0af65f..8a7028c 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -207,10 +207,21 @@ def effective_area( return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) elif callable(flux): flux_pdgids = [pdgid.value for pdgid in flux.pdgids] + def flux_func(energy): - return sum(flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) for i_species in range(nspecies)) + return sum( + flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) + for i_species in range(nspecies) + ) + from scipy.integrate import quad - flux_integrals = np.asarray([quad(flux_func, energy_bins[bin_index], energy_bins[bin_index+1])[0] for bin_index in range(len(energy_bins)-1)]) + + flux_integrals = np.asarray( + [ + quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0] + for bin_index in range(len(energy_bins) - 1) + ] + ) e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) return np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) else: From 09eb1a1418eb00e5577fe14904c09dda1a834c33 Mon Sep 17 00:00:00 2001 From: Julian Saffer <31139467+jsaffer@users.noreply.github.com> Date: Fri, 13 Dec 2024 14:51:45 +0100 Subject: [PATCH 06/10] More fixes of pylint / ruff complaints Adding type annotations for flux_func function, and moving the return outside of the it/elif/else condition. --- src/simweights/_weighter.py | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index 8a7028c..e17322f 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -204,29 +204,19 @@ def effective_area( assert np.array_equal(czbin, cos_zenith_bins) if np.isscalar(flux): e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin)) - return np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) + effective_area_binned = np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) elif callable(flux): flux_pdgids = [pdgid.value for pdgid in flux.pdgids] - - def flux_func(energy): - return sum( - flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) - for i_species in range(nspecies) - ) - + def flux_func(energy: ArrayLike) -> NDArray[np.float64]: + return sum(flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) for i_species in range(nspecies)) from scipy.integrate import quad - - flux_integrals = np.asarray( - [ - quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0] - for bin_index in range(len(energy_bins) - 1) - ] - ) + flux_integrals = np.asarray([quad(flux_func, energy_bins[bin_index], energy_bins[bin_index+1])[0] for bin_index in range(len(energy_bins)-1)]) e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) - return np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) + effective_area_binned = np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) else: mesg = f"flux of type {type(flux)} is supplied but only scalar flux or cosmic ray flux models are supported so far" raise TypeError(mesg) + return effective_area_binned def __add__(self: Weighter, other: Weighter | int) -> Weighter: if other == 0: From 68f0c5e582a68dc5e7d181d9e7112a4990f62bdb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 13 Dec 2024 13:52:14 +0000 Subject: [PATCH 07/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/simweights/_weighter.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index e17322f..aa73d07 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -207,10 +207,21 @@ def effective_area( effective_area_binned = np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) elif callable(flux): flux_pdgids = [pdgid.value for pdgid in flux.pdgids] + def flux_func(energy: ArrayLike) -> NDArray[np.float64]: - return sum(flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) for i_species in range(nspecies)) + return sum( + flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) + for i_species in range(nspecies) + ) + from scipy.integrate import quad - flux_integrals = np.asarray([quad(flux_func, energy_bins[bin_index], energy_bins[bin_index+1])[0] for bin_index in range(len(energy_bins)-1)]) + + flux_integrals = np.asarray( + [ + quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0] + for bin_index in range(len(energy_bins) - 1) + ] + ) e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) effective_area_binned = np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) else: From 0056dd0059b23e16ce6fd43e157b8fff25b31d4d Mon Sep 17 00:00:00 2001 From: Kevin Meagher <11620178+kjmeagher@users.noreply.github.com> Date: Mon, 16 Dec 2024 18:20:09 -0600 Subject: [PATCH 08/10] clean up flux weighting for effective_area --- pyproject.toml | 2 +- src/simweights/_weighter.py | 28 +++++++++------------------- 2 files changed, 10 insertions(+), 20 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 374d6de..a6d23e6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -78,7 +78,7 @@ warn_unreachable = true max-line-length = "128" [tool.pylint.messages_control] -disable = "C0114,R0902,R0913,R0917" +disable = "C0114,R0902,R0913,R0917,R0914" [tool.pytest.ini_options] addopts = ["-ra", "--strict-config", "--strict-markers", "--cov=simweights", "-W ignore"] diff --git a/src/simweights/_weighter.py b/src/simweights/_weighter.py index aa73d07..df44892 100644 --- a/src/simweights/_weighter.py +++ b/src/simweights/_weighter.py @@ -9,6 +9,7 @@ from typing import TYPE_CHECKING, Any, Callable, Iterable import numpy as np +from scipy.integrate import quad # pylint: disable=import-error from ._utils import get_column, get_table @@ -140,7 +141,7 @@ def effective_area( energy_bins: ArrayLike, cos_zenith_bins: ArrayLike, mask: ArrayLike | None = None, - flux: Any = 1e-4, # default is 1 GeV^-1 m^-2 sr^-1 flux + flux: Any = 1, # default is 1 GeV^-1 cm^-2 sr^-1 flux ) -> NDArray[np.float64]: r"""Calculate The effective area for the given energy and zenith bins. @@ -175,6 +176,7 @@ def effective_area( is the number of cos(zenith) bins. """ + cm2_to_m2 = 1e-4 energy_bins = np.array(energy_bins) cos_zenith_bins = np.array(cos_zenith_bins) @@ -198,36 +200,24 @@ def effective_area( bins=[cos_zenith_bins, energy_bins], ) - nspecies = len(np.unique(self.get_weight_column("pdgid")[maska])) - assert np.array_equal(enbin, energy_bins) assert np.array_equal(czbin, cos_zenith_bins) + pdgids = np.unique(self.get_weight_column("pdgid")[maska]) + if np.isscalar(flux): - e_width, z_width = np.meshgrid(np.ediff1d(enbin), np.ediff1d(czbin)) - effective_area_binned = np.asarray(hist_val / (e_width * 2 * np.pi * z_width * nspecies), dtype=np.float64) + flux_integrals = len(pdgids) * np.ediff1d(enbin) elif callable(flux): - flux_pdgids = [pdgid.value for pdgid in flux.pdgids] - - def flux_func(energy: ArrayLike) -> NDArray[np.float64]: - return sum( - flux._funcs[flux_pdgids.index(np.unique(self.get_weight_column("pdgid")[maska])[i_species])](energy) - for i_species in range(nspecies) - ) - - from scipy.integrate import quad - flux_integrals = np.asarray( [ - quad(flux_func, energy_bins[bin_index], energy_bins[bin_index + 1])[0] + sum(quad(flux, energy_bins[bin_index], energy_bins[bin_index + 1], args=(p,))[0] for p in pdgids) for bin_index in range(len(energy_bins) - 1) ] ) - e_width, z_width = np.meshgrid(flux_integrals, np.ediff1d(czbin)) - effective_area_binned = np.asarray(1e-4 * hist_val / (e_width * 2 * np.pi * z_width), dtype=np.float64) else: mesg = f"flux of type {type(flux)} is supplied but only scalar flux or cosmic ray flux models are supported so far" raise TypeError(mesg) - return effective_area_binned + e_int, z_int = np.meshgrid(flux_integrals, np.ediff1d(czbin)) + return np.asarray(cm2_to_m2 * hist_val / (e_int * 2 * np.pi * z_int), dtype=np.float64) def __add__(self: Weighter, other: Weighter | int) -> Weighter: if other == 0: From 3f99af3a5346838dbe4991797645bf44c2e674ee Mon Sep 17 00:00:00 2001 From: Julian Saffer Date: Wed, 18 Dec 2024 04:10:40 -0600 Subject: [PATCH 09/10] Adding tests for weighted effective area, also covering the ValueError of invalid flux argument and the TypeError of a valid flux argument that is not supported yet --- tests/test_weighter.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/test_weighter.py b/tests/test_weighter.py index 4433188..9f4457a 100755 --- a/tests/test_weighter.py +++ b/tests/test_weighter.py @@ -207,6 +207,23 @@ def test_effective_area(self): self.weighter1.effective_area(np.linspace(5e5, 5e6, self.N1 + 1), [0, 1]), [np.full(self.N1, self.c1.etendue / 2e4 / np.pi)], ) + self.assertAlmostEqual( + self.weighter1.effective_area([5e5, 5e6], [0, 1], flux=TIG1996())[0][0], + 149998.97822505102, + 6, + ) + self.assertAlmostEqual( + self.weighter1.effective_area([5e5, 5e6], [0, 1], np.ones(self.N1, dtype=bool), flux=TIG1996())[0][0], + 149998.97822505102, + 6, + ) + + with self.assertRaises(ValueError): + self.weighter1.effective_area([5e5, 5e6], [0, 1], flux="flux") + + with self.assertRaises(TypeError): + self.weighter1.effective_area([5e5, 5e6], [0, 1], flux=list(range(self.N1))) + def test_weighter_addition(self): weighter_sum = self.weighter1 + self.weighter1 From 4ad2c984a3fb4167b224e7aefcbf32c822318438 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 18 Dec 2024 10:13:45 +0000 Subject: [PATCH 10/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_weighter.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_weighter.py b/tests/test_weighter.py index 9f4457a..3987313 100755 --- a/tests/test_weighter.py +++ b/tests/test_weighter.py @@ -223,7 +223,6 @@ def test_effective_area(self): with self.assertRaises(TypeError): self.weighter1.effective_area([5e5, 5e6], [0, 1], flux=list(range(self.N1))) - def test_weighter_addition(self): weighter_sum = self.weighter1 + self.weighter1