Add support for external custom data

docs/source/changes/4.added.rst

@@ -0,0 +1 @@
+Added support for custom external performance data.

docs/source/index.rst

@@ -34,10 +34,15 @@ telescope system given a spectral shape.
 The signal significances of each spectral point are computed according to
 Eq. 17 definition from [2]_.
 
-The currently available performance data shipped with the package is:
-
-- MAGIC *low* zenith angle (0 to 30 degrees) in the range 40GeV-16TeV from [1]_,
-- MAGIC *mid* zenith angle (30 to 45 degrees) in the range 40GeV-16TeV from [1]_.
+The currently available performance data publicly shipped with the package is
+summarized by this table:
+
+==========  ===============  ============  ============  ==========
+Istrument   Zenith range ID  Zenith range  Energy range  References
+==========  ===============  ============  ============  ==========
+MAGIC       low              0 to 30 deg   40GeV-16TeV   [1]_
+MAGIC       mid              30 to 45 deg  40GeV-16TeV   [1]_
+==========  ===============  ============  ============  ==========
 
 Current caveats
 ---------------

docs/source/userguide.rst

@@ -32,6 +32,59 @@ for a source named "my_source",
 
     iact-estimator --config config.yml --source-name my_source
 
+Using custom performance data
+=============================
+
+If you want to use your own performance data,
+either for testing or because it is still not published,
+you can load it from an ECSV with the format shown below.
+
+Producing such a table is very easy with astropy
+starting from the array quantities containing the data,
+for example,
+
+.. code-block:: python
+
+    from astropy.table import QTable
+    import astropy.units as u
+    import numpy as np
+
+    table = QTable([min_energy, max_energy, gamma_rate, bkg_rate],
+                    names=("min_energy", "max_energy", "gamma_rate", "background_rate"),
+                    meta={"name":"some_descriptive_title"})
+    table.write("my_performance.ecsv")
+
+this will result in a data file similar to this,
+
+.. code-block::
+
+    # %ECSV 1.0
+    # ---
+    # datatype:
+    # - {name: min_energy, unit: GeV, datatype: float64}
+    # - {name: max_energy, unit: GeV, datatype: float64}
+    # - {name: gamma_rate, datatype: float64}
+    # - {name: background_rate, datatype: float64}
+    # meta:
+    #   __serialized_columns__:
+    #     max_energy:
+    #       __class__: astropy.units.quantity.Quantity
+    #       unit: &id001 !astropy.units.Unit {unit: GeV}
+    #       value: !astropy.table.SerializedColumn {name: max_energy}
+    #     min_energy:
+    #       __class__: astropy.units.quantity.Quantity
+    #       unit: *id001
+    #       value: !astropy.table.SerializedColumn {name: min_energy}
+    #   name: some_descriptive_title
+    #   offset_degradation: false
+    # schema: astropy-2.0
+    min_energy max_energy gamma_rate background_rate
+    20 40 7.2 359.0
+    ...etc
+
+You can then load it using the `--performance` flag of :ref:`iact-estimator`
+to tell the command-line tool where to find the data file.
+
 Output
 ======
 

src/iact_estimator/core.py

@@ -4,7 +4,6 @@
 import importlib
 
 import astropy.units as u
-from gammapy.stats import WStatCountsStatistic
 import numpy as np
 from scipy import interpolate
 from scipy.integrate import quad
@@ -15,14 +14,14 @@
 )
 from .io import load_ebl
 from .spectral import crab_nebula_spectrum
+from .statistics import probability_to_sigma, sigma_to_probability, significance_li_ma
 
 __all__ = [
     "setup_logging",
     "check_input_configuration",
     "initialize_model",
     "observed_flux",
     "get_sed",
-    "significance_li_ma",
     "prepare_data",
     "source_detection",
     "calculate",
@@ -133,7 +132,7 @@ def get_horizon_stereo_profile(M1_data, M2_data):
     return az * u.deg, zd * u.deg
 
 
-def check_input_configuration(config):
+def check_input_configuration(config, performance_data):
     """
     Import and initialize a spectral model.
 
@@ -151,6 +150,9 @@ def check_input_configuration(config):
     # Assume a valid configuration
     is_valid = True
 
+    if performance_data:
+        performance_metadata = performance_data.meta
+
     extension = u.Quantity(config["extension"]).to_value("deg")
 
     if extension > 1:
@@ -167,9 +169,13 @@ def check_input_configuration(config):
             " region (n_off_regions) should be used."
         )
         is_valid = False
-    if config["sum_trigger"] and (config["zenith_performance"] == "mid"):
-        logger.warning("This functionality has not been yet implemented.")
-        # raise NotImplementedError("This functionality has not been yet implemented.")
+    if config["sum_trigger"] and (config["zenith_range"] == "mid"):
+        is_valid = False
+        raise NotImplementedError(
+            "MAGIC SUM trigger at mid zenith range has not been yet implemented."
+        )
+    if ("LST" in performance_metadata) and config["sum_trigger"]:
+        logger.warning("LST mode is not compatible with SUMT")
         is_valid = False
     if config["offset_degradation_factor"] > 1.00001:
         logger.warning(
@@ -231,7 +237,7 @@ def initialize_model(config):
 
 
 @u.quantity_input(energy=u.TeV)
-def observed_flux(energy, redshift, flux_int):
+def observed_flux(energy, redshift, flux_int, ebl_file_path):
     """
     Get the attenuated flux of the source.
 
@@ -249,7 +255,7 @@ def observed_flux(energy, redshift, flux_int):
     valid : `bool`
         Initialized instance of a spectral model.
     """
-    ebl_redshift, ebl_energy, ebl_taus = load_ebl()
+    ebl_redshift, ebl_energy, ebl_taus = load_ebl(ebl_file_path)
     ftau = interpolate.interp2d(
         ebl_redshift, ebl_energy.to_value(energy.unit), ebl_taus, kind="cubic"
     )
@@ -280,43 +286,7 @@ def get_sed(energy, flux):
     return sed
 
 
-def significance_li_ma(n_on, n_off, alpha, mu_sig=None):
-    """
-    Get the Li & Ma significance.
-
-    This is equivalent to eq.17 of [1]_.
-
-    Parameters
-    ----------
-    n_on : `int`
-        Measured counts in ON region.
-    n_off : `int`
-        Measured counts in OFF region.
-    alpha : `float`
-        Acceptance ratio of ON and OFF measurements.
-    mu_sig : `float`
-        Expected signal counts in ON region.
-
-    Returns
-    -------
-    sqrt_ts : `float``
-        Significance as the square root of the Test Statistic.
-
-    Notes
-    -----
-    The implementation uses `gammapy.stats.WStatCountsStatistic`
-    and takes the square root of the Test Statistic.
-
-    References
-    ----------
-    .. [1] Li, T.-P. & Ma, Y.-Q., ApJ, 1983, 272, 317, 10.1086/161295.
-    """
-    statistics = WStatCountsStatistic(n_on, n_off, alpha, mu_sig)
-    sqrt_ts = statistics.sqrt_ts
-    return sqrt_ts
-
-
-def prepare_data(config):
+def prepare_data(config, performance_data=None):
     """
     Extract the performance data.
 
@@ -335,22 +305,28 @@ def prepare_data(config):
         Rate of background events from performance data.
     """
 
-    performance_data = {
-        "low": LOW_ZENITH_PERFORMANCE,
-        "mid": MID_ZENITH_PERFORMANCE,
-    }
+    if not performance_data:
+        if config["sum_trigger"] and config["zenith_range"] == "mid":
+            message = "MAGIC Mid zenith performance with the SUM trigger is not currently available."
+            logger.critical(message)
+            raise NotImplementedError(message)
 
-    if config["sum_trigger"] and config["zenith_performance"] == "low":
-        message = "Low zenith performance with the SUM trigger"
-        " is not currently available."
-        logger.critical(message)
-        raise NotImplementedError(message)
+        # use packaged (public) datasets
+        available_datasets = {
+            "low": LOW_ZENITH_PERFORMANCE,
+            "mid": MID_ZENITH_PERFORMANCE,
+        }
 
-    min_energy = performance_data[config["zenith_performance"]]["min_energy"]
-    max_energy = performance_data[config["zenith_performance"]]["max_energy"]
+        performance_data = available_datasets[config["zenith_range"]]
+
+    min_energy = performance_data["min_energy"]
+    max_energy = performance_data["max_energy"]
     energy_bins = np.append(min_energy.value, max_energy[-1].value) * min_energy.unit
-    gamma_rate = performance_data[config["zenith_performance"]]["gamma_rate"]
-    background_rate = performance_data[config["zenith_performance"]]["background_rate"]
+    gamma_rate = performance_data["gamma_rate"]
+    background_rate = performance_data["background_rate"]
+
+    gamma_rate *= config["offset_degradation_factor"]
+    background_rate *= config["offset_degradation_factor"]
 
     return energy_bins, gamma_rate, background_rate
 
@@ -375,12 +351,21 @@ def source_detection(sigmas, observation_time):
 
     time = observation_time.to("h")
 
+    combined_probability = 1
+    combined_significance_text = ""
     if len(sigmas) > 0:
-        combined_significance = sum(sigmas) / np.sqrt(len(sigmas))
+        combined_probability = np.prod(sigma_to_probability(sigmas))
+        combined_significance = probability_to_sigma(combined_probability)
+        combined_significance_text = "{0:.2f}".format(combined_significance)
+        if (
+            combined_probability < 1.0e-307
+        ):  # numerical accuracy problem, but significance will be either way large
+            combined_significance = 38  # or more ...
+            combined_significance_text = ">38"
 
         print(
             f"Combined significance (using the {len(sigmas):d} data points"
-            f" shown in the SED) = {combined_significance:.1f}"
+            f" shown in the SED) = {combined_significance_text}"
         )
     else:
         print(f"The source will not be detected in {time}.")

src/iact_estimator/io.py

@@ -3,6 +3,7 @@
 import logging
 from pathlib import Path
 
+from astropy.table import QTable
 from numpy import loadtxt
 from yaml import load
 
@@ -11,7 +12,7 @@
 except ImportError:
     from yaml import Loader
 
-__all__ = ["read_yaml", "load_ebl"]
+__all__ = ["read_yaml", "load_ebl", "load_performance_ecsv"]
 
 logger = logging.getLogger(__name__)
 
@@ -69,3 +70,28 @@ def load_ebl(ebl_file_path):
     taus = ebl_body[:, 1:]
     if len(taus > 0):
         return zz, energies, taus
+
+
+def load_performance_ecsv(input_file_path):
+    """
+    Load performance data from an ECSV file as a dictionary.
+
+    Parameters
+    ----------
+    input_file_path : `str`
+        Path to the input ECSV file.
+
+    Returns
+    -------
+    table : `~astropy.table.QTable`
+        Contents of the YAML file in form
+        of a Python dictionary.
+    """
+    input_file_path = Path(input_file_path).resolve()
+    try:
+        table = QTable.read(input_file_path)
+    except FileNotFoundError:
+        logger.exception(
+            "Performance file not found at %s", input_file_path, exc_info=True
+        )
+    return table

src/iact_estimator/plots.py

@@ -149,7 +149,7 @@ def plot_sed(
     ax.set(
         xlabel=rf"Energy [ {energy_unit} ]",
         ylabel=rf"$E^{2}$ $\dfrac{{dN}}{{dE dA dt}}$ [ {energy_flux_unit.to_string('latex', fraction=False)} ]",
-        title=rf"$\Sigma\sigma_i / \sqrt{{{len(sigmas)}}} = {{{round(sig_comb.value, 1)}}}$",
+        title=rf"$\Sigma\sigma_i / \sqrt{{{len(sigmas)}}} = {{{round(sig_comb, 1)}}}$",
     )
     energy = np.logspace(
         np.log10(min_energy.to_value(energy_unit)),

src/iact_estimator/resources/config.yml

@@ -13,7 +13,13 @@ observation_datetime: "2024-06-15 18:00"
 extension: 0.0 deg
 redshift: -1 # redshift of the source (for the EBL absorption), if -1 then no absorption
 sum_trigger: False
-zenith_performance: "low" # "low"=0-30deg, "mid"=30-45deg
+
+# see documentation for a table of publicly available
+# values
+# For private data, these values should come with the
+# metadata of your file
+instrument: "MAGIC"
+zenith_range: "low"
 # you can check visibility of your source e.g. here: http://www.magic.iac.es/scheduler/
 
 ebl_file_path:

src/iact_estimator/scripts/estimator.py

@@ -10,7 +10,7 @@
 from astropy.visualization import quantity_support
 import matplotlib.pyplot as plt
 
-from ..io import read_yaml
+from ..io import read_yaml, load_performance_ecsv
 from ..core import (
     setup_logging,
     initialize_model,
@@ -32,6 +32,12 @@
     type=str,
     help="Name of the source to estimate.",
 )
+parser.add_argument(
+    "--performance",
+    default="",
+    type=str,
+    help="Custom performance data.",
+)
 parser.add_argument(
     "--output-path",
     default=None,
@@ -67,8 +73,14 @@ def main():
     logger.info("Loading configuration file")
     config = read_yaml(args.config)
 
+    performance_data = None
+    if args.performance:
+        performance_data_file = Path(args.performance).resolve()
+        logger.info("Loading performance data from %s", performance_data_file)
+        performance_data = load_performance_ecsv(performance_data_file)
+
     logging.info("Validating input configuration")
-    if not check_input_configuration(config):
+    if not check_input_configuration(config, performance_data):
         logging.critical("One or more invalid configuration settings have been found.")
         parser.exit(status=1, message="iact-estimator terminated abnormally.")
 
@@ -101,7 +113,10 @@ def main():
             label=source_name,
         )
 
-    energy_bins, gamma_rate, background_rate = prepare_data(config)
+    if not performance_data:
+        energy_bins, gamma_rate, background_rate = prepare_data(config)
+    else:
+        energy_bins, gamma_rate, background_rate = prepare_data(performance_data)
 
     en, sed, dsed, sigmas, detected = calculate(
         energy_bins, gamma_rate, background_rate, config, assumed_spectrum