From 7abc89aa1a3d96a9d2ebb1619234d1a9b17e622a Mon Sep 17 00:00:00 2001 From: Julian Saffer <31139467+jsaffer@users.noreply.github.com> Date: Wed, 21 May 2025 17:25:28 +0200 Subject: [PATCH 1/4] Correct the 4-component parameterization of GST The original GST paper (http://doi.org/10.1007/s11467-013-0319-7) gives a parameterization of the CR flux in several populations based on p, He, C, O and Fe mass groups (Table 3). I see that this doesn't match the 5-component model that is used in some working groups (that's why I didn't touch the 5-component model) but for the 4-component version it would be appropriate to use the correct Z values and combine C and O. --- src/simweights/_fluxes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py index 974125a..537dfc9 100644 --- a/src/simweights/_fluxes.py +++ b/src/simweights/_fluxes.py @@ -327,7 +327,7 @@ class GlobalFitGST(CosmicRayFlux): class GlobalFitGST_IT(CosmicRayFlux): # pylint: disable=invalid-name r"""GlobalFitGST for four components [p, He, O, Fe]. - The Oxygen group is the sum of Nitrogen and Aluminum groups of GlobalFitGST. + The Oxygen group is the sum of Carbon and Oxygen groups. """ pdgids = PDGID_4COMP From 08955363ec3962082802cb1ca80a147eca89f25c Mon Sep 17 00:00:00 2001 From: Julian Saffer <31139467+jsaffer@users.noreply.github.com> Date: Thu, 22 May 2025 23:16:52 +0200 Subject: [PATCH 2/4] Update GST_IT Updating the 4-component GST model to the 4-population version, having C and O combined with the correct Z and Fe combined with the super heavy component --- src/simweights/_fluxes.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py index 537dfc9..c1e3761 100644 --- a/src/simweights/_fluxes.py +++ b/src/simweights/_fluxes.py @@ -325,9 +325,10 @@ class GlobalFitGST(CosmicRayFlux): class GlobalFitGST_IT(CosmicRayFlux): # pylint: disable=invalid-name - r"""GlobalFitGST for four components [p, He, O, Fe]. + r"""GlobalFitGST (4 populations) for four components [p, He, O, Fe]. The Oxygen group is the sum of Carbon and Oxygen groups. + The Iron group is the sum of Iron, Tellurium and Mercury groups. """ pdgids = PDGID_4COMP From e5c7c05dde5ceaea8475d4652fa52c676c954354 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 23 May 2025 20:24:49 +0000 Subject: [PATCH 3/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/simweights/_fluxes.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py index c1e3761..fbcf9dd 100644 --- a/src/simweights/_fluxes.py +++ b/src/simweights/_fluxes.py @@ -333,19 +333,20 @@ class GlobalFitGST_IT(CosmicRayFlux): # pylint: disable=invalid-name pdgids = PDGID_4COMP _funcs = ( - lambda E: 0.7 * E**-2.66 * exp(-E / 1.2e5) + 0.015 * E**-2.4 * exp(-E / 4e6) + 0.0014 * E**-2.4 * exp(-E / 1.3e9), - lambda E: 0.32 * E**-2.58 * exp(-E / 1.2e5 / 2) + 0.0065 * E**-2.3 * exp(-E / 4e6 / 2), - lambda E: ( - 0.01 * E**-2.40 * exp(-E / 1.2e5 / 7) - + 0.0006 * E**-2.3 * exp(-E / 4e6 / 7) - + 0.013 * E**-2.40 * exp(-E / 1.2e5 / 13) - + 0.0007 * E**-2.3 * exp(-E / 4e6 / 13) - ), - lambda E: ( - 0.006 * E**-2.30 * exp(-E / 1.2e5 / 26) - + 0.00023 * E**-2.2 * exp(-E / 4e6 / 26) - + 0.0000025 * E**-2.2 * exp(-E / 1.3e9 / 26) - ), + lambda E: 0.7000 * E**-2.66 * exp(-E / 1.2e5) + + 0.0150 * E**-2.4 * exp(-E / 4e6) + + 0.0012 * E**-2.4 * exp(-E / 1.5e9) + + 0.00012 * E**-2.4 * exp(-E / 40e9), + lambda E: 0.3200 * E**-2.58 * exp(-E / 1.2e5 / 2) + 0.0065 * E**-2.3 * exp(-E / 4e6 / 2), + lambda E: 0.0100 * E**-2.40 * exp(-E / 1.2e5 / 6) + + 0.0006 * E**-2.3 * exp(-E / 4e6 / 6) + + 0.0130 * E**-2.40 * exp(-E / 1.2e5 / 8) + + 0.0007 * E**-2.3 * exp(-E / 4e6 / 8), + lambda E: 0.0060 * E**-2.30 * exp(-E / 1.2e5 / 26) + + 0.00021 * E**-2.2 * exp(-E / 4e6 / 26) + + 0.0000011 * E**-2.2 * exp(-E / 1.5e9 / 26) + + 0.00001 * E**-2.2 * exp(-E / 4e6 / 52) + + 0.000053 * E**-2.2 * exp(-E / 4e6 / 80), ) From 138a5db95146a943e7f24e7884ef685852d838b0 Mon Sep 17 00:00:00 2001 From: Kevin Meagher <11620178+kjmeagher@users.noreply.github.com> Date: Wed, 29 Apr 2026 14:23:45 -0500 Subject: [PATCH 4/4] Rename new flux model `GlobalFitGST4Comp` to keep _IT model the same --- src/simweights/__init__.py | 2 + src/simweights/_fluxes.py | 61 +++++++++++++----- tests/test_fluxes.py | 7 +- .../test_flux_model_GlobalFitGST4Comp_.npz | Bin 0 -> 1079 bytes 4 files changed, 49 insertions(+), 21 deletions(-) create mode 100644 tests/test_fluxes/test_flux_model_GlobalFitGST4Comp_.npz diff --git a/src/simweights/__init__.py b/src/simweights/__init__.py index 261c920..cb6a56c 100644 --- a/src/simweights/__init__.py +++ b/src/simweights/__init__.py @@ -28,6 +28,7 @@ "GenieSurface", "GenieWeighter", "GlobalFitGST", + "GlobalFitGST4Comp", "GlobalFitGST_IT", "GlobalSplineFit", "GlobalSplineFit5Comp", @@ -57,6 +58,7 @@ GaisserH4a_IT, GaisserHillas, GlobalFitGST, + GlobalFitGST4Comp, GlobalFitGST_IT, GlobalSplineFit, GlobalSplineFit5Comp, diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py index fbcf9dd..8ae67c1 100644 --- a/src/simweights/_fluxes.py +++ b/src/simweights/_fluxes.py @@ -87,7 +87,7 @@ def _condition( def __call__(self: CosmicRayFlux, energy: ArrayLike, pdgid: ArrayLike) -> NDArray[float64]: energy_arr, pdgid_arr = broadcast_arrays(energy, pdgid) pcond = self._condition(energy_arr, pdgid_arr) - return piecewise(energy, pcond, self._funcs) + return piecewise(energy, pcond, self._funcs) # type: ignore[no-any-return,call-overload] class Hoerandel(CosmicRayFlux): @@ -325,28 +325,59 @@ class GlobalFitGST(CosmicRayFlux): class GlobalFitGST_IT(CosmicRayFlux): # pylint: disable=invalid-name - r"""GlobalFitGST (4 populations) for four components [p, He, O, Fe]. + r"""GlobalFitGST for four components [p, He, O, Fe]. + + The Oxygen group is the sum of Nitrogen and Aluminum groups of GlobalFitGST. + """ + + pdgids = PDGID_4COMP + _funcs = ( + lambda E: 0.7 * E**-2.66 * exp(-E / 1.2e5) + 0.015 * E**-2.4 * exp(-E / 4e6) + 0.0014 * E**-2.4 * exp(-E / 1.3e9), + lambda E: 0.32 * E**-2.58 * exp(-E / 1.2e5 / 2) + 0.0065 * E**-2.3 * exp(-E / 4e6 / 2), + lambda E: ( + 0.01 * E**-2.40 * exp(-E / 1.2e5 / 7) + + 0.0006 * E**-2.3 * exp(-E / 4e6 / 7) + + 0.013 * E**-2.40 * exp(-E / 1.2e5 / 13) + + 0.0007 * E**-2.3 * exp(-E / 4e6 / 13) + ), + lambda E: ( + 0.006 * E**-2.30 * exp(-E / 1.2e5 / 26) + + 0.00023 * E**-2.2 * exp(-E / 4e6 / 26) + + 0.0000025 * E**-2.2 * exp(-E / 1.3e9 / 26) + ), + ) + + +class GlobalFitGST4Comp(CosmicRayFlux): # pylint: disable=invalid-name + r"""GlobalFitGST4Comp (4 populations) for four components [p, He, O, Fe]. The Oxygen group is the sum of Carbon and Oxygen groups. The Iron group is the sum of Iron, Tellurium and Mercury groups. + Provided by J. Saffer. """ pdgids = PDGID_4COMP _funcs = ( - lambda E: 0.7000 * E**-2.66 * exp(-E / 1.2e5) - + 0.0150 * E**-2.4 * exp(-E / 4e6) - + 0.0012 * E**-2.4 * exp(-E / 1.5e9) - + 0.00012 * E**-2.4 * exp(-E / 40e9), + lambda E: ( + 0.7000 * E**-2.66 * exp(-E / 1.2e5) + + 0.0150 * E**-2.4 * exp(-E / 4e6) + + 0.0012 * E**-2.4 * exp(-E / 1.5e9) + + 0.00012 * E**-2.4 * exp(-E / 40e9) + ), lambda E: 0.3200 * E**-2.58 * exp(-E / 1.2e5 / 2) + 0.0065 * E**-2.3 * exp(-E / 4e6 / 2), - lambda E: 0.0100 * E**-2.40 * exp(-E / 1.2e5 / 6) - + 0.0006 * E**-2.3 * exp(-E / 4e6 / 6) - + 0.0130 * E**-2.40 * exp(-E / 1.2e5 / 8) - + 0.0007 * E**-2.3 * exp(-E / 4e6 / 8), - lambda E: 0.0060 * E**-2.30 * exp(-E / 1.2e5 / 26) - + 0.00021 * E**-2.2 * exp(-E / 4e6 / 26) - + 0.0000011 * E**-2.2 * exp(-E / 1.5e9 / 26) - + 0.00001 * E**-2.2 * exp(-E / 4e6 / 52) - + 0.000053 * E**-2.2 * exp(-E / 4e6 / 80), + lambda E: ( + 0.0100 * E**-2.40 * exp(-E / 1.2e5 / 6) + + 0.0006 * E**-2.3 * exp(-E / 4e6 / 6) + + 0.0130 * E**-2.40 * exp(-E / 1.2e5 / 8) + + 0.0007 * E**-2.3 * exp(-E / 4e6 / 8) + ), + lambda E: ( + 0.0060 * E**-2.30 * exp(-E / 1.2e5 / 26) + + 0.00021 * E**-2.2 * exp(-E / 4e6 / 26) + + 0.0000011 * E**-2.2 * exp(-E / 1.5e9 / 26) + + 0.00001 * E**-2.2 * exp(-E / 4e6 / 52) + + 0.000053 * E**-2.2 * exp(-E / 4e6 / 80) + ), ) diff --git a/tests/test_fluxes.py b/tests/test_fluxes.py index aab813a..e3d2757 100755 --- a/tests/test_fluxes.py +++ b/tests/test_fluxes.py @@ -29,6 +29,7 @@ simweights.TIG1996(), simweights.GlobalFitGST(), simweights.GlobalFitGST_IT(), + simweights.GlobalFitGST4Comp(), simweights.GlobalSplineFit(), simweights.GlobalSplineFit5Comp(), simweights.GlobalSplineFit_IT(), @@ -39,12 +40,6 @@ @pytest.mark.parametrize("flux", flux_models, ids=[x.__class__.__name__ for x in flux_models]) def test_flux_model(flux, ndarrays_regression): - # this is the old regression test it can stick around for a bit but will be deleted at a certain point - for pdgid in flux.pdgids: - v1 = flux(E, pdgid) - v2 = np.array(flux_values[flux.__class__.__name__][str(int(pdgid))]) / 1e4 - assert v1 == pytest.approx(v2, rel=1e-13) - ndarrays_regression.check({pdgid.name: flux(E, pdgid) for pdgid in flux.pdgids}, default_tolerance={"rtol": 1e-13}) # make sure you get zero for non CR primaries assert flux(E, 22) == pytest.approx(0) diff --git a/tests/test_fluxes/test_flux_model_GlobalFitGST4Comp_.npz b/tests/test_fluxes/test_flux_model_GlobalFitGST4Comp_.npz new file mode 100644 index 0000000000000000000000000000000000000000..119df0981d5ceccfe4a6cb4ba5e091fbe19fda5c GIT binary patch literal 1079 zcmWIWW@gc4U|`??Vnv36lzWr@Ljfm)2tzlr9 zVfraG&1#kv*H$hOhIK2?9hG^>_xz3C{>M)*?|zc;{J8#m&0WdexyN3s`|1A@&5f42 z?peG4Vs`P{;^!)dFP~p~dy?LoTU!+G?CP;s^@$^}N&j+BffRZTXFJO%TD`Ql<@;-Xw%Rd^Frazmf4yTVvRC{K z&2V^S0``v+KzZKrw zd$UV!%cb_OzjjeOCGGl7>FC>krv5$~w>UFVsb*VITE+$$H2-XUa3KxmAAVrexTTun z@(>l1ntr&v?xveJvcA8HdGJj1pH0~&?J}0RxpPu}M|z(qyBk<|d&NF4#ld(KHD!{1N09glP)uAwgEX6ltVb7aR8zk z7(paF6QOGZrBP7o18L&}sf2