Merge pull request #242 from icecube/CutoffPL_fix

Implement CutoffPowerLaw Fluxes with public data

Author: ssclafani949

skyllh/analyses/i3/publicdata_ps/time_integrated_ps.py

@@ -42,6 +42,7 @@
     SpatialBoxEventSelectionMethod,
 )
 from skyllh.core.flux_model import (
+    CutoffPowerLawEnergyFluxProfile,
     PowerLawEnergyFluxProfile,
     SteadyPointlikeFFM,
 )
@@ -134,6 +135,7 @@ def create_analysis(
         refplflux_Phi0=1,
         refplflux_E0=1e3,
         refplflux_gamma=2.0,
+        refplflux_Ec=np.inf,
         ns_seed=10.0,
         ns_min=0.,
         ns_max=1e3,
@@ -171,6 +173,8 @@ def create_analysis(
         The reference energy to use for the reference power law flux model.
     refplflux_gamma : float
         The spectral index to use for the reference power law flux model.
+    refplflux_Ec: float,
+        The cutoff energy for the cutoff power law flux model.
     ns_seed : float
         Value to seed the minimizer with for the ns fit.
     ns_min : float
@@ -269,15 +273,28 @@ def create_analysis(
         dtc_except_fields = ['mcweight', 'time']
 
     # Define the flux model.
-    fluxmodel = SteadyPointlikeFFM(
-        Phi0=refplflux_Phi0,
-        energy_profile=PowerLawEnergyFluxProfile(
-            E0=refplflux_E0,
-            gamma=refplflux_gamma,
+    if refplflux_Ec == np.inf:
+        fluxmodel = SteadyPointlikeFFM(
+            Phi0=refplflux_Phi0,
+            energy_profile=PowerLawEnergyFluxProfile(
+                E0=refplflux_E0,
+                gamma=refplflux_gamma,
+                cfg=cfg,
+            ),
             cfg=cfg,
-        ),
-        cfg=cfg,
-    )
+        )
+    else:
+        fluxmodel = SteadyPointlikeFFM(
+            Phi0=refplflux_Phi0,
+            energy_profile=CutoffPowerLawEnergyFluxProfile(
+                E0=refplflux_E0,
+                gamma=refplflux_gamma,
+                Ecut=refplflux_Ec, #in GeV
+                cfg=cfg,
+            ),
+            cfg=cfg,
+        )
+
 
     # Define the fit parameter ns.
     param_ns = Parameter(

skyllh/core/flux_model.py

@@ -771,6 +771,60 @@ def __call__(
         values *= np.exp(-E / self._Ecut)
 
         return values
+    def get_integral(
+            self,
+            E1,
+            E2,
+            unit=None,
+    ):
+        """This falls back to the default implementation for the CutoffPowerLawEnergyFluxProfile
+        to avoid using the inhereted function
+        .. note::
+
+            This implementation utilizes the ``scipy.integrate.quad`` function
+            to perform a generic numeric integration. Hence, this implementation
+            is slow and should be reimplemented by the derived class if an
+            analytic integral form is available.
+
+        Parameters
+        ----------
+        E1 : float | 1d numpy ndarray of float
+            The lower energy bound of the integration.
+        E2 : float | 1d numpy ndarray of float
+            The upper energy bound of the integration.
+        unit : instance of astropy.units.UnitBase | None
+            The unit of the given energies.
+            If set to ``None``, the set energy unit of this EnergyFluxProfile
+            instance is assumed.
+
+        Returns
+        -------
+        integral : instance of ndarray
+            The (n,)-shaped numpy ndarray holding the integral values of the
+            given integral ranges.
+        """
+        E1 = np.atleast_1d(E1)
+        E2 = np.atleast_1d(E2)
+
+        if (unit is not None) and (unit != self._energy_unit):
+            time_unit_conv_factor = unit.to(self._energy_unit)
+            E1 = E1 * time_unit_conv_factor
+            E2 = E2 * time_unit_conv_factor
+
+        integral = np.empty((len(E1),), dtype=np.float64)
+
+        for (i, (E1_i, E2_i)) in enumerate(zip(E1, E2)):
+
+            # integration for log(energy) for hopefully better numerical stability
+            tmp_e = np.linspace(np.log10(E1_i), np.log10(E2_i), 5000)
+            tmp_int = np.trapz(np.log(10) * self(10**tmp_e) * 10**tmp_e, tmp_e) 
+
+            # make sure it is always positive (probably not an issue any more with np.trapz. 
+            # used to be an issue using the spline integrate self.function.integrate)
+            integral[i] = max(0.0, tmp_int)
+
+        return integral
+
 
 
 class LogParabolaPowerLawEnergyFluxProfile(