Cr flux 4comp fix (#28)

* fixing 4-comp fluxes for GST and GSF * fixing 4-comp fluxes for GST and GSF * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * progress on #28 * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * progress on #28 * progress on #28 * progress on #28 * progress on #28 * progress on #28 * progress on #28 * progress on #28 * progress on #28 --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
icecube · Jan 25, 2024 · 01198a9 · 01198a9
1 parent d0f7fda
commit 01198a9
Show file tree

Hide file tree

Showing 4 changed files with 88 additions and 0 deletions.
diff --git a/src/simweights/__init__.py b/src/simweights/__init__.py
@@ -27,7 +27,9 @@
     "GaisserH4a_IT",
     "GaisserHillas",
     "GlobalFitGST",
+    "GlobalFitGST_IT",
     "GlobalSplineFit5Comp",
+    "GlobalSplineFit_IT",
     "Hoerandel",
     "Hoerandel5",
     "Hoerandel_IT",
@@ -51,7 +53,9 @@
     GaisserH4a_IT,
     GaisserHillas,
     GlobalFitGST,
+    GlobalFitGST_IT,
     GlobalSplineFit5Comp,
+    GlobalSplineFit_IT,
     Hoerandel,
     Hoerandel5,
     Hoerandel_IT,

diff --git a/src/simweights/_fluxes.py b/src/simweights/_fluxes.py
@@ -310,6 +310,26 @@ 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.
+    """
+
+    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 GlobalSplineFitBase(CosmicRayFlux):
     r"""Data-driven spline fit of the cosmic ray spectrum by Dembinski et. al. \ [#GSFDembinski].
 
@@ -352,6 +372,20 @@ def __init__(self: GlobalSplineFit5Comp) -> None:
         super().__init__()
 
 
+class GlobalSplineFit_IT(GlobalSplineFitBase):  # pylint: disable=invalid-name
+    r"""Sum of the flux of the GSF model for the standard 4 components injected by IceCube.
+
+    [(H),  (He),  (Li, Be, B, C, N, O, F, Ne),  (Na, Mg, Al, Si, P, S, Cl, Ar, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni)]
+    GSF is a Data-driven spline fit of the cosmic ray spectrum by Dembinski et. al. \ [#GSFDembinski].
+    """
+
+    pdgids = PDGID_4COMP
+    groups = ((1, 1), (2, 2), (3, 10), (11, 28))
+
+    def __init__(self: GlobalSplineFit_IT) -> None:
+        super().__init__()
+
+
 class FixedFractionFlux(CosmicRayFlux):
     """Total energy per particle flux flux split among mass groups with a constant fraction.
 

diff --git a/tests/flux_values.json b/tests/flux_values.json
@@ -129,6 +129,28 @@
       2.2713129512457712e-18, 3.871036212755032e-22, 1.8597326553441085e-24
     ]
   },
+  "GlobalFitGST_IT": {
+    "2212": [
+      0.03607536338982165, 8.303604540195788e-5, 1.8863983153209037e-7,
+      3.1276596162531363e-10, 5.209459034660495e-13, 4.1532872002068553e-16,
+      8.179387207124814e-19, 1.6295037238885218e-21, 6.3885346081344595e-27
+    ],
+    "1000020040": [
+      0.023762074349477715, 6.617029782920994e-5, 1.878720144259654e-7,
+      4.685738192425947e-10, 9.255593536028597e-13, 1.4792625219449634e-15,
+      9.643447774808006e-23, 6.700436854786504e-74, 0.0
+    ],
+    "1000080160": [
+      0.003971477965627884, 1.6135814814332014e-5, 6.546771674170436e-8,
+      2.5170413550201714e-10, 5.942565425469719e-13, 7.956149842340593e-16,
+      4.744600882605979e-19, 6.212710347575392e-29, 2.1228946532983858e-106
+    ],
+    "1000260560": [
+      0.0015996433832794234, 8.135139888332144e-6, 4.1420823134587994e-8,
+      2.0697967478899172e-10, 8.354760133167253e-13, 1.0349313641570426e-15,
+      2.2713129512457712e-18, 3.871036212755032e-22, 1.8597326553441085e-24
+    ]
+  },
   "Hoerandel": {
     "2212": [
       0.04477279875079687, 8.72999889028403e-5, 1.7021971482865675e-7,
@@ -534,5 +556,27 @@
       1.5683732598520708e-10, 2.262722571198835e-13, 2.9271639177164597e-16,
       5.2555897942362e-19, 1.8275958248865348e-22, 6.139956685254725e-28
     ]
+  },
+  "GlobalSplineFit_IT": {
+    "2212": [
+      0.04585449431423036, 7.789783588789647e-5, 1.9061451942412818e-7,
+      3.913540332159069e-10, 9.171060628426788e-13, 7.933042844441958e-16,
+      6.150944574738882e-19, 1.1644454409750712e-21, 3.79633922979983e-25
+    ],
+    "1000020040": [
+      0.03192055738163634, 6.467091543445148e-5, 1.7585497031971908e-7,
+      4.1712912655120846e-10, 8.035956231065464e-13, 9.732708568816786e-16,
+      5.416659594465701e-19, 5.335274557933634e-22, 1.3445890398933065e-24
+    ],
+    "1000080160": [
+      0.018496228648837695, 4.259693993066384e-5, 8.770466269694322e-8,
+      2.6164362595183446e-10, 6.022544255316352e-13, 7.62886294617275e-16,
+      9.38795530909692e-19, 3.8207698885288526e-22, 1.1392141934461967e-25
+    ],
+    "1000260560": [
+      0.012020219553048266, 4.9354088622501694e-5, 9.361411920367124e-8,
+      2.9805554501291447e-10, 4.2376209855749804e-13, 5.757781836050234e-16,
+      9.852668748367386e-19, 2.4213191949850687e-22, 1.0847379951951749e-26
+    ]
   }
 }
diff --git a/tests/test_fluxes.py b/tests/test_fluxes.py
@@ -60,12 +60,18 @@ def test_TIG1996(self):
     def test_GlobalFitGST(self):
         self.flux_cmp("GlobalFitGST")
 
+    def test_GlobalFitGST_IT(self):
+        self.flux_cmp("GlobalFitGST_IT")
+
     def test_GlobalSplineFit(self):
         self.flux_cmp("GlobalSplineFit")
 
     def test_GlobalSplineFit5Comp(self):
         self.flux_cmp("GlobalSplineFit5Comp")
 
+    def test_GlobalSplineFit_IT(self):
+        self.flux_cmp("GlobalSplineFit_IT")
+
     def test_FixedFractionFlux(self):
         self.flux_cmp("FixedFractionFlux", {2212: 0.1, 1000020040: 0.2, 1000080160: 0.3, 1000260560: 0.4})
         self.flux_cmp(