From df47f71c23fcd3ed4d8fb5ea8b983452396a6716 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 2 Sep 2024 18:11:38 +0200 Subject: [PATCH] Add more tests for irf makers --- src/ctapipe/irf/tests/test_irfs.py | 145 ++++++++++++++++++++++++- src/ctapipe/irf/tests/test_optimize.py | 18 +-- 2 files changed, 151 insertions(+), 12 deletions(-) diff --git a/src/ctapipe/irf/tests/test_irfs.py b/src/ctapipe/irf/tests/test_irfs.py index ee8d696bd83..394149fb824 100644 --- a/src/ctapipe/irf/tests/test_irfs.py +++ b/src/ctapipe/irf/tests/test_irfs.py @@ -2,8 +2,9 @@ import numpy as np import pytest from astropy.table import QTable, vstack +from pyirf.simulations import SimulatedEventsInfo -from ctapipe.irf import BackgroundRate2dMaker, EventPreProcessor +from ctapipe.irf import EventPreProcessor @pytest.fixture(scope="session") @@ -21,14 +22,150 @@ def irf_events_table(): e_tab = QTable(data=empty.T, names=tab.colnames, units=units) # Setting values following pyirf test in pyirf/irf/tests/test_background.py e_tab["reco_energy"] = np.append(np.full(N1, 1), np.full(N2, 2)) * u.TeV - e_tab["reco_source_fov_offset"] = (np.zeros(N) * u.deg,) + e_tab["true_energy"] = np.append(np.full(N1, 0.9), np.full(N2, 2.1)) * u.TeV + e_tab["reco_source_fov_offset"] = ( + np.append(np.full(N1, 0.1), np.full(N2, 0.05)) * u.deg + ) + e_tab["true_source_fov_offset"] = ( + np.append(np.full(N1, 0.11), np.full(N2, 0.04)) * u.deg + ) ev = vstack([e_tab, tab], join_type="exact", metadata_conflicts="silent") return ev def test_make_2d_bkg(irf_events_table): - bkgMkr = BackgroundRate2dMaker() + from ctapipe.irf import BackgroundRate2dMaker + + bkgMkr = BackgroundRate2dMaker( + fov_offset_n_bins=3, + fov_offset_max=3 * u.deg, + reco_energy_n_bins_per_decade=7, + reco_energy_max=155 * u.TeV, + ) bkg_hdu = bkgMkr.make_bkg_hdu(events=irf_events_table, obs_time=1 * u.s) - assert bkg_hdu.data["BKG"].shape == (1, 1, 20) + # min 7 bins per decade between 0.015 TeV and 155 TeV -> 7 * 4 + 1 = 29 bins + assert bkg_hdu.data["BKG"].shape == (1, 3, 29) + + for col, val in zip(["THETA_LO", "ENERG_LO"], [0 * u.deg, 0.015 * u.TeV]): + assert u.isclose( + u.Quantity(bkg_hdu.data[col][0][0], bkg_hdu.columns[col].unit), val + ) + + for col, val in zip(["THETA_HI", "ENERG_HI"], [3 * u.deg, 155 * u.TeV]): + assert u.isclose( + u.Quantity(bkg_hdu.data[col][0][-1], bkg_hdu.columns[col].unit), val + ) + + +def test_make_2d_energy_migration(irf_events_table): + from ctapipe.irf import EnergyMigration2dMaker + + migMkr = EnergyMigration2dMaker( + fov_offset_n_bins=3, + fov_offset_max=3 * u.deg, + true_energy_n_bins_per_decade=7, + true_energy_max=155 * u.TeV, + energy_migration_n_bins=20, + energy_migration_min=0.1, + energy_migration_max=10, + ) + mig_hdu = migMkr.make_edisp_hdu(events=irf_events_table, point_like=False) + # min 7 bins per decade between 0.015 TeV and 155 TeV -> 7 * 4 + 1 = 29 bins + assert mig_hdu.data["MATRIX"].shape == (1, 3, 20, 29) + + for col, val in zip( + ["THETA_LO", "ENERG_LO", "MIGRA_LO"], [0 * u.deg, 0.015 * u.TeV, 0.1] + ): + assert u.isclose( + u.Quantity(mig_hdu.data[col][0][0], mig_hdu.columns[col].unit), val + ) + + for col, val in zip( + ["THETA_HI", "ENERG_HI", "MIGRA_HI"], [3 * u.deg, 155 * u.TeV, 10] + ): + assert u.isclose( + u.Quantity(mig_hdu.data[col][0][-1], mig_hdu.columns[col].unit), val + ) + + +def test_make_2d_eff_area(irf_events_table): + from ctapipe.irf import EffectiveArea2dMaker + + effAreaMkr = EffectiveArea2dMaker( + fov_offset_n_bins=3, + fov_offset_max=3 * u.deg, + true_energy_n_bins_per_decade=7, + true_energy_max=155 * u.TeV, + ) + sim_info = SimulatedEventsInfo( + n_showers=3000, + energy_min=0.01 * u.TeV, + energy_max=10 * u.TeV, + max_impact=1000 * u.m, + spectral_index=-1.9, + viewcone_min=0 * u.deg, + viewcone_max=10 * u.deg, + ) + eff_area_hdu = effAreaMkr.make_aeff_hdu( + events=irf_events_table, + point_like=False, + signal_is_point_like=False, + sim_info=sim_info, + ) + # min 7 bins per decade between 0.015 TeV and 155 TeV -> 7 * 4 + 1 = 29 bins + assert eff_area_hdu.data["EFFAREA"].shape == (1, 3, 29) + + for col, val in zip(["THETA_LO", "ENERG_LO"], [0 * u.deg, 0.015 * u.TeV]): + assert u.isclose( + u.Quantity(eff_area_hdu.data[col][0][0], eff_area_hdu.columns[col].unit), + val, + ) + + for col, val in zip(["THETA_HI", "ENERG_HI"], [3 * u.deg, 155 * u.TeV]): + assert u.isclose( + u.Quantity(eff_area_hdu.data[col][0][-1], eff_area_hdu.columns[col].unit), + val, + ) + + # point like data -> only 1 fov offset bin + eff_area_hdu = effAreaMkr.make_aeff_hdu( + events=irf_events_table, + point_like=False, + signal_is_point_like=True, + sim_info=sim_info, + ) + assert eff_area_hdu.data["EFFAREA"].shape == (1, 1, 29) + + +def test_make_3d_psf(irf_events_table): + from ctapipe.irf import Psf3dMaker + + psfMkr = Psf3dMaker( + fov_offset_n_bins=3, + fov_offset_max=3 * u.deg, + true_energy_n_bins_per_decade=7, + true_energy_max=155 * u.TeV, + source_offset_n_bins=110, + source_offset_max=2 * u.deg, + ) + psf_hdu = psfMkr.make_psf_hdu(events=irf_events_table) + # min 7 bins per decade between 0.015 TeV and 155 TeV -> 7 * 4 + 1 = 29 bins + assert psf_hdu.data["RPSF"].shape == (1, 110, 3, 29) + + for col, val in zip( + ["THETA_LO", "ENERG_LO", "RAD_LO"], [0 * u.deg, 0.015 * u.TeV, 0 * u.deg] + ): + assert u.isclose( + u.Quantity(psf_hdu.data[col][0][0], psf_hdu.columns[col].unit), + val, + ) + + for col, val in zip( + ["THETA_HI", "ENERG_HI", "RAD_HI"], [3 * u.deg, 155 * u.TeV, 2 * u.deg] + ): + assert u.isclose( + u.Quantity(psf_hdu.data[col][0][-1], psf_hdu.columns[col].unit), + val, + ) diff --git a/src/ctapipe/irf/tests/test_optimize.py b/src/ctapipe/irf/tests/test_optimize.py index e554dec343b..ab744cf7843 100644 --- a/src/ctapipe/irf/tests/test_optimize.py +++ b/src/ctapipe/irf/tests/test_optimize.py @@ -51,10 +51,11 @@ def test_optimization_result_store(tmp_path, irf_events_loader_test_config): def test_gh_percentile_cut_calculator(): from ctapipe.irf import GhPercentileCutCalculator - calc = GhPercentileCutCalculator() - calc.target_percentile = 75 - calc.min_counts = 1 - calc.smoothing = -1 + calc = GhPercentileCutCalculator( + target_percentile=75, + min_counts=1, + smoothing=-1, + ) cuts = calc.calculate_gh_cut( gammaness=np.array([0.1, 0.6, 0.45, 0.98, 0.32, 0.95, 0.25, 0.87]), reco_energy=[0.17, 0.36, 0.47, 0.22, 1.2, 5, 4.2, 9.1] * u.TeV, @@ -69,10 +70,11 @@ def test_gh_percentile_cut_calculator(): def test_theta_percentile_cut_calculator(): from ctapipe.irf import ThetaPercentileCutCalculator - calc = ThetaPercentileCutCalculator() - calc.target_percentile = 75 - calc.min_counts = 1 - calc.smoothing = -1 + calc = ThetaPercentileCutCalculator( + target_percentile=75, + min_counts=1, + smoothing=-1, + ) cuts = calc.calculate_theta_cut( theta=[0.1, 0.07, 0.21, 0.4, 0.03, 0.08, 0.11, 0.18] * u.deg, reco_energy=[0.17, 0.36, 0.47, 0.22, 1.2, 5, 4.2, 9.1] * u.TeV,