From 38e0baa42437a44dfd6d2034701108672c2aed0e Mon Sep 17 00:00:00 2001 From: Gabriel Emery Date: Wed, 12 Feb 2025 10:46:48 +0100 Subject: [PATCH 1/2] Implement custom energy-dependent cuts for IRFs and DL3 creations. --- magicctapipe/io/__init__.py | 2 + magicctapipe/io/io.py | 49 ++++++ magicctapipe/io/tests/test_io.py | 40 +++++ magicctapipe/resources/config.yaml | 20 ++- magicctapipe/resources/test_config.yaml | 21 ++- .../lst1_magic/lst1_magic_create_irf.py | 154 ++++++++++++------ 6 files changed, 227 insertions(+), 59 deletions(-) diff --git a/magicctapipe/io/__init__.py b/magicctapipe/io/__init__.py index e1cf9fcc..bfdc7aee 100644 --- a/magicctapipe/io/__init__.py +++ b/magicctapipe/io/__init__.py @@ -15,6 +15,7 @@ from .io import ( check_input_list, format_object, + get_custom_cuts, get_dl2_mean, get_stereo_events, get_stereo_events_old, @@ -40,6 +41,7 @@ "create_gti_hdu", "create_pointing_hdu", "format_object", + "get_custom_cuts", "get_dl2_mean", "get_stereo_events", "get_stereo_events_old", diff --git a/magicctapipe/io/io.py b/magicctapipe/io/io.py index 6009a12a..fceaffa6 100644 --- a/magicctapipe/io/io.py +++ b/magicctapipe/io/io.py @@ -22,6 +22,7 @@ from pyirf.binning import join_bin_lo_hi from pyirf.simulations import SimulatedEventsInfo from pyirf.utils import calculate_source_fov_offset, calculate_theta +from scipy.interpolate import interp1d try: from importlib.resources import files @@ -33,6 +34,7 @@ __all__ = [ "check_input_list", "format_object", + "get_custom_cuts", "get_dl2_mean", "get_stereo_events", "get_stereo_events_old", @@ -249,6 +251,53 @@ def format_object(input_object): return string +def get_custom_cuts(config, energy_bins): + """ + Extract and interpolate custom, energy-dependent cuts for the IRFs creation. + + Parameters + ---------- + config : dict + Dictionary containing the custom cuts and the interpolation tragedy. + Required entries are a sub-dictionary 'custom_cuts' and a field + 'interpolate_kind' used in scipy.interpolate.interp1d. + energy_bins : astropy.units.Quantity + Energy bin edges targeted in the IRFs + + Returns + ------- + astropy.table.QTable + Custom cut table in target 'energy_bins' + """ + cut_table = QTable() + interpolate_kind = config["interpolate_kind"] + energy_bins_cuts = config["custom_cuts"]["energy_bins"] * u.Unit( + config["custom_cuts"]["energy_unit"] + ) + energy_bins_cuts_center = 0.5 * (energy_bins_cuts[:-1] + energy_bins_cuts[1:]).to( + energy_bins.unit + ) + cuts = config["custom_cuts"]["cut"] * u.Unit(config["custom_cuts"]["cut_unit"]) + + cut_table["low"] = energy_bins[:-1] + cut_table["high"] = energy_bins[1:] + cut_table["center"] = 0.5 * (cut_table["low"] + cut_table["high"]) + f_cut = interp1d( + energy_bins_cuts_center, + cuts, + kind=interpolate_kind, + bounds_error=False, + fill_value=(cuts[0], cuts[-1]), + assume_sorted=True, + ) + + cut_table["cut"] = f_cut(cut_table["center"]) * u.Unit( + config["custom_cuts"]["cut_unit"] + ) + + return cut_table + + def get_stereo_events_old( event_data, quality_cuts=None, group_index=["obs_id", "event_id"] ): diff --git a/magicctapipe/io/tests/test_io.py b/magicctapipe/io/tests/test_io.py index 80505e85..5c2c2659 100644 --- a/magicctapipe/io/tests/test_io.py +++ b/magicctapipe/io/tests/test_io.py @@ -1,13 +1,17 @@ import glob +import astropy.table +import astropy.units as u import numpy as np import pandas as pd import pytest +import yaml from ctapipe_io_lst import REFERENCE_LOCATION from magicctapipe.io.io import ( check_input_list, format_object, + get_custom_cuts, get_dl2_mean, get_stereo_events, load_dl2_data_file, @@ -152,6 +156,42 @@ def test_check_input_list(self): } ) + def test_get_custom_cuts(self, config): + with open(config, "rb") as f: + config = yaml.safe_load(f) + # Create energy bins + config_eng_bins = config["create_irf"]["energy_bins"] + eng_bins_start = u.Quantity(config_eng_bins["start"]) + eng_bins_stop = u.Quantity(config_eng_bins["stop"]) + + energy_bins = u.TeV * np.geomspace( + start=eng_bins_start.to_value("TeV").round(3), + stop=eng_bins_stop.to_value("TeV").round(3), + num=config_eng_bins["n_edges"], + ) + + # Load the cuts configs. 'cut_type' will be ignored + config_gh_irf = config["create_irf"]["gammaness"] + config_theta_irf = config["create_irf"]["theta"] + + # Get the custom cuts without interpolating on energy_bins + cut_table_gh = get_custom_cuts(config_gh_irf, energy_bins) + cut_table_theta = get_custom_cuts(config_theta_irf, energy_bins) + assert isinstance(cut_table_gh, astropy.table.QTable) + assert isinstance(cut_table_theta, astropy.table.QTable) + assert cut_table_gh["low"].unit == "TeV" + assert cut_table_gh["high"].unit == "TeV" + assert cut_table_gh["center"].unit == "TeV" + assert cut_table_gh["cut"].unit == u.dimensionless_unscaled + assert cut_table_gh["cut"][15] == 0.5 + assert cut_table_theta["cut"].unit == u.deg + assert cut_table_theta["cut"][15] == 0.3 * u.deg + + # Get the custom cuts with linear interpolating on energy_bins + config_gh_irf["interpolate_kind"] = "linear" + cut_table_gh = get_custom_cuts(config_gh_irf, energy_bins) + assert np.isclose(cut_table_gh["cut"][15].value, 0.5074, atol=0.0001) + def test_telescope_combinations(self, config_gen, config_gen_4lst): """ Simple check on telescope combinations diff --git a/magicctapipe/resources/config.yaml b/magicctapipe/resources/config.yaml index 8b4e14c6..a6870710 100644 --- a/magicctapipe/resources/config.yaml +++ b/magicctapipe/resources/config.yaml @@ -256,18 +256,34 @@ create_irf: n_edges: 21 gammaness: - cut_type: "dynamic" # select "global" or "dynamic" + cut_type: "dynamic" # select "global" or "dynamic" or "custom" global_cut_value: 0.8 # used for the global cut efficiency: 0.9 # used for the dynamic cuts min_cut: 0.05 # used for the dynamic cuts max_cut: 0.85 # used for the dynamic cuts + custom_cuts: { # used for custom cuts + energy_bins: [0.01, 0.1, 1, 10, 100], + energy_unit: "TeV", + cut: [0.4, 0.45, 0.50, 0.55], + cut_unit: "", + } + interpolate_kind: "nearest" # used for custom cuts, option for scipy.interpolate.interp1d + # Standard values are "nearest", "linear" and "cubic" theta: - cut_type: "dynamic" # select "global" or "dynamic" + cut_type: "dynamic" # select "global" or "dynamic" or "custom" global_cut_value: "0.2 deg" # used for the global cut efficiency: 0.75 # used for the dynamic cuts min_cut: "0.1 deg" # used for the dynamic cuts max_cut: "0.25 deg" # used for the dynamic cuts + custom_cuts: { # used for custom cuts + energy_bins: [0.01, 0.1, 1, 10, 100], + energy_unit: "TeV", + cut: [0.4, 0.35, 0.30, 0.25], + cut_unit: "deg", + } + interpolate_kind: "nearest" # used for custom cuts, option for scipy.interpolate.interp1d + # Standard values are "nearest", "linear" and "cubic" dl2_to_dl3: diff --git a/magicctapipe/resources/test_config.yaml b/magicctapipe/resources/test_config.yaml index 04b04a73..3b80c2ee 100644 --- a/magicctapipe/resources/test_config.yaml +++ b/magicctapipe/resources/test_config.yaml @@ -226,19 +226,34 @@ create_irf: n_edges: 21 gammaness: - cut_type: "dynamic" # select "global" or "dynamic" + cut_type: "dynamic" # select "global" or "dynamic" or "custom" global_cut_value: 0.8 # used for the global cut efficiency: 0.9 # used for the dynamic cuts min_cut: 0.05 # used for the dynamic cuts max_cut: 0.85 # used for the dynamic cuts + custom_cuts: { # used for custom cuts + energy_bins: [0.01, 0.1, 1, 10, 100], + energy_unit: "TeV", + cut: [0.4, 0.45, 0.50, 0.55], + cut_unit: "", + } + interpolate_kind: "nearest" # used for custom cuts, option for scipy.interpolate.interp1d + # Standard values are "nearest", "linear" and "cubic" theta: - cut_type: "dynamic" # select "global" or "dynamic" + cut_type: "dynamic" # select "global" or "dynamic" or "custom" global_cut_value: "0.2 deg" # used for the global cut efficiency: 0.75 # used for the dynamic cuts min_cut: "0.1 deg" # used for the dynamic cuts max_cut: "0.3 deg" # used for the dynamic cuts - + custom_cuts: { # used for custom cuts + energy_bins: [0.01, 0.1, 1, 10, 100], + energy_unit: "TeV", + cut: [0.4, 0.35, 0.30, 0.25], + cut_unit: "deg", + } + interpolate_kind: "nearest" # used for custom cuts, option for scipy.interpolate.interp1d + # Standard values are "nearest", "linear" and "cubic" dl2_to_dl3: interpolation_method: "nearest" # select "nearest", "linear" or "cubic" diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py index 25e81920..5749e9e9 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py @@ -23,14 +23,17 @@ The "software(_only_3tel)" types are supposed for the software event coincidence with LST-mono and MAGIC-stereo observations, allowing for only the events triggering both M1 and M2. The "software" type allows -for the events of the any 2-tel combinations except the MAGIC-stereo +for the events of any 2-tel combinations except the MAGIC-stereo combination at the moment. The "software_only_3tel" type allows for only the events of the 3-tel combination. The "magic_only" type allows for only the events of the MAGIC-stereo combination. -There are two types of gammaness and theta cuts, "global" and "dynamic". +There are three types of gammaness and theta cuts, "global", "dynamic" and +"custom". In case of the dynamic cuts, the optimal cut satisfying a given efficiency will be calculated for every energy bin. +In case of custom cuts, the provided cuts with be applied to the +target energy binning using the selected scipy.interpolate.interpd1d method. Usage: $ python lst1_magic_create_irf.py @@ -74,7 +77,12 @@ calculate_event_weights, ) -from magicctapipe.io import create_gh_cuts_hdu, format_object, load_mc_dl2_data_file +from magicctapipe.io import ( + create_gh_cuts_hdu, + format_object, + get_custom_cuts, + load_mc_dl2_data_file, +) __all__ = ["create_irf"] @@ -110,6 +118,9 @@ def create_irf( If the input type of gammaness or theta cut is not known """ + dynamic_keys = ["efficiency", "min_cut", "max_cut"] + custom_keys = ["custom_cuts", "interpolate_kind"] + config_irf = config["create_irf"] quality_cuts = config_irf["quality_cuts"] @@ -335,34 +346,47 @@ def create_irf( mask_gh = event_table_bkg["gammaness"] > cut_value_gh event_table_bkg = event_table_bkg[mask_gh] - elif cut_type_gh == "dynamic": - config_gh_cuts.pop("global_cut_value", None) - - logger.info("\nDynamic gammaness cuts:") - logger.info(format_object(config_gh_cuts)) + elif cut_type_gh == "dynamic" or "custom": + if cut_type_gh == "dynamic": + logger.info("\nDynamic gammaness cuts:") + logger.info( + format_object( + {k: config_gh_cuts[k] for k in dynamic_keys if k in config_gh_cuts} + ) + ) - gh_efficiency = config_gh_cuts["efficiency"] - gh_cut_min = config_gh_cuts["min_cut"] - gh_cut_max = config_gh_cuts["max_cut"] + gh_efficiency = config_gh_cuts["efficiency"] + gh_cut_min = config_gh_cuts["min_cut"] + gh_cut_max = config_gh_cuts["max_cut"] - extra_header["GH_EFF"] = gh_efficiency - extra_header["GH_MIN"] = gh_cut_min - extra_header["GH_MAX"] = gh_cut_max + extra_header["GH_EFF"] = gh_efficiency + extra_header["GH_MIN"] = gh_cut_min + extra_header["GH_MAX"] = gh_cut_max - output_suffix = f"gh_dyn{gh_efficiency}" + output_suffix = f"gh_dyn{gh_efficiency}" - # Calculate dynamic gammaness cuts - gh_percentile = 100 * (1 - gh_efficiency) + # Calculate dynamic gammaness cuts + gh_percentile = 100 * (1 - gh_efficiency) - cut_table_gh = calculate_percentile_cut( - values=event_table_gamma["gammaness"], - bin_values=event_table_gamma["reco_energy"], - bins=energy_bins, - fill_value=gh_cut_min, - percentile=gh_percentile, - min_value=gh_cut_min, - max_value=gh_cut_max, - ) + cut_table_gh = calculate_percentile_cut( + values=event_table_gamma["gammaness"], + bin_values=event_table_gamma["reco_energy"], + bins=energy_bins, + fill_value=gh_cut_min, + percentile=gh_percentile, + min_value=gh_cut_min, + max_value=gh_cut_max, + ) + if cut_type_gh == "custom": + cut_table_gh = get_custom_cuts(config_gh_cuts, energy_bins) + extra_header["GH_CUTS"] = "Custom energy-dependent" + output_suffix = "gh_custom" + logger.info("\nCustom gammaness cuts.") + logger.info( + format_object( + {k: config_gh_cuts[k] for k in custom_keys if k in config_gh_cuts} + ) + ) logger.info(f"\nGammaness-cut table:\n\n{cut_table_gh}") @@ -422,34 +446,56 @@ def create_irf( mask_theta = event_table_gamma["theta"].to_value("deg") < cut_value_theta event_table_gamma = event_table_gamma[mask_theta] - elif cut_type_theta == "dynamic": - config_theta_cuts.pop("global_cut_value", None) - - logger.info("\nDynamic theta cuts:") - logger.info(format_object(config_theta_cuts)) - - theta_efficiency = config_theta_cuts["efficiency"] - theta_cut_min = u.Quantity(config_theta_cuts["min_cut"]) - theta_cut_max = u.Quantity(config_theta_cuts["max_cut"]) - - extra_header["TH_EFF"] = theta_efficiency - extra_header["TH_MIN"] = (theta_cut_min.to_value("deg"), "deg") - extra_header["TH_MAX"] = (theta_cut_max.to_value("deg"), "deg") - - output_suffix += f"_theta_dyn{theta_efficiency}" - - # Calculate dynamic theta cuts - theta_percentile = 100 * theta_efficiency - - cut_table_theta = calculate_percentile_cut( - values=event_table_gamma["theta"], - bin_values=event_table_gamma["reco_energy"], - bins=energy_bins, - fill_value=theta_cut_max, - percentile=theta_percentile, - min_value=theta_cut_min, - max_value=theta_cut_max, - ) + elif cut_type_theta == "dynamic" or "custom": + if cut_type_theta == "dynamic": + + logger.info("\nDynamic theta cuts:") + logger.info( + format_object( + { + k: config_theta_cuts[k] + for k in dynamic_keys + if k in config_theta_cuts + } + ) + ) + + theta_efficiency = config_theta_cuts["efficiency"] + theta_cut_min = u.Quantity(config_theta_cuts["min_cut"]) + theta_cut_max = u.Quantity(config_theta_cuts["max_cut"]) + + extra_header["TH_EFF"] = theta_efficiency + extra_header["TH_MIN"] = (theta_cut_min.to_value("deg"), "deg") + extra_header["TH_MAX"] = (theta_cut_max.to_value("deg"), "deg") + + output_suffix += f"_theta_dyn{theta_efficiency}" + + # Calculate dynamic theta cuts + theta_percentile = 100 * theta_efficiency + + cut_table_theta = calculate_percentile_cut( + values=event_table_gamma["theta"], + bin_values=event_table_gamma["reco_energy"], + bins=energy_bins, + fill_value=theta_cut_max, + percentile=theta_percentile, + min_value=theta_cut_min, + max_value=theta_cut_max, + ) + if cut_type_theta == "custom": + cut_table_theta = get_custom_cuts(config_theta_cuts, energy_bins) + extra_header["THETA_CUTS"] = "Custom energy-dependent" + output_suffix += "_theta_custom" + logger.info("\nCustom theta cuts.") + logger.info( + format_object( + { + k: config_theta_cuts[k] + for k in custom_keys + if k in config_theta_cuts + } + ) + ) logger.info(f"\nTheta-cut table:\n\n{cut_table_theta}") From 9814e2b29a1c8de1a7b2df86e6ff06a9421711be Mon Sep 17 00:00:00 2001 From: Gabriel Emery Date: Fri, 30 May 2025 16:33:12 +0200 Subject: [PATCH 2/2] Reduce nested if level. Fix typo. Reduce duplicated code in tests. --- magicctapipe/io/io.py | 2 +- magicctapipe/io/tests/test_io.py | 11 +- .../lst1_magic/lst1_magic_create_irf.py | 165 +++++++++--------- 3 files changed, 88 insertions(+), 90 deletions(-) diff --git a/magicctapipe/io/io.py b/magicctapipe/io/io.py index fceaffa6..29fbd170 100644 --- a/magicctapipe/io/io.py +++ b/magicctapipe/io/io.py @@ -258,7 +258,7 @@ def get_custom_cuts(config, energy_bins): Parameters ---------- config : dict - Dictionary containing the custom cuts and the interpolation tragedy. + Dictionary containing the custom cuts and the interpolation strategy. Required entries are a sub-dictionary 'custom_cuts' and a field 'interpolate_kind' used in scipy.interpolate.interp1d. energy_bins : astropy.units.Quantity diff --git a/magicctapipe/io/tests/test_io.py b/magicctapipe/io/tests/test_io.py index 5c2c2659..43b5d30c 100644 --- a/magicctapipe/io/tests/test_io.py +++ b/magicctapipe/io/tests/test_io.py @@ -177,14 +177,13 @@ def test_get_custom_cuts(self, config): # Get the custom cuts without interpolating on energy_bins cut_table_gh = get_custom_cuts(config_gh_irf, energy_bins) cut_table_theta = get_custom_cuts(config_theta_irf, energy_bins) - assert isinstance(cut_table_gh, astropy.table.QTable) - assert isinstance(cut_table_theta, astropy.table.QTable) - assert cut_table_gh["low"].unit == "TeV" - assert cut_table_gh["high"].unit == "TeV" - assert cut_table_gh["center"].unit == "TeV" + for cut_table in (cut_table_gh, cut_table_theta): + assert isinstance(cut_table, astropy.table.QTable) + assert cut_table["low"].unit == "TeV" + assert cut_table["high"].unit == "TeV" + assert cut_table["center"].unit == "TeV" assert cut_table_gh["cut"].unit == u.dimensionless_unscaled assert cut_table_gh["cut"][15] == 0.5 - assert cut_table_theta["cut"].unit == u.deg assert cut_table_theta["cut"][15] == 0.3 * u.deg # Get the custom cuts with linear interpolating on energy_bins diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py index 5749e9e9..c76153ba 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py +++ b/magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py @@ -346,48 +346,51 @@ def create_irf( mask_gh = event_table_bkg["gammaness"] > cut_value_gh event_table_bkg = event_table_bkg[mask_gh] - elif cut_type_gh == "dynamic" or "custom": - if cut_type_gh == "dynamic": - logger.info("\nDynamic gammaness cuts:") - logger.info( - format_object( - {k: config_gh_cuts[k] for k in dynamic_keys if k in config_gh_cuts} - ) + elif cut_type_gh == "dynamic": + logger.info("\nDynamic gammaness cuts:") + logger.info( + format_object( + {k: config_gh_cuts[k] for k in dynamic_keys if k in config_gh_cuts} ) + ) - gh_efficiency = config_gh_cuts["efficiency"] - gh_cut_min = config_gh_cuts["min_cut"] - gh_cut_max = config_gh_cuts["max_cut"] + gh_efficiency = config_gh_cuts["efficiency"] + gh_cut_min = config_gh_cuts["min_cut"] + gh_cut_max = config_gh_cuts["max_cut"] - extra_header["GH_EFF"] = gh_efficiency - extra_header["GH_MIN"] = gh_cut_min - extra_header["GH_MAX"] = gh_cut_max + extra_header["GH_EFF"] = gh_efficiency + extra_header["GH_MIN"] = gh_cut_min + extra_header["GH_MAX"] = gh_cut_max - output_suffix = f"gh_dyn{gh_efficiency}" + output_suffix = f"gh_dyn{gh_efficiency}" - # Calculate dynamic gammaness cuts - gh_percentile = 100 * (1 - gh_efficiency) + # Calculate dynamic gammaness cuts + gh_percentile = 100 * (1 - gh_efficiency) - cut_table_gh = calculate_percentile_cut( - values=event_table_gamma["gammaness"], - bin_values=event_table_gamma["reco_energy"], - bins=energy_bins, - fill_value=gh_cut_min, - percentile=gh_percentile, - min_value=gh_cut_min, - max_value=gh_cut_max, - ) - if cut_type_gh == "custom": - cut_table_gh = get_custom_cuts(config_gh_cuts, energy_bins) - extra_header["GH_CUTS"] = "Custom energy-dependent" - output_suffix = "gh_custom" - logger.info("\nCustom gammaness cuts.") - logger.info( - format_object( - {k: config_gh_cuts[k] for k in custom_keys if k in config_gh_cuts} - ) + cut_table_gh = calculate_percentile_cut( + values=event_table_gamma["gammaness"], + bin_values=event_table_gamma["reco_energy"], + bins=energy_bins, + fill_value=gh_cut_min, + percentile=gh_percentile, + min_value=gh_cut_min, + max_value=gh_cut_max, + ) + + elif cut_type_gh == "custom": + cut_table_gh = get_custom_cuts(config_gh_cuts, energy_bins) + extra_header["GH_CUTS"] = "Custom energy-dependent" + output_suffix = "gh_custom" + logger.info("\nCustom gammaness cuts.") + logger.info( + format_object( + {k: config_gh_cuts[k] for k in custom_keys if k in config_gh_cuts} ) + ) + else: + raise ValueError(f"Unknown gammaness-cut type '{cut_type_gh}'.") + if cut_type_gh == "dynamic" or "custom": logger.info(f"\nGammaness-cut table:\n\n{cut_table_gh}") # Apply the dynamic gammaness cuts @@ -425,9 +428,6 @@ def create_irf( irf_hdus.append(hdu_gh_cuts) - else: - raise ValueError(f"Unknown gammaness-cut type '{cut_type_gh}'.") - if is_point_like: # Apply the theta cut config_theta_cuts = config_irf["theta"] @@ -446,57 +446,59 @@ def create_irf( mask_theta = event_table_gamma["theta"].to_value("deg") < cut_value_theta event_table_gamma = event_table_gamma[mask_theta] - elif cut_type_theta == "dynamic" or "custom": - if cut_type_theta == "dynamic": - - logger.info("\nDynamic theta cuts:") - logger.info( - format_object( - { - k: config_theta_cuts[k] - for k in dynamic_keys - if k in config_theta_cuts - } - ) + elif cut_type_theta == "dynamic": + + logger.info("\nDynamic theta cuts:") + logger.info( + format_object( + { + k: config_theta_cuts[k] + for k in dynamic_keys + if k in config_theta_cuts + } ) + ) - theta_efficiency = config_theta_cuts["efficiency"] - theta_cut_min = u.Quantity(config_theta_cuts["min_cut"]) - theta_cut_max = u.Quantity(config_theta_cuts["max_cut"]) + theta_efficiency = config_theta_cuts["efficiency"] + theta_cut_min = u.Quantity(config_theta_cuts["min_cut"]) + theta_cut_max = u.Quantity(config_theta_cuts["max_cut"]) - extra_header["TH_EFF"] = theta_efficiency - extra_header["TH_MIN"] = (theta_cut_min.to_value("deg"), "deg") - extra_header["TH_MAX"] = (theta_cut_max.to_value("deg"), "deg") + extra_header["TH_EFF"] = theta_efficiency + extra_header["TH_MIN"] = (theta_cut_min.to_value("deg"), "deg") + extra_header["TH_MAX"] = (theta_cut_max.to_value("deg"), "deg") - output_suffix += f"_theta_dyn{theta_efficiency}" + output_suffix += f"_theta_dyn{theta_efficiency}" - # Calculate dynamic theta cuts - theta_percentile = 100 * theta_efficiency + # Calculate dynamic theta cuts + theta_percentile = 100 * theta_efficiency - cut_table_theta = calculate_percentile_cut( - values=event_table_gamma["theta"], - bin_values=event_table_gamma["reco_energy"], - bins=energy_bins, - fill_value=theta_cut_max, - percentile=theta_percentile, - min_value=theta_cut_min, - max_value=theta_cut_max, - ) - if cut_type_theta == "custom": - cut_table_theta = get_custom_cuts(config_theta_cuts, energy_bins) - extra_header["THETA_CUTS"] = "Custom energy-dependent" - output_suffix += "_theta_custom" - logger.info("\nCustom theta cuts.") - logger.info( - format_object( - { - k: config_theta_cuts[k] - for k in custom_keys - if k in config_theta_cuts - } - ) + cut_table_theta = calculate_percentile_cut( + values=event_table_gamma["theta"], + bin_values=event_table_gamma["reco_energy"], + bins=energy_bins, + fill_value=theta_cut_max, + percentile=theta_percentile, + min_value=theta_cut_min, + max_value=theta_cut_max, + ) + elif cut_type_theta == "custom": + cut_table_theta = get_custom_cuts(config_theta_cuts, energy_bins) + extra_header["THETA_CUTS"] = "Custom energy-dependent" + output_suffix += "_theta_custom" + logger.info("\nCustom theta cuts.") + logger.info( + format_object( + { + k: config_theta_cuts[k] + for k in custom_keys + if k in config_theta_cuts + } ) + ) + else: + raise ValueError(f"Unknown theta-cut type '{cut_type_theta}'.") + if cut_type_theta == "dynamic" or "custom": logger.info(f"\nTheta-cut table:\n\n{cut_table_theta}") # Apply the dynamic theta cuts @@ -525,9 +527,6 @@ def create_irf( irf_hdus.append(hdu_rad_max) - else: - raise ValueError(f"Unknown theta-cut type '{cut_type_theta}'.") - # Create an effective-area HDU logger.info("\nCreating an effective-area HDU...")