From def7f65a286d34498290638ed5554e5d16bd85b4 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Sun, 17 May 2026 18:17:05 +0200 Subject: [PATCH 1/7] Correct limits --- .../production_derive_corsika_limits.py | 45 +++--- .../derive_corsika_limits.py | 142 +++++++++++++----- src/simtools/sim_events/histograms.py | 3 + ...n_derive_corsika_limits_hdf5_db_arrays.yml | 5 +- 4 files changed, 131 insertions(+), 64 deletions(-) diff --git a/src/simtools/applications/production_derive_corsika_limits.py b/src/simtools/applications/production_derive_corsika_limits.py index 615f79d43f..16a6eca32e 100644 --- a/src/simtools/applications/production_derive_corsika_limits.py +++ b/src/simtools/applications/production_derive_corsika_limits.py @@ -12,8 +12,7 @@ Broad-range simulations in this context are simulation sets generated with wide-ranging definitions for above parameters. -Limits are computed based on a configurable maximum event loss fraction -and a minimum number of events lost after cuts. +Limits are computed from configurable per-axis allowed-loss settings. Results are provided as a table with the following columns: +---------------------+-----------+--------+-----------------------------------------------+ @@ -58,11 +57,11 @@ multi-production processing. telescope_ids (str, optional) Custom array layout file containing telescope IDs. -loss_fraction (float, required) - Maximum event-loss fraction for limit computation. -loss_min_events (int, optional) - Minimum number of events that must be lost after applying a limit. - Default: 10. +allowed_losses (str, required, repeatable) + Per-axis allowed-loss tuple in the form + ``axis,fraction,min_events``. + Use once per axis (energy, core_distance, angular_distance), or use ``all`` + to set all axes and optionally override selected axes with additional entries. plot_histograms (bool, optional) Plot histograms of the event data. output_file (str, optional) @@ -83,8 +82,9 @@ simtools-production-derive-corsika-limits \\ --event_data_file event_dat_file.hdf5 \\ --array_layout_name alpha,beta \\ - --loss_fraction 1e-6 \\ - --loss_min_events 10 \\ + --allowed_losses energy,1e-6,10 \ + --allowed_losses core_distance,1e-6,10 \ + --allowed_losses angular_distance,1e-6,10 \ --plot_histograms \\ --output_file corsika_simulation_limits.ecsv @@ -95,8 +95,7 @@ simtools-production-derive-corsika-limits \\ --event_data_file event_dat_file.hdf5 \\ --telescope_ids path/to/telescope_configs.yaml \\ - --loss_fraction 1e-6 \\ - --loss_min_events 10 \\ + --allowed_losses all,1e-6,10 \ --plot_histograms \\ --output_file corsika_simulation_limits.ecsv @@ -108,8 +107,7 @@ --event_data_file pattern_1_*.hdf5 \\ --event_data_file pattern_2_*.hdf5 \\ --array_layout_name alpha \\ - --loss_fraction 1e-6 \\ - --loss_min_events 10 \\ + --allowed_losses all,1e-6,10 \ --plot_histograms \\ --n_workers 4 \\ --output_file corsika_simulation_limits.ecsv @@ -140,17 +138,18 @@ def _add_arguments(parser): required=True, ) parser.add_argument( - "--loss_fraction", - type=float, + "--allowed_losses", + type=str, required=True, - help="Maximum event-loss fraction for limit computation.", - ) - parser.add_argument( - "--loss_min_events", - type=int, - required=False, - default=10, - help="Minimum number of events that must be lost after applying a limit.", + nargs="+", + action="extend", + metavar="AXIS,FRACTION,MIN_EVENTS", + help=( + "Per-axis allowed losses as axis,fraction,min_events. Repeat this argument " + "for each axis using canonical names energy, core_distance, angular_distance, " + "or use all. " + "Example: --allowed_losses energy,1e-6,10" + ), ) parser.add_argument( "--plot_histograms", diff --git a/src/simtools/production_configuration/derive_corsika_limits.py b/src/simtools/production_configuration/derive_corsika_limits.py index cd9979d5e7..fabee5aebe 100644 --- a/src/simtools/production_configuration/derive_corsika_limits.py +++ b/src/simtools/production_configuration/derive_corsika_limits.py @@ -21,6 +21,7 @@ _logger = logging.getLogger(__name__) FILE_INFO_KEYS = ("primary_particle", "zenith", "azimuth", "nsb_level") +LOSS_AXES = ("energy", "core_distance", "angular_distance") RESULT_COLUMNS = [ "production_index", "event_data_file", @@ -122,8 +123,7 @@ def _execute_production_job(job_spec): production_pattern = job_spec["production_pattern"] array_name = job_spec["array_name"] telescope_ids = job_spec["telescope_ids"] - loss_fraction = job_spec["loss_fraction"] - loss_min_events = job_spec.get("loss_min_events", 10) + allowed_losses = job_spec["allowed_losses"] plot_histograms = job_spec["plot_histograms"] output_subdir = job_spec.get("output_subdir") differential_loss_bins_per_decade = job_spec.get("differential_loss_bins_per_decade", 0) @@ -137,8 +137,7 @@ def _execute_production_job(job_spec): production_pattern, array_name, telescope_ids, - loss_fraction, - loss_min_events, + allowed_losses, plot_histograms, output_subdir=output_subdir, differential_loss_bins_per_decade=differential_loss_bins_per_decade, @@ -175,6 +174,72 @@ def _resolve_telescope_configs(args_dict): ) +def _parse_allowed_losses(allowed_losses_args): + """ + Parse repeatable --allowed_losses values into per-axis settings. + + Parameters + ---------- + allowed_losses_args : list[str] + List of values in the form "axis,fraction,min_events". + + Returns + ------- + dict + Mapping of axis name to dict with keys "loss_fraction" and "loss_min_events". + """ + if not allowed_losses_args: + raise ValueError( + "No allowed-loss configuration provided. Use --allowed_losses axis,fraction,min_events" + ) + + parsed = {} + for raw_value in allowed_losses_args: + parts = [part.strip() for part in raw_value.split(",")] + if len(parts) != 3: + raise ValueError( + f"Invalid --allowed_losses value '{raw_value}'. " + "Expected format: axis,fraction,min_events" + ) + + axis_raw, fraction_raw, min_events_raw = parts + try: + fraction = float(fraction_raw) + min_events = int(min_events_raw) + except (TypeError, ValueError) as exc: + raise ValueError( + f"Invalid --allowed_losses value '{raw_value}': " + "fraction must be float and min_events must be int" + ) from exc + + axis_name = axis_raw.strip().lower() + if axis_name == "all": + for axis_name in LOSS_AXES: + parsed[axis_name] = { + "loss_fraction": fraction, + "loss_min_events": min_events, + } + continue + + if axis_name not in LOSS_AXES: + raise ValueError( + "Invalid axis for --allowed_losses. Allowed axes: " + "energy, core_distance, angular_distance, all." + ) + parsed[axis_name] = { + "loss_fraction": fraction, + "loss_min_events": min_events, + } + + missing_axes = [axis_name for axis_name in LOSS_AXES if axis_name not in parsed] + if missing_axes: + raise ValueError( + f"Missing --allowed_losses entries for axis/axes: {', '.join(missing_axes)}" + ) + + return parsed + + def _build_production_subdirectories(production_patterns, output_dir, is_multi_production): """Build and create per-production output subdirectories when needed.""" if not is_multi_production: @@ -205,6 +270,7 @@ def generate_corsika_limits_grid(args_dict): Dictionary containing command line arguments. """ production_patterns = _normalize_event_data_file(args_dict["event_data_file"]) + allowed_losses = _parse_allowed_losses(args_dict.get("allowed_losses")) differential_loss_bins_per_decade = int(args_dict.get("differential_loss_bins_per_decade", 0)) n_productions = len(production_patterns) is_multi_production = n_productions > 1 @@ -236,8 +302,7 @@ def generate_corsika_limits_grid(args_dict): "production_pattern": production_pattern, "array_name": array_name, "telescope_ids": telescope_ids, - "loss_fraction": args_dict["loss_fraction"], - "loss_min_events": int(args_dict.get("loss_min_events", 10)), + "allowed_losses": allowed_losses, "plot_histograms": args_dict["plot_histograms"], "output_subdir": output_subdir, "differential_loss_bins_per_decade": differential_loss_bins_per_decade, @@ -252,15 +317,14 @@ def generate_corsika_limits_grid(args_dict): max_workers=n_workers, ) - write_results(results, args_dict) + write_results(results, args_dict, allowed_losses) def _process_file( file_path, array_name, telescope_ids, - loss_fraction, - loss_min_events=10, + allowed_losses, plot_histograms=False, output_subdir=None, differential_loss_bins_per_decade=0, @@ -278,10 +342,8 @@ def _process_file( Name of the telescope array configuration. telescope_ids : list[str] List of telescope IDs (array-element names) to filter the events. - loss_fraction : float - Fraction of events to be lost. - loss_min_events : int - Minimum number of events to be lost after applying a derived limit. + allowed_losses : dict + Per-axis loss settings for energy/core_distance/angular_distance. plot_histograms : bool Whether to plot histograms. output_subdir : Path or None, optional @@ -306,29 +368,28 @@ def _process_file( limits = { "lower_energy_limit": compute_lower_energy_limit( histograms, - loss_fraction, - 0, # No minimum event loss applied for energy limit (not differential) + allowed_losses["energy"]["loss_fraction"], + allowed_losses["energy"]["loss_min_events"], ), } if differential_loss_bins_per_decade > 0: limits.update( _compute_differential_limits( histograms, - loss_fraction, - loss_min_events, + allowed_losses, differential_loss_bins_per_decade, ) ) else: limits["upper_radius_limit"] = compute_upper_radius_limit( histograms, - loss_fraction, - loss_min_events, + allowed_losses["core_distance"]["loss_fraction"], + allowed_losses["core_distance"]["loss_min_events"], ) limits["viewcone_radius"] = compute_viewcone( histograms, - loss_fraction, - loss_min_events, + allowed_losses["angular_distance"]["loss_fraction"], + allowed_losses["angular_distance"]["loss_min_events"], ) limits.update({key: histograms.file_info.get(key) for key in FILE_INFO_KEYS}) @@ -345,7 +406,11 @@ def _process_file( return limits -def _compute_differential_limits(histograms, loss_fraction, loss_min_events, bins_per_decade): +def _compute_differential_limits( + histograms, + allowed_losses, + bins_per_decade, +): """Compute core and viewcone limits per energy bin and return max limits.""" low = int(np.floor(np.log10(np.min(histograms.energy_bins)))) high = int(np.ceil(np.log10(np.max(histograms.energy_bins)))) @@ -356,8 +421,7 @@ def _compute_differential_limits(histograms, loss_fraction, loss_min_events, bin histograms.core_distance_bins, histograms.energy_bins, diff_e_bins, - loss_fraction, - loss_min_events, + allowed_losses["core_distance"], "core_scatter", "m", ) @@ -366,8 +430,7 @@ def _compute_differential_limits(histograms, loss_fraction, loss_min_events, bin histograms.view_cone_bins, histograms.energy_bins, diff_e_bins, - loss_fraction, - loss_min_events, + allowed_losses["angular_distance"], "viewcone", "deg", ) @@ -401,8 +464,7 @@ def _differential_upper_limits( x_bins, y_bins, diff_e_bins, - loss_fraction, - loss_min_events, + allowed_loss, name, unit, ): @@ -420,8 +482,8 @@ def _differential_upper_limits( limit = _compute_limits( projected, x_bins, - loss_fraction, - loss_min_events, + allowed_loss["loss_fraction"], + allowed_loss["loss_min_events"], limit_type="upper", ) keep = np.searchsorted(x_bins, limit, side="left") @@ -437,7 +499,7 @@ def _differential_upper_limits( ) -def write_results(results, args_dict): +def write_results(results, args_dict, allowed_losses): """ Write the computed limits as astropy table to file. @@ -450,8 +512,7 @@ def write_results(results, args_dict): """ table = _create_results_table( results, - args_dict["loss_fraction"], - args_dict.get("loss_min_events", 10), + allowed_losses, ) output_dir = io_handler.IOHandler().get_output_directory() @@ -463,7 +524,7 @@ def write_results(results, args_dict): MetadataCollector.dump(args_dict, output_file) -def _create_results_table(results, loss_fraction, loss_min_events=10): +def _create_results_table(results, allowed_losses): """ Convert list of simulation results to an astropy Table with metadata. @@ -473,10 +534,8 @@ def _create_results_table(results, loss_fraction, loss_min_events=10): ---------- results : list[dict] Computed limits per file and telescope configuration. - loss_fraction : float - Fraction of lost events (added to metadata). - loss_min_events : int, optional - Minimum number of events to keep after applying a limit. + allowed_losses : dict + Per-axis loss settings added to metadata. Returns ------- @@ -498,10 +557,13 @@ def _create_results_table(results, loss_fraction, loss_min_events=10): { "created": datetime.datetime.now().isoformat(), "description": "Lookup table for CORSIKA limits computed from simulations.", - "loss_fraction": loss_fraction, - "loss_min_events": int(loss_min_events), } ) + for axis_name in LOSS_AXES: + table.meta[f"loss_fraction_{axis_name}"] = allowed_losses[axis_name]["loss_fraction"] + table.meta[f"loss_min_events_{axis_name}"] = int( + allowed_losses[axis_name]["loss_min_events"] + ) return table diff --git a/src/simtools/sim_events/histograms.py b/src/simtools/sim_events/histograms.py index 226896bd1b..e640782c8a 100644 --- a/src/simtools/sim_events/histograms.py +++ b/src/simtools/sim_events/histograms.py @@ -200,12 +200,14 @@ def _define_histograms(self, event_data, triggered_data, shower_data): "event_data": event_data, "bin_edges": self.core_distance_bins, "axis_titles": ["Core Distance (m)", event_count_axis_title], + "plot_scales": {"x": "log", "y": "log"}, }, "angular_distance": { "event_data_column": "angular_distance", "event_data": triggered_data, "bin_edges": self.view_cone_bins, "axis_titles": ["Angular Distance (deg)", event_count_axis_title], + "plot_scales": {"x": "log", "y": "log"}, }, "x_core_shower_vs_y_core_shower": { "event_data_column": ("x_core_shower", "y_core_shower"), @@ -274,6 +276,7 @@ def get_histogram_definition( "1d": is_1d, "bin_edges": bin_edges, "title": title, + "title_fontsize": "x-small", "axis_titles": axis_titles, "suffix": suffix, "plot_scales": plot_scales, diff --git a/tests/integration_tests/config/production_derive_corsika_limits_hdf5_db_arrays.yml b/tests/integration_tests/config/production_derive_corsika_limits_hdf5_db_arrays.yml index 7d588aa8f2..5edd77be5c 100644 --- a/tests/integration_tests/config/production_derive_corsika_limits_hdf5_db_arrays.yml +++ b/tests/integration_tests/config/production_derive_corsika_limits_hdf5_db_arrays.yml @@ -11,7 +11,10 @@ applications: - LSTN-01 - CTAO-North-Alpha event_data_file: tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_*.hdf5 - loss_fraction: 0.001 + allowed_losses: + - energy,0.0001,10 + - core_distance,0.001,10 + - angular_distance,0.001,10 model_version: 6.0.2 output_file: corsika_simulation_limits.ecsv output_path: simtools-output From bfbc52e1ebb74549d5e801777d650899a1544484 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Sun, 17 May 2026 18:29:49 +0200 Subject: [PATCH 2/7] tests --- .../derive_corsika_limits.py | 15 +- .../test_derive_corsika_limits.py | 169 +++++++++++++++--- 2 files changed, 157 insertions(+), 27 deletions(-) diff --git a/src/simtools/production_configuration/derive_corsika_limits.py b/src/simtools/production_configuration/derive_corsika_limits.py index fabee5aebe..4fe54e4e86 100644 --- a/src/simtools/production_configuration/derive_corsika_limits.py +++ b/src/simtools/production_configuration/derive_corsika_limits.py @@ -436,7 +436,7 @@ def _compute_differential_limits( ) upper_radius_limit = _is_close( - core_max * u.m, + _core_distance_ground_to_shower(core_max * u.m, histograms.file_info.get("zenith")), histograms.file_info["core_scatter_max"].to("m") if "core_scatter_max" in histograms.file_info else None, @@ -697,6 +697,13 @@ def _is_close(value, reference, warning_text): return value +def _core_distance_ground_to_shower(core_distance, zenith): + """Convert core distance from ground to shower coordinates using zenith.""" + if zenith is None: + return core_distance + return core_distance * np.cos(zenith.to("rad").value) + + def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): """ Compute the upper radial distance based on the event loss fraction. @@ -715,7 +722,7 @@ def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): astropy.units.Quantity Upper radial distance in m. """ - radius_limit = ( + radius_limit_ground = ( _compute_limits( histograms.histograms["core_distance"]["histogram"], histograms.core_distance_bins, @@ -725,6 +732,10 @@ def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): ) * u.m ) + radius_limit = _core_distance_ground_to_shower( + radius_limit_ground, + histograms.file_info.get("zenith"), + ) return _is_close( radius_limit, histograms.file_info["core_scatter_max"].to("m") diff --git a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py index 7bfa9ebf63..ce11890cf2 100644 --- a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py @@ -17,6 +17,11 @@ ) COMPUTE_VIEWCONE_PATH = "simtools.production_configuration.derive_corsika_limits.compute_viewcone" MOCK_FILE_PATH = "mock_file.fits" +DEFAULT_ALLOWED_LOSSES = { + "energy": {"loss_fraction": 0.2, "loss_min_events": 10}, + "core_distance": {"loss_fraction": 0.2, "loss_min_events": 10}, + "angular_distance": {"loss_fraction": 0.2, "loss_min_events": 10}, +} def _pool_result( @@ -56,7 +61,7 @@ def test_process_file_passes_event_data_patterns_through(mocker): "input/*.h5", "array_name", [1, 2], - 0.2, + DEFAULT_ALLOWED_LOSSES, plot_histograms=False, ) @@ -75,7 +80,7 @@ def test_write_results(mocker, mock_args_dict, mock_results, tmp_test_directory) mock_dump = mocker.patch("simtools.data_model.metadata_collector.MetadataCollector.dump") - derive_corsika_limits.write_results(mock_results, mock_args_dict) + derive_corsika_limits.write_results(mock_results, mock_args_dict, DEFAULT_ALLOWED_LOSSES) # Verify metadata was written mock_dump.assert_called_once() @@ -85,14 +90,19 @@ def test_write_results(mocker, mock_args_dict, mock_results, tmp_test_directory) def test_create_results_table(mock_results): """Test _create_results_table function.""" - table = derive_corsika_limits._create_results_table(mock_results, loss_fraction=0.2) + table = derive_corsika_limits._create_results_table(mock_results, DEFAULT_ALLOWED_LOSSES) table.info() assert isinstance(table, Table) assert len(table) == 1 assert "zenith" in table.colnames assert table["zenith"].unit == u.deg - assert table.meta["loss_fraction"] == pytest.approx(0.2) + assert table.meta["loss_fraction_energy"] == pytest.approx(0.2) + assert table.meta["loss_min_events_energy"] == 10 + assert table.meta["loss_fraction_core_distance"] == pytest.approx(0.2) + assert table.meta["loss_min_events_core_distance"] == 10 + assert table.meta["loss_fraction_angular_distance"] == pytest.approx(0.2) + assert table.meta["loss_min_events_angular_distance"] == 10 assert isinstance(table.meta["created"], str) assert "description" in table.meta @@ -405,7 +415,7 @@ def test_process_file_with_mocked_histograms(mocker): file_path=MOCK_FILE_PATH, array_name="MockArray", telescope_ids=[1, 2], - loss_fraction=0.2, + allowed_losses=DEFAULT_ALLOWED_LOSSES, plot_histograms=False, ) @@ -426,7 +436,7 @@ def test_process_file_with_mocked_histograms(mocker): energy_bins_per_decade=10, ) mock_histograms.fill.assert_called_once() - mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2, 0) + mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2, 10) mock_compute_upper_radius_limit.assert_called_once_with(mock_histograms, 0.2, 10) mock_compute_viewcone.assert_called_once_with(mock_histograms, 0.2, 10) @@ -462,7 +472,7 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): file_path=MOCK_FILE_PATH, array_name="MockArray", telescope_ids=[1, 2], - loss_fraction=0.2, + allowed_losses=DEFAULT_ALLOWED_LOSSES, plot_histograms=False, differential_loss_bins_per_decade=6, ) @@ -473,10 +483,10 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): assert result["core_vs_energy_curve"] == {"x": [100.0, 120.0], "y": [0.1, 1.0]} assert result["angular_distance_vs_energy_curve"] == {"x": [2.5, 3.0], "y": [0.1, 1.0]} - mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2, 0) + mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2, 10) mock_compute_upper_radius_limit.assert_not_called() mock_compute_viewcone.assert_not_called() - mock_differential.assert_called_once_with(mock_histograms, 0.2, 10, 6) + mock_differential.assert_called_once_with(mock_histograms, DEFAULT_ALLOWED_LOSSES, 6) @pytest.mark.parametrize( @@ -516,15 +526,17 @@ def test_compute_differential_limits( side_effect=[125.0 * u.m, 3.25 * u.deg], ) - derive_corsika_limits._compute_differential_limits(histograms, 0.2, 10, 2) + derive_corsika_limits._compute_differential_limits(histograms, DEFAULT_ALLOWED_LOSSES, 2) expected_diff_bins = np.logspace(0, 1, 3) np.testing.assert_allclose(mock_diff_limits.call_args_list[0].args[3], expected_diff_bins) np.testing.assert_allclose(mock_diff_limits.call_args_list[1].args[3], expected_diff_bins) assert mock_diff_limits.call_args_list[0].args[0] == "core-hist" - assert mock_diff_limits.call_args_list[0].args[6:] == ("core_scatter", "m") + assert mock_diff_limits.call_args_list[0].args[5:] == ("core_scatter", "m") assert mock_diff_limits.call_args_list[1].args[0] == "viewcone-hist" - assert mock_diff_limits.call_args_list[1].args[6:] == ("viewcone", "deg") + assert mock_diff_limits.call_args_list[1].args[5:] == ("viewcone", "deg") + assert mock_diff_limits.call_args_list[0].args[4] == DEFAULT_ALLOWED_LOSSES["core_distance"] + assert mock_diff_limits.call_args_list[1].args[4] == DEFAULT_ALLOWED_LOSSES["angular_distance"] assert mock_is_close.call_args_list[0].args[0].value == pytest.approx(120.0) assert mock_is_close.call_args_list[0].args[1] == expected_core_scatter_max @@ -555,7 +567,7 @@ def test_process_file_passes_energy_bins_per_decade_to_histograms(mocker): file_path=MOCK_FILE_PATH, array_name="MockArray", telescope_ids=[1, 2], - loss_fraction=0.2, + allowed_losses=DEFAULT_ALLOWED_LOSSES, plot_histograms=False, differential_loss_bins_per_decade=6, ) @@ -581,8 +593,7 @@ def test_differential_upper_limits(mocker): x_bins=np.array([0.0, 1.0, 2.0, 3.0]), y_bins=np.array([1.0, 2.0, 4.0]), diff_e_bins=np.array([1.0, 2.0, 2.5, 3.0]), - loss_fraction=0.2, - loss_min_events=10, + allowed_loss=DEFAULT_ALLOWED_LOSSES["core_distance"], name="core_scatter", unit="m", ) @@ -611,8 +622,7 @@ def test_differential_upper_limits_falls_back_to_last_bin_edge(mocker): x_bins=np.array([0.0, 1.0, 2.0, 3.0]), y_bins=np.array([1.0, 2.0, 4.0]), diff_e_bins=np.array([1.0, 2.0, 3.0]), - loss_fraction=0.2, - loss_min_events=10, + allowed_loss=DEFAULT_ALLOWED_LOSSES["angular_distance"], name="viewcone", unit="deg", ) @@ -659,7 +669,7 @@ def test_process_file_with_plot_histograms(mocker, tmp_test_directory): file_path=MOCK_FILE_PATH, array_name="MockArray", telescope_ids=[1, 2], - loss_fraction=0.2, + allowed_losses=DEFAULT_ALLOWED_LOSSES, plot_histograms=True, ) @@ -738,6 +748,103 @@ def test_get_production_directory_name_appends_uuid_on_collision(mocker): mock_uuid.assert_called_once() +def test_parse_allowed_losses_explicit_axes(): + """Test _parse_allowed_losses with explicit per-axis entries.""" + result = derive_corsika_limits._parse_allowed_losses( + [ + "energy,1e-6,10", + "core_distance,2e-6,20", + "angular_distance,3e-6,30", + ] + ) + + assert result["energy"]["loss_fraction"] == pytest.approx(1e-6) + assert result["energy"]["loss_min_events"] == 10 + assert result["core_distance"]["loss_fraction"] == pytest.approx(2e-6) + assert result["core_distance"]["loss_min_events"] == 20 + assert result["angular_distance"]["loss_fraction"] == pytest.approx(3e-6) + assert result["angular_distance"]["loss_min_events"] == 30 + + +def test_parse_allowed_losses_all_and_override(): + """Test _parse_allowed_losses supports all plus later axis override.""" + result = derive_corsika_limits._parse_allowed_losses( + [ + "all,1e-6,10", + "energy,5e-7,5", + ] + ) + + assert result["energy"]["loss_fraction"] == pytest.approx(5e-7) + assert result["energy"]["loss_min_events"] == 5 + assert result["core_distance"]["loss_fraction"] == pytest.approx(1e-6) + assert result["core_distance"]["loss_min_events"] == 10 + assert result["angular_distance"]["loss_fraction"] == pytest.approx(1e-6) + assert result["angular_distance"]["loss_min_events"] == 10 + + +def test_parse_allowed_losses_missing_axis_raises(): + """Test _parse_allowed_losses raises when required axes are missing.""" + with pytest.raises(ValueError, match="Missing --allowed_losses entries"): + derive_corsika_limits._parse_allowed_losses( + [ + "energy,1e-6,10", + "core_distance,1e-6,10", + ] + ) + + +def test_parse_allowed_losses_invalid_axis_raises(): + """Test _parse_allowed_losses rejects invalid axis names.""" + with pytest.raises(ValueError, match="Invalid axis for --allowed_losses"): + derive_corsika_limits._parse_allowed_losses( + [ + "energy,1e-6,10", + "core_distance,1e-6,10", + "viewcone,1e-6,10", + ] + ) + + +def test_build_production_subdirectories_non_multi_returns_empty(tmp_test_directory): + """Test _build_production_subdirectories returns empty dict for single production.""" + result = derive_corsika_limits._build_production_subdirectories( + ["pattern_1_*.hdf5"], + tmp_test_directory, + is_multi_production=False, + ) + assert result == {} + + +def test_build_production_subdirectories_creates_dirs(tmp_test_directory): + """Test _build_production_subdirectories creates per-production directories.""" + patterns = ["pattern_1_*.hdf5", "pattern_2_*.hdf5"] + result = derive_corsika_limits._build_production_subdirectories( + patterns, + tmp_test_directory, + is_multi_production=True, + ) + + assert set(result.keys()) == set(patterns) + for output_subdir in result.values(): + assert output_subdir.exists() + assert output_subdir.isdir() + + +def test_core_distance_ground_to_shower_converts_with_zenith(): + """Test _core_distance_ground_to_shower scales by cos(zenith).""" + result = derive_corsika_limits._core_distance_ground_to_shower(100.0 * u.m, 60.0 * u.deg) + assert result.unit == u.m + assert result.value == pytest.approx(50.0) + + +def test_core_distance_ground_to_shower_no_zenith_returns_input(): + """Test _core_distance_ground_to_shower returns input unchanged without zenith.""" + input_distance = 123.0 * u.m + result = derive_corsika_limits._core_distance_ground_to_shower(input_distance, None) + assert result == input_distance + + def test_execute_production_job_single_job(mocker): """Test _execute_production_job executes one job correctly.""" mock_histograms = mocker.MagicMock() @@ -767,7 +874,7 @@ def test_execute_production_job_single_job(mocker): "production_pattern": "pattern_*.hdf5", "array_name": "LST", "telescope_ids": ["LSTN-01"], - "loss_fraction": 0.2, + "allowed_losses": DEFAULT_ALLOWED_LOSSES, "plot_histograms": False, "output_subdir": None, } @@ -821,7 +928,11 @@ def test_generate_corsika_limits_grid_multi_production(mocker, tmp_test_director args_dict = { "event_data_file": ["pattern_1_*.hdf5", "pattern_2_*.hdf5"], "telescope_ids": "telescope_ids.yml", - "loss_fraction": 0.2, + "allowed_losses": [ + "energy,0.2,10", + "core_distance,0.2,10", + "angular_distance,0.2,10", + ], "plot_histograms": False, "output_file": "test_output.ecsv", "n_workers": 2, @@ -870,7 +981,11 @@ def test_generate_corsika_limits_grid_single_production_uses_pool(mocker, tmp_te args_dict = { "event_data_file": "pattern_*.hdf5", # Single string, not list "telescope_ids": "telescope_ids.yml", - "loss_fraction": 0.2, + "allowed_losses": [ + "energy,0.2,10", + "core_distance,0.2,10", + "angular_distance,0.2,10", + ], "plot_histograms": False, "output_file": "test_output.ecsv", "n_workers": 0, @@ -891,7 +1006,7 @@ def test_create_results_table_with_production_columns(mock_results): res["production_index"] = i res["event_data_file"] = f"pattern_{i}_*.hdf5" - table = derive_corsika_limits._create_results_table(mock_results, loss_fraction=0.2) + table = derive_corsika_limits._create_results_table(mock_results, DEFAULT_ALLOWED_LOSSES) # Should include production-origin columns assert "production_index" in table.colnames @@ -905,7 +1020,7 @@ def test_create_results_table_with_production_columns(mock_results): def test_create_results_table_without_production_columns(mock_results): """Test _create_results_table with missing production metadata values.""" # Results without production metadata (old format) - table = derive_corsika_limits._create_results_table(mock_results, loss_fraction=0.2) + table = derive_corsika_limits._create_results_table(mock_results, DEFAULT_ALLOWED_LOSSES) # Production-origin columns are included and filled with None if missing assert "production_index" in table.colnames @@ -952,7 +1067,7 @@ def test_process_file_with_output_subdir(mocker, tmp_test_directory): file_path=MOCK_FILE_PATH, array_name="MockArray", telescope_ids=["LSTN-01"], - loss_fraction=0.2, + allowed_losses=DEFAULT_ALLOWED_LOSSES, plot_histograms=True, output_subdir=output_subdir, ) @@ -972,7 +1087,11 @@ def mock_args_dict(): "ignore_runtime_environment": False, "event_data_file": "dummy_event_data.h5", "output_file": "corsika_limits.ecsv", - "loss_fraction": 0.2, + "allowed_losses": [ + "energy,0.2,10", + "core_distance,0.2,10", + "angular_distance,0.2,10", + ], "plot_histograms": False, "n_workers": 1, "array_layout_name": None, From d521b4034f001e1d929df3616e2b471acfe856e5 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 18 May 2026 14:35:08 +0200 Subject: [PATCH 3/7] plotting improvement --- src/simtools/sim_events/histograms.py | 6 +- .../plot_simtel_event_histograms.py | 80 +---------- .../test_plot_simtel_event_histograms.py | 128 +----------------- 3 files changed, 8 insertions(+), 206 deletions(-) diff --git a/src/simtools/sim_events/histograms.py b/src/simtools/sim_events/histograms.py index e640782c8a..4d64340bf1 100644 --- a/src/simtools/sim_events/histograms.py +++ b/src/simtools/sim_events/histograms.py @@ -200,14 +200,14 @@ def _define_histograms(self, event_data, triggered_data, shower_data): "event_data": event_data, "bin_edges": self.core_distance_bins, "axis_titles": ["Core Distance (m)", event_count_axis_title], - "plot_scales": {"x": "log", "y": "log"}, + "plot_scales": {"y": "log"}, }, "angular_distance": { "event_data_column": "angular_distance", "event_data": triggered_data, "bin_edges": self.view_cone_bins, "axis_titles": ["Angular Distance (deg)", event_count_axis_title], - "plot_scales": {"x": "log", "y": "log"}, + "plot_scales": {"y": "log"}, }, "x_core_shower_vs_y_core_shower": { "event_data_column": ("x_core_shower", "y_core_shower"), @@ -276,7 +276,7 @@ def get_histogram_definition( "1d": is_1d, "bin_edges": bin_edges, "title": title, - "title_fontsize": "x-small", + "title_fontsize": "xx-small", "axis_titles": axis_titles, "suffix": suffix, "plot_scales": plot_scales, diff --git a/src/simtools/visualization/plot_simtel_event_histograms.py b/src/simtools/visualization/plot_simtel_event_histograms.py index 0e652a3085..c6a6104ea5 100644 --- a/src/simtools/visualization/plot_simtel_event_histograms.py +++ b/src/simtools/visualization/plot_simtel_event_histograms.py @@ -6,12 +6,10 @@ import numpy as np from matplotlib.colors import LogNorm -from simtools.sim_events.histograms import EventDataHistograms - _logger = logging.getLogger(__name__) -def plot(histograms, output_path=None, limits=None, rebin_factor=2, array_name=None): +def plot(histograms, output_path=None, limits=None, array_name=None): """ Plot simtel event histograms. @@ -26,16 +24,13 @@ def plot(histograms, output_path=None, limits=None, rebin_factor=2, array_name=N - "upper_radius_limit": Upper limit for core distance - "lower_energy_limit": Lower limit for energy - "viewcone_radius": Radius for the viewcone - rebin_factor: int, optional - Factor by which to reduce the number of bins in 2D histograms for re-binned plots. - Default is 2 (merge every 2 bins). Set to 0 or 1 to disable re-binning. array_name: str, optional Name of the telescope array configuration. """ _logger.info(f"Plotting histograms written to {output_path}") plots = _generate_plot_configurations(histograms, limits) - _execute_plotting_loop(plots, output_path, rebin_factor, array_name) + _execute_plotting_loop(plots, output_path, array_name) def _get_limits(name, limits): @@ -155,7 +150,7 @@ def _create_2d_plot_config(histogram, name, plot_params, limits): } -def _execute_plotting_loop(plots, output_path, rebin_factor, array_name): +def _execute_plotting_loop(plots, output_path, array_name): """Execute the main plotting loop for all plot configurations.""" for plot_key, plot_args in plots.items(): plot_filename = plot_args.pop("filename") @@ -169,14 +164,7 @@ def _execute_plotting_loop(plots, output_path, rebin_factor, array_name): filename = _build_plot_filename(plot_filename, array_name) output_file = output_path / filename if output_path else None - result = _create_plot(**plot_args, output_file=output_file) - - # Skip re-binned plot if main plot failed - if result is None: - continue - - if _should_create_rebinned_plot(rebin_factor, plot_args, plot_key): - _create_rebinned_plot(plot_args, filename, output_path, rebin_factor) + _create_plot(**plot_args, output_file=output_file) def _build_plot_filename(base_filename, array_name=None): @@ -198,66 +186,6 @@ def _build_plot_filename(base_filename, array_name=None): return f"{base_filename}_{array_name}.png" if array_name else f"{base_filename}.png" -def _should_create_rebinned_plot(rebin_factor, plot_args, plot_key): - """ - Check if a re-binned version of the plot should be created. - - Parameters - ---------- - rebin_factor : int - Factor by which to rebin the energy axis - plot_args : dict - Plot arguments - plot_key : str - Key identifying the plot type - - Returns - ------- - bool - True if a re-binned plot should be created, False otherwise - """ - return ( - rebin_factor > 1 - and plot_args["plot_type"] == "histogram2d" - and plot_key.endswith("_cumulative") - and plot_args.get("plot_params", {}).get("norm") == "linear" - ) - - -def _create_rebinned_plot(plot_args, filename, output_path, rebin_factor): - """ - Create a re-binned version of a 2D histogram plot. - - Parameters - ---------- - plot_args : dict - Plot arguments for the original plot - filename : str - Filename of the original plot - output_path : Path or None - Path to save the plot to, or None - rebin_factor : int - Factor by which to rebin the energy axis - """ - data = plot_args["data"] - bins = plot_args["bins"] - - rebinned_data, rebinned_x_bins, rebinned_y_bins = EventDataHistograms.rebin_2d_histogram( - data, bins[0], bins[1], rebin_factor - ) - - rebinned_plot_args = plot_args.copy() - rebinned_plot_args["data"] = rebinned_data - rebinned_plot_args["bins"] = [rebinned_x_bins, rebinned_y_bins] - - if rebinned_plot_args.get("labels", {}).get("title"): - rebinned_plot_args["labels"]["title"] += f" (Energy rebinned {rebin_factor}x)" - - rebinned_filename = f"{filename.replace('.png', '')}_rebinned.png" - rebinned_output_file = output_path / rebinned_filename if output_path else None - _create_plot(**rebinned_plot_args, output_file=rebinned_output_file) - - def _create_plot( data, bins=None, diff --git a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py index 3765d34681..e7c6304bcf 100644 --- a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py +++ b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py @@ -11,7 +11,6 @@ _build_plot_filename, _create_2d_histogram_plot, _create_plot, - _create_rebinned_plot, _execute_plotting_loop, _get_limits, plot, @@ -29,8 +28,6 @@ PATCH_CONTOUR = f"{MOD}.plt.contour" PATCH_COLORBAR = f"{MOD}.plt.colorbar" PATCH_CREATE_PLOT = f"{MOD}._create_plot" -PATCH_CREATE_REBINNED = f"{MOD}._create_rebinned_plot" -PATCH_REBIN = f"{MOD}.EventDataHistograms.rebin_2d_histogram" PATCH_HAS_DATA = f"{MOD}._has_data" PATCH_BUILD_FILENAME = f"{MOD}._build_plot_filename" @@ -356,79 +353,6 @@ def test_create_plot_early_return_when_no_data(): mock_subplots.assert_not_called() -@pytest.mark.parametrize("output_path_is_none", [True, False]) -def test_create_rebinned_plot(output_path_is_none, tmp_path): - data = np.array([[1, 2], [3, 4]]) - bins = [np.array([0, 1, 2]), np.array([0, 1, 2])] - rebin_factor = 2 - plot_args = { - "data": data, - "bins": bins, - "plot_type": "histogram2d", - "plot_params": {"norm": "linear", "cmap": "viridis"}, - "labels": {"title": "Original Plot"}, - } - filename = "test_plot.png" - output_path = None if output_path_is_none else tmp_path - - rebinned_data = np.array([[10]]) - rebinned_x_bins = np.array([0, 2]) - rebinned_y_bins = np.array([0, 2]) - - with ( - patch( - PATCH_REBIN, - return_value=(rebinned_data, rebinned_x_bins, rebinned_y_bins), - ) as mock_rebin, - patch(PATCH_CREATE_PLOT) as mock_create_plot, - ): - _create_rebinned_plot(plot_args, filename, output_path, rebin_factor) - - mock_rebin.assert_called_once_with(data, bins[0], bins[1], rebin_factor) - mock_create_plot.assert_called_once() - rebinned_plot_args = mock_create_plot.call_args[1] - np.testing.assert_array_equal(rebinned_plot_args["data"], rebinned_data) - if output_path is None: - assert rebinned_plot_args["output_file"] is None - else: - assert rebinned_plot_args["output_file"].name == "test_plot_rebinned.png" - - -def test_should_create_rebinned_plot(): - plot_args = { - "plot_type": "histogram2d", - "plot_params": {"norm": "linear"}, - } - - # Test case: rebin_factor > 1, plot_type is histogram2d, ends with _cumulative, norm is linear - assert plot_simtel_event_histograms._should_create_rebinned_plot( - rebin_factor=2, plot_args=plot_args, plot_key="test_cumulative" - ) - - # Test case: rebin_factor <= 1 - assert not plot_simtel_event_histograms._should_create_rebinned_plot( - rebin_factor=1, plot_args=plot_args, plot_key="test_cumulative" - ) - - # Test case: plot_type is not histogram2d - plot_args["plot_type"] = "histogram" - assert not plot_simtel_event_histograms._should_create_rebinned_plot( - rebin_factor=2, plot_args=plot_args, plot_key="test_cumulative" - ) - - # Test case: plot_key does not end with _cumulative - plot_args["plot_type"] = "histogram2d" - assert not plot_simtel_event_histograms._should_create_rebinned_plot( - rebin_factor=2, plot_args=plot_args, plot_key="test" - ) - - # Test case: norm is not linear - plot_args["plot_params"]["norm"] = "log" - assert not plot_simtel_event_histograms._should_create_rebinned_plot( - rebin_factor=2, plot_args=plot_args, plot_key="test_cumulative" - ) - - def test_build_plot_filename(): # Test without array_name base_filename = "test_plot" @@ -463,17 +387,15 @@ def test_execute_plotting_loop(): }, } output_path = MagicMock() - rebin_factor = 2 array_name = "test_array" with ( patch(PATCH_CREATE_PLOT) as mock_create_plot, - patch(PATCH_CREATE_REBINNED) as mock_create_rebinned_plot, patch(PATCH_BUILD_FILENAME, side_effect=lambda base, array: f"{base}_{array}.png"), ): mock_create_plot.return_value = MagicMock() # Simulate successful plot creation - _execute_plotting_loop(plots, output_path, rebin_factor, array_name) + _execute_plotting_loop(plots, output_path, array_name) # Ensure exactly one plot was created assert mock_create_plot.call_count == 1 @@ -495,54 +417,6 @@ def test_execute_plotting_loop(): # Ensure the second plot was skipped (still only one call) mock_create_plot.assert_called_once() - # For histogram (not 2D cumulative), no rebinned plot expected - mock_create_rebinned_plot.assert_not_called() - - -def test_execute_plotting_loop_rebin_and_failed_plot(): - """Cover branches where a plot returns None and where a rebinned plot is created.""" - plots = { - # Meets all conditions for rebin (2D cumulative, linear norm, rebin_factor > 1) - "plotA_cumulative": { - "data": np.array([[1, 2], [3, 4]]), - "bins": (np.array([0, 1, 2]), np.array([0, 1, 2])), - "plot_type": "histogram2d", - "plot_params": {"norm": "linear"}, - "labels": {"title": "Rebin Test"}, - "scales": {}, - "filename": "plotA_cumulative", - }, - # This plot will return None from _create_plot to exercise the early-continue branch - "plotB": { - "data": np.array([1, 2, 3]), - "bins": np.array([0, 1, 2, 3]), - "plot_type": "histogram", - "plot_params": {"color": "blue"}, - "labels": {"title": "Should Skip"}, - "scales": {}, - "filename": "plotB", - }, - } - output_path = MagicMock() - rebin_factor = 2 - array_name = None - - with ( - patch(PATCH_CREATE_PLOT, side_effect=[MagicMock(), None]) as mock_create_plot, - patch(PATCH_CREATE_REBINNED) as mock_create_rebinned_plot, - patch(PATCH_BUILD_FILENAME, side_effect=lambda base, array: f"{base}.png"), - ): - _execute_plotting_loop(plots, output_path, rebin_factor, array_name) - - # First call produced a figure, second returned None - assert mock_create_plot.call_count == 2 - # Rebinned plot should be created exactly once for plotA_cumulative - mock_create_rebinned_plot.assert_called_once() - args, _ = mock_create_rebinned_plot.call_args - assert args[1] == "plotA_cumulative.png" - assert args[2] is output_path - assert args[3] == rebin_factor - def test_create_2d_plot_config(): histograms = { From bda34631cdd3fea761ca8cbfc718f2540326bf3e Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 18 May 2026 14:37:52 +0200 Subject: [PATCH 4/7] unit test --- .../visualization/test_plot_simtel_event_histograms.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py index e7c6304bcf..8b75b2267e 100644 --- a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py +++ b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py @@ -641,7 +641,6 @@ def mock_histograms(): def test_plot_with_output_path(mock_histograms): output_path = Path("/mock/output/path") limits = {"upper_radius_limit": MagicMock(value=100)} - rebin_factor = 2 array_name = "test_array" with ( @@ -654,19 +653,17 @@ def test_plot_with_output_path(mock_histograms): histograms=mock_histograms, output_path=output_path, limits=limits, - rebin_factor=rebin_factor, array_name=array_name, ) mock_generate_configs.assert_called_once_with(mock_histograms, limits) mock_execute_loop.assert_called_once_with( - {"mock_plot": "mock_config"}, output_path, rebin_factor, array_name + {"mock_plot": "mock_config"}, output_path, array_name ) def test_plot_without_output_path(mock_histograms): limits = None - rebin_factor = 1 array_name = None with ( @@ -679,14 +676,11 @@ def test_plot_without_output_path(mock_histograms): histograms=mock_histograms, output_path=None, limits=limits, - rebin_factor=rebin_factor, array_name=array_name, ) mock_generate_configs.assert_called_once_with(mock_histograms, limits) - mock_execute_loop.assert_called_once_with( - {"mock_plot": "mock_config"}, None, rebin_factor, array_name - ) + mock_execute_loop.assert_called_once_with({"mock_plot": "mock_config"}, None, array_name) def test_get_axis_title(): From 0411baf54d7859804805591b9e415e60571f916a Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 18 May 2026 14:58:14 +0200 Subject: [PATCH 5/7] ground vs shower coordinates --- .../derive_corsika_limits.py | 50 +++++++---- .../plot_simtel_event_histograms.py | 13 +-- .../production_derive_corsika_limits_fits.yml | 5 +- .../test_derive_corsika_limits.py | 90 +++++++++++++------ 4 files changed, 111 insertions(+), 47 deletions(-) diff --git a/src/simtools/production_configuration/derive_corsika_limits.py b/src/simtools/production_configuration/derive_corsika_limits.py index 4fe54e4e86..7657c3ef9f 100644 --- a/src/simtools/production_configuration/derive_corsika_limits.py +++ b/src/simtools/production_configuration/derive_corsika_limits.py @@ -32,7 +32,8 @@ "azimuth", "nsb_level", "lower_energy_limit", - "upper_radius_limit", + "upper_radius_limit_ground", + "upper_radius_limit_shower", "viewcone_radius", ] @@ -359,9 +360,9 @@ def _process_file( """ histograms = EventDataHistograms( file_path, - array_name=array_name, - telescope_list=telescope_ids, - energy_bins_per_decade=differential_loss_bins_per_decade or 10, + array_name, + telescope_ids, + differential_loss_bins_per_decade or 10, ) histograms.fill() @@ -381,11 +382,12 @@ def _process_file( ) ) else: - limits["upper_radius_limit"] = compute_upper_radius_limit( + radius_limits = compute_upper_radius_limit( histograms, allowed_losses["core_distance"]["loss_fraction"], allowed_losses["core_distance"]["loss_min_events"], ) + limits.update(radius_limits) limits["viewcone_radius"] = compute_viewcone( histograms, allowed_losses["angular_distance"]["loss_fraction"], @@ -435,8 +437,12 @@ def _compute_differential_limits( "deg", ) - upper_radius_limit = _is_close( - _core_distance_ground_to_shower(core_max * u.m, histograms.file_info.get("zenith")), + upper_radius_limit_ground = core_max * u.m + upper_radius_limit_shower = _core_distance_ground_to_shower( + upper_radius_limit_ground, histograms.file_info.get("zenith") + ) + upper_radius_limit_shower = _is_close( + upper_radius_limit_shower, histograms.file_info["core_scatter_max"].to("m") if "core_scatter_max" in histograms.file_info else None, @@ -449,10 +455,16 @@ def _compute_differential_limits( else None, "Upper viewcone limit is equal to the maximum viewcone distance of", ) - _logger.info(f"Differential upper_radius_limit (max over bins): {upper_radius_limit}") + _logger.info( + f"Differential upper_radius_limit_ground (max over bins): {upper_radius_limit_ground}" + ) + _logger.info( + f"Differential upper_radius_limit_shower (max over bins): {upper_radius_limit_shower}" + ) _logger.info(f"Differential viewcone_radius (max over bins): {viewcone_radius}") return { - "upper_radius_limit": upper_radius_limit, + "upper_radius_limit_ground": upper_radius_limit_ground, + "upper_radius_limit_shower": upper_radius_limit_shower, "viewcone_radius": viewcone_radius, "core_vs_energy_curve": {"x": core_x, "y": core_y}, "angular_distance_vs_energy_curve": {"x": vc_x, "y": vc_y}, @@ -589,7 +601,7 @@ def _round_value(key, val): """Round value based on key type.""" if key == "lower_energy_limit": return np.floor(val * 1e3) / 1e3 - if key == "upper_radius_limit": + if key in ("upper_radius_limit_ground", "upper_radius_limit_shower"): return np.ceil(val / 25) * 25 if key == "viewcone_radius": return np.ceil(val / 0.25) * 0.25 @@ -706,7 +718,7 @@ def _core_distance_ground_to_shower(core_distance, zenith): def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): """ - Compute the upper radial distance based on the event loss fraction. + Compute the upper radial distance in both ground and shower coordinates. Parameters ---------- @@ -719,8 +731,10 @@ def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): Returns ------- - astropy.units.Quantity - Upper radial distance in m. + dict + Dictionary containing: + - "upper_radius_limit_ground": Upper radial distance in ground coordinates (m). + - "upper_radius_limit_shower": Upper radial distance in shower coordinates (m). """ radius_limit_ground = ( _compute_limits( @@ -732,17 +746,21 @@ def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): ) * u.m ) - radius_limit = _core_distance_ground_to_shower( + radius_limit_shower = _core_distance_ground_to_shower( radius_limit_ground, histograms.file_info.get("zenith"), ) - return _is_close( - radius_limit, + radius_limit_shower = _is_close( + radius_limit_shower, histograms.file_info["core_scatter_max"].to("m") if "core_scatter_max" in histograms.file_info else None, "Upper radius limit is equal to the maximum core scatter distance of", ) + return { + "upper_radius_limit_ground": radius_limit_ground, + "upper_radius_limit_shower": radius_limit_shower, + } def compute_viewcone(histograms, loss_fraction, loss_min_events=10): diff --git a/src/simtools/visualization/plot_simtel_event_histograms.py b/src/simtools/visualization/plot_simtel_event_histograms.py index c6a6104ea5..74a477d214 100644 --- a/src/simtools/visualization/plot_simtel_event_histograms.py +++ b/src/simtools/visualization/plot_simtel_event_histograms.py @@ -21,9 +21,11 @@ def plot(histograms, output_path=None, limits=None, array_name=None): Directory to save plots. If None, plots will be displayed. limits: dict, optional Dictionary containing limits for plotting. Keys can include: - - "upper_radius_limit": Upper limit for core distance + - "upper_radius_limit_ground": Upper limit for core distance (ground coordinates) + - "upper_radius_limit_shower": Upper limit for core distance (shower coordinates) - "lower_energy_limit": Lower limit for energy - "viewcone_radius": Radius for the viewcone + Note: Plotting always uses ground coordinates (upper_radius_limit_ground). array_name: str, optional Name of the telescope array configuration. """ @@ -38,6 +40,7 @@ def _get_limits(name, limits): Extract limits from the provided dictionary for plotting. Fine tuned to expected histograms to be plotted. + Always uses ground coordinates for upper_radius_limit. """ def _safe_value(limits, key): @@ -46,15 +49,15 @@ def _safe_value(limits, key): mapping = { "energy": {"x": _safe_value(limits, "lower_energy_limit")}, - "core_distance": {"x": _safe_value(limits, "upper_radius_limit")}, + "core_distance": {"x": _safe_value(limits, "upper_radius_limit_ground")}, "angular_distance": {"x": _safe_value(limits, "viewcone_radius")}, "core_vs_energy": { - "x": _safe_value(limits, "upper_radius_limit"), + "x": _safe_value(limits, "upper_radius_limit_ground"), "y": _safe_value(limits, "lower_energy_limit"), "curve": limits.get("core_vs_energy_curve"), }, "core_vs_energy_cumulative": { - "x": _safe_value(limits, "upper_radius_limit"), + "x": _safe_value(limits, "upper_radius_limit_ground"), "y": _safe_value(limits, "lower_energy_limit"), "curve": limits.get("core_vs_energy_curve"), }, @@ -68,7 +71,7 @@ def _safe_value(limits, key): "y": _safe_value(limits, "lower_energy_limit"), "curve": limits.get("angular_distance_vs_energy_curve"), }, - "x_core_shower_vs_y_core_shower": {"r": _safe_value(limits, "upper_radius_limit")}, + "x_core_shower_vs_y_core_shower": {"r": _safe_value(limits, "upper_radius_limit_ground")}, } return mapping.get(name) diff --git a/tests/integration_tests/config/production_derive_corsika_limits_fits.yml b/tests/integration_tests/config/production_derive_corsika_limits_fits.yml index 4e27ffb13e..c277abc10e 100644 --- a/tests/integration_tests/config/production_derive_corsika_limits_fits.yml +++ b/tests/integration_tests/config/production_derive_corsika_limits_fits.yml @@ -3,7 +3,10 @@ applications: - application: simtools-production-derive-corsika-limits configuration: event_data_file: tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_reduced_event_data.fits.gz - loss_fraction: 0.001 + allowed_losses: + - energy,0.0001,10 + - core_distance,0.001,10 + - angular_distance,0.001,10 output_file: corsika_simulation_limits.ecsv output_path: simtools-output plot_histograms: true diff --git a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py index ce11890cf2..9933478f1a 100644 --- a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py @@ -30,7 +30,8 @@ def _pool_result( array_name="LST", telescope_ids=None, lower_energy_limit=0.5 * u.TeV, - upper_radius_limit=400.0 * u.m, + upper_radius_limit_ground=400.0 * u.m, + upper_radius_limit_shower=380.0 * u.m, viewcone_radius=5.0 * u.deg, ): """Build a standard mocked pool result row for grid execution tests.""" @@ -40,7 +41,8 @@ def _pool_result( "array_name": array_name, "telescope_ids": telescope_ids or ["LSTN-01"], "lower_energy_limit": lower_energy_limit, - "upper_radius_limit": upper_radius_limit, + "upper_radius_limit_ground": upper_radius_limit_ground, + "upper_radius_limit_shower": upper_radius_limit_shower, "viewcone_radius": viewcone_radius, "primary_particle": "gamma", "zenith": 20.0 * u.deg, @@ -54,7 +56,13 @@ def test_process_file_passes_event_data_patterns_through(mocker): mock_histograms = mocker.MagicMock() mock_histogram_class = mocker.patch(SIM_EVENTS_HISTOGRAMS_PATH, return_value=mock_histograms) mocker.patch(COMPUTE_LOWER_ENERGY_LIMIT_PATH, return_value=1.0 * u.TeV) - mocker.patch(COMPUTE_UPPER_RADIUS_LIMIT_PATH, return_value=100.0 * u.m) + mocker.patch( + COMPUTE_UPPER_RADIUS_LIMIT_PATH, + return_value={ + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, + }, + ) mocker.patch(COMPUTE_VIEWCONE_PATH, return_value=2.0 * u.deg) derive_corsika_limits._process_file( @@ -115,11 +123,13 @@ def test_round_value(): assert derive_corsika_limits._round_value("lower_energy_limit", 0.9876) == pytest.approx(0.987) assert derive_corsika_limits._round_value("lower_energy_limit", 2.0) == pytest.approx(2.0) - # Test upper_radius_limit rounding - assert derive_corsika_limits._round_value("upper_radius_limit", 123.4) == 125 - assert derive_corsika_limits._round_value("upper_radius_limit", 100.0) == 100 - assert derive_corsika_limits._round_value("upper_radius_limit", 101.0) == 125 - assert derive_corsika_limits._round_value("upper_radius_limit", 75.0) == 75 + # Test upper_radius_limit_ground and upper_radius_limit_shower rounding + assert derive_corsika_limits._round_value("upper_radius_limit_ground", 123.4) == 125 + assert derive_corsika_limits._round_value("upper_radius_limit_ground", 100.0) == 100 + assert derive_corsika_limits._round_value("upper_radius_limit_ground", 101.0) == 125 + assert derive_corsika_limits._round_value("upper_radius_limit_ground", 75.0) == 75 + assert derive_corsika_limits._round_value("upper_radius_limit_shower", 123.4) == 125 + assert derive_corsika_limits._round_value("upper_radius_limit_shower", 100.0) == 100 # Test viewcone_radius rounding assert derive_corsika_limits._round_value("viewcone_radius", 1.1) == pytest.approx(1.25) @@ -361,14 +371,18 @@ def test_compute_upper_radius_limit(mocker): result = derive_corsika_limits.compute_upper_radius_limit(mock_histograms, 0.2) - assert isinstance(result, u.Quantity) - assert result.unit == u.m - assert result.value > 0 + assert isinstance(result, dict) + assert "upper_radius_limit_ground" in result + assert "upper_radius_limit_shower" in result + assert result["upper_radius_limit_ground"].unit == u.m + assert result["upper_radius_limit_shower"].unit == u.m + assert result["upper_radius_limit_ground"].value > 0 + assert result["upper_radius_limit_shower"].value > 0 expected = ( derive_corsika_limits._compute_limits(mock_hist, mock_bins, 0.2, limit_type="upper") * u.m ) - assert result == expected + assert result["upper_radius_limit_ground"] == expected def test_is_close(caplog): @@ -404,7 +418,10 @@ def test_process_file_with_mocked_histograms(mocker): ) mock_compute_upper_radius_limit = mocker.patch( COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, + return_value={ + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, + }, ) mock_compute_viewcone = mocker.patch( COMPUTE_VIEWCONE_PATH, @@ -425,7 +442,8 @@ def test_process_file_with_mocked_histograms(mocker): "azimuth": None, "nsb_level": None, "lower_energy_limit": 1.0 * u.TeV, - "upper_radius_limit": 100.0 * u.m, + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, "viewcone_radius": 2.0 * u.deg, } @@ -456,12 +474,19 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): COMPUTE_LOWER_ENERGY_LIMIT_PATH, return_value=1.0 * u.TeV, ) - mock_compute_upper_radius_limit = mocker.patch(COMPUTE_UPPER_RADIUS_LIMIT_PATH) + mock_compute_upper_radius_limit = mocker.patch( + COMPUTE_UPPER_RADIUS_LIMIT_PATH, + return_value={ + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, + }, + ) mock_compute_viewcone = mocker.patch(COMPUTE_VIEWCONE_PATH) mock_differential = mocker.patch( "simtools.production_configuration.derive_corsika_limits._compute_differential_limits", return_value={ - "upper_radius_limit": 120.0 * u.m, + "upper_radius_limit_ground": 120.0 * u.m, + "upper_radius_limit_shower": 114.0 * u.m, "viewcone_radius": 3.0 * u.deg, "core_vs_energy_curve": {"x": [100.0, 120.0], "y": [0.1, 1.0]}, "angular_distance_vs_energy_curve": {"x": [2.5, 3.0], "y": [0.1, 1.0]}, @@ -478,7 +503,8 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): ) assert result["lower_energy_limit"].value == pytest.approx(1.0) - assert result["upper_radius_limit"].value == pytest.approx(120.0) + assert result["upper_radius_limit_ground"].value == pytest.approx(120.0) + assert result["upper_radius_limit_shower"].value == pytest.approx(114.0) assert result["viewcone_radius"].value == pytest.approx(3.0) assert result["core_vs_energy_curve"] == {"x": [100.0, 120.0], "y": [0.1, 1.0]} assert result["angular_distance_vs_energy_curve"] == {"x": [2.5, 3.0], "y": [0.1, 1.0]} @@ -556,7 +582,8 @@ def test_process_file_passes_energy_bins_per_decade_to_histograms(mocker): mocker.patch( "simtools.production_configuration.derive_corsika_limits._compute_differential_limits", return_value={ - "upper_radius_limit": 120.0 * u.m, + "upper_radius_limit_ground": 120.0 * u.m, + "upper_radius_limit_shower": 114.0 * u.m, "viewcone_radius": 3.0 * u.deg, "core_vs_energy_curve": {"x": [100.0], "y": [1.0]}, "angular_distance_vs_energy_curve": {"x": [3.0], "y": [1.0]}, @@ -654,7 +681,10 @@ def test_process_file_with_plot_histograms(mocker, tmp_test_directory): ) mocker.patch( COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, + return_value={ + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, + }, ) mocker.patch( COMPUTE_VIEWCONE_PATH, @@ -684,7 +714,8 @@ def test_process_file_with_plot_histograms(mocker, tmp_test_directory): "azimuth": None, "nsb_level": None, "lower_energy_limit": 1.0 * u.TeV, - "upper_radius_limit": 100.0 * u.m, + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, "viewcone_radius": 2.0 * u.deg, } assert kwargs["array_name"] == "MockArray" @@ -862,7 +893,10 @@ def test_execute_production_job_single_job(mocker): ) mocker.patch( COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, + return_value={ + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, + }, ) mocker.patch( COMPUTE_VIEWCONE_PATH, @@ -886,7 +920,8 @@ def test_execute_production_job_single_job(mocker): assert result["event_data_file"] == "pattern_*.hdf5" assert result["array_name"] == "LST" assert "lower_energy_limit" in result - assert "upper_radius_limit" in result + assert "upper_radius_limit_ground" in result + assert "upper_radius_limit_shower" in result assert "viewcone_radius" in result @@ -906,7 +941,8 @@ def test_generate_corsika_limits_grid_multi_production(mocker, tmp_test_director production_index=1, event_data_file="pattern_2_*.hdf5", lower_energy_limit=0.6 * u.TeV, - upper_radius_limit=450.0 * u.m, + upper_radius_limit_ground=450.0 * u.m, + upper_radius_limit_shower=428.0 * u.m, viewcone_radius=5.5 * u.deg, ), ] @@ -1050,7 +1086,10 @@ def test_process_file_with_output_subdir(mocker, tmp_test_directory): ) mocker.patch( COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, + return_value={ + "upper_radius_limit_ground": 100.0 * u.m, + "upper_radius_limit_shower": 95.0 * u.m, + }, ) mocker.patch( COMPUTE_VIEWCONE_PATH, @@ -1112,7 +1151,8 @@ def mock_results(): "azimuth": 180.0 * u.deg, "nsb_level": 1.0, "lower_energy_limit": 0.5 * u.TeV, - "upper_radius_limit": 400.0 * u.m, + "upper_radius_limit_ground": 400.0 * u.m, + "upper_radius_limit_shower": 380.0 * u.m, "viewcone_radius": 5.0 * u.deg, } ] From 068eb77c2230068d5e96b40fc4c84d63e5ab7f95 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Mon, 18 May 2026 15:43:30 +0200 Subject: [PATCH 6/7] tests --- .../test_derive_corsika_limits.py | 18 +++++++++--------- .../test_plot_simtel_event_histograms.py | 12 ++++++------ 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py index 9933478f1a..cfc5a8deae 100644 --- a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py @@ -75,9 +75,9 @@ def test_process_file_passes_event_data_patterns_through(mocker): mock_histogram_class.assert_called_once_with( "input/*.h5", - array_name="array_name", - telescope_list=[1, 2], - energy_bins_per_decade=10, + "array_name", + [1, 2], + 10, ) @@ -449,9 +449,9 @@ def test_process_file_with_mocked_histograms(mocker): mock_histogram_class.assert_called_once_with( MOCK_FILE_PATH, - array_name="MockArray", - telescope_list=[1, 2], - energy_bins_per_decade=10, + "MockArray", + [1, 2], + 10, ) mock_histograms.fill.assert_called_once() mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2, 10) @@ -601,9 +601,9 @@ def test_process_file_passes_energy_bins_per_decade_to_histograms(mocker): mock_event_histograms.assert_called_once_with( MOCK_FILE_PATH, - array_name="MockArray", - telescope_list=[1, 2], - energy_bins_per_decade=6, + "MockArray", + [1, 2], + 6, ) diff --git a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py index 8b75b2267e..8fdf6ba6f9 100644 --- a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py +++ b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py @@ -437,7 +437,7 @@ def test_create_2d_plot_config(): "event_type": TRIGGERED, } limits = { - "upper_radius_limit": MagicMock(value=100), + "upper_radius_limit_ground": MagicMock(value=100), "lower_energy_limit": MagicMock(value=0.1), "viewcone_radius": MagicMock(value=5), } @@ -488,7 +488,7 @@ def test_create_2d_plot_config_core_xy(): "event_type": TRIGGERED, } limits = { - "upper_radius_limit": MagicMock(value=100), + "upper_radius_limit_ground": MagicMock(value=100), "lower_energy_limit": MagicMock(value=0.1), "viewcone_radius": MagicMock(value=5), } @@ -583,7 +583,7 @@ def test_get_limits(): # Test with limits containing all required keys limits = { "lower_energy_limit": MagicMock(value=42), - "upper_radius_limit": MagicMock(value=100), + "upper_radius_limit_ground": MagicMock(value=100), "viewcone_radius": MagicMock(value=5), } result = _get_limits("energy", limits) @@ -592,7 +592,7 @@ def test_get_limits(): # Test with partial limits (should not raise, but will return x only) limits = { "lower_energy_limit": MagicMock(value=42), - "upper_radius_limit": MagicMock(value=100), + "upper_radius_limit_ground": MagicMock(value=100), "viewcone_radius": MagicMock(value=5), } result = _get_limits("core_distance", limits) @@ -600,7 +600,7 @@ def test_get_limits(): # Test with all limits provided limits = { - "upper_radius_limit": MagicMock(value=100), + "upper_radius_limit_ground": MagicMock(value=100), "lower_energy_limit": MagicMock(value=0.1), "viewcone_radius": MagicMock(value=5), } @@ -640,7 +640,7 @@ def mock_histograms(): def test_plot_with_output_path(mock_histograms): output_path = Path("/mock/output/path") - limits = {"upper_radius_limit": MagicMock(value=100)} + limits = {"upper_radius_limit_ground": MagicMock(value=100)} array_name = "test_array" with ( From 0cc5a4cde6ff35db17002f46eca5f0f59975ead0 Mon Sep 17 00:00:00 2001 From: Gernot Maier Date: Wed, 20 May 2026 18:40:47 +0200 Subject: [PATCH 7/7] fix ground-to-shower transformation [skip ci] --- .../production_derive_corsika_limits.py | 57 ++-- .../derive_corsika_limits.py | 274 +++++++----------- src/simtools/sim_events/histograms.py | 11 +- src/simtools/sim_events/reader.py | 4 +- src/simtools/utils/geometry.py | 24 +- .../plot_simtel_event_histograms.py | 21 +- .../plot_simulated_event_distributions.yml | 2 +- .../test_derive_corsika_limits.py | 77 ++++- .../unit_tests/sim_events/test_histograms.py | 12 +- .../test_plot_simtel_event_histograms.py | 22 +- 10 files changed, 254 insertions(+), 250 deletions(-) diff --git a/src/simtools/applications/production_derive_corsika_limits.py b/src/simtools/applications/production_derive_corsika_limits.py index 16a6eca32e..4d40f790d5 100644 --- a/src/simtools/applications/production_derive_corsika_limits.py +++ b/src/simtools/applications/production_derive_corsika_limits.py @@ -15,29 +15,40 @@ Limits are computed from configurable per-axis allowed-loss settings. Results are provided as a table with the following columns: -+---------------------+-----------+--------+-----------------------------------------------+ -| Field | Data Type | Units | Description | -+=====================+===========+========+===============================================+ -| primary_particle | string | | Particle type (e.g., gamma, proton). | -+---------------------+-----------+--------+-----------------------------------------------+ -| array_name | string | | Array name (custom or as defined in | -| | | | 'array_layouts'). | -+---------------------+-----------+--------+-----------------------------------------------+ -| telescope_ids | string | | Comma-separated list of telescope IDs | -| | | | of this array. | -+---------------------+-----------+--------+-----------------------------------------------+ -| zenith | float64 | deg | Direction of array pointing zenith. | -+---------------------+-----------+--------+-----------------------------------------------+ -| azimuth | float64 | deg | Direction of array pointing azimuth. | -+---------------------+-----------+--------+-----------------------------------------------+ -| nsb_level | float64 | | Night sky background level. | -+---------------------+-----------+--------+-----------------------------------------------+ -| lower_energy_limit | float64 | TeV | Derived lower energy limit (**ERANGE**) | -+---------------------+-----------+--------+-----------------------------------------------+ -| upper_radius_limit | float64 | m | Derived upper core distance limit (**CSCAT**) | -+---------------------+-----------+--------+-----------------------------------------------+ -| viewcone_radius | float64 | deg | Derived viewcone radius limit (**VIEWCONE**) | -+---------------------+-----------+--------+-----------------------------------------------+ ++---------------------------+-----------+--------+----------------------------------------------+ +| Field | Data Type | Units | Description | ++===========================+===========+========+==============================================+ +| production_index | int64 | | Production index for multi-production runs. | ++---------------------------+-----------+--------+----------------------------------------------+ +| event_data_file | string | | Input event-data pattern for this row. | ++---------------------------+-----------+--------+----------------------------------------------+ +| primary_particle | string | | Particle type (e.g., gamma, proton). | ++---------------------------+-----------+--------+----------------------------------------------+ +| array_name | string | | Array name (custom or as defined in | +| | | | array_layouts). | ++---------------------------+-----------+--------+----------------------------------------------+ +| zenith | float64 | deg | Direction of array pointing zenith. | ++---------------------------+-----------+--------+----------------------------------------------+ +| azimuth | float64 | deg | Direction of array pointing azimuth. | ++---------------------------+-----------+--------+----------------------------------------------+ +| nsb_level | float64 | | Night sky background level. | ++---------------------------+-----------+--------+----------------------------------------------+ +| lower_energy_limit | float64 | TeV | Derived lower energy limit (**ERANGE**). | ++---------------------------+-----------+--------+----------------------------------------------+ +| upper_radius_limit | float64 | m | Derived upper core distance | +| | | | (**CSCAT**). | ++---------------------------+-----------+--------+----------------------------------------------+ +| viewcone_radius | float64 | deg | Derived viewcone radius (**VIEWCONE**). | ++---------------------------+-----------+--------+----------------------------------------------+ +| br_energy_min | float64 | TeV | Energy min from broad-range simulations. | ++---------------------------+-----------+--------+----------------------------------------------+ +| br_energy_max | float64 | TeV | Energy max from broad-range simulations. | ++---------------------------+-----------+--------+----------------------------------------------+ +| br_core_scatter_max | float64 | m | Core scatter max from broad-range | +| | | | simulations. | ++---------------------------+-----------+--------+----------------------------------------------+ +| br_viewcone_max | float64 | deg | Viewcone max from broad-range simulations. | ++---------------------------+-----------+--------+----------------------------------------------+ The input event data files are generated using the application simtools-generate-simtel-event-data and are required for each point in the observational parameter space (e.g., array pointing diff --git a/src/simtools/production_configuration/derive_corsika_limits.py b/src/simtools/production_configuration/derive_corsika_limits.py index 7657c3ef9f..2c6a1885b9 100644 --- a/src/simtools/production_configuration/derive_corsika_limits.py +++ b/src/simtools/production_configuration/derive_corsika_limits.py @@ -21,20 +21,34 @@ _logger = logging.getLogger(__name__) FILE_INFO_KEYS = ("primary_particle", "zenith", "azimuth", "nsb_level") +BROAD_RANGE_FILE_INFO_KEYS = { + "br_energy_min": "energy_min", + "br_energy_max": "energy_max", + "br_core_scatter_max": "core_scatter_max", + "br_viewcone_max": "viewcone_max", +} +COLUMN_DESCRIPTIONS = { + "br_energy_min": "Energy min from broad-range simulations.", + "br_energy_max": "Energy max from broad-range simulations.", + "br_core_scatter_max": "Core scatter max from broad-range simulations.", + "br_viewcone_max": "Viewcone max from broad-range simulations.", +} LOSS_AXES = ("energy", "core_distance", "angular_distance") RESULT_COLUMNS = [ "production_index", "event_data_file", "primary_particle", "array_name", - "telescope_ids", "zenith", "azimuth", "nsb_level", "lower_energy_limit", - "upper_radius_limit_ground", - "upper_radius_limit_shower", + "upper_radius_limit", "viewcone_radius", + "br_energy_min", + "br_energy_max", + "br_core_scatter_max", + "br_viewcone_max", ] @@ -373,28 +387,20 @@ def _process_file( allowed_losses["energy"]["loss_min_events"], ), } - if differential_loss_bins_per_decade > 0: - limits.update( - _compute_differential_limits( - histograms, - allowed_losses, - differential_loss_bins_per_decade, - ) - ) - else: - radius_limits = compute_upper_radius_limit( - histograms, - allowed_losses["core_distance"]["loss_fraction"], - allowed_losses["core_distance"]["loss_min_events"], - ) - limits.update(radius_limits) - limits["viewcone_radius"] = compute_viewcone( + limits.update( + _compute_limits( histograms, - allowed_losses["angular_distance"]["loss_fraction"], - allowed_losses["angular_distance"]["loss_min_events"], + allowed_losses, + differential_loss_bins_per_decade, ) - + ) limits.update({key: histograms.file_info.get(key) for key in FILE_INFO_KEYS}) + limits.update( + { + output_key: histograms.file_info.get(file_info_key) + for output_key, file_info_key in BROAD_RANGE_FILE_INFO_KEYS.items() + } + ) if plot_histograms: plot_output_path = output_subdir or io_handler.IOHandler().get_output_directory() @@ -408,41 +414,64 @@ def _process_file( return limits -def _compute_differential_limits( - histograms, - allowed_losses, - bins_per_decade, -): - """Compute core and viewcone limits per energy bin and return max limits.""" - low = int(np.floor(np.log10(np.min(histograms.energy_bins)))) - high = int(np.ceil(np.log10(np.max(histograms.energy_bins)))) - diff_e_bins = np.logspace(low, high, (high - low) * bins_per_decade + 1) - - core_max, core_x, core_y = _differential_upper_limits( - histograms.histograms["core_vs_energy"]["histogram"], - histograms.core_distance_bins, - histograms.energy_bins, - diff_e_bins, - allowed_losses["core_distance"], - "core_scatter", - "m", - ) - vc_max, vc_x, vc_y = _differential_upper_limits( - histograms.histograms["angular_distance_vs_energy"]["histogram"], - histograms.view_cone_bins, - histograms.energy_bins, - diff_e_bins, - allowed_losses["angular_distance"], - "viewcone", - "deg", - ) +def _compute_limits(histograms, allowed_losses, bins_per_decade): + """ + Compute core and viewcone limits per energy bin and return max limits. - upper_radius_limit_ground = core_max * u.m - upper_radius_limit_shower = _core_distance_ground_to_shower( - upper_radius_limit_ground, histograms.file_info.get("zenith") - ) - upper_radius_limit_shower = _is_close( - upper_radius_limit_shower, + Apply the allowed loss criteria per energy bin and return the maximum limit + """ + energy_range = [float(histograms.energy_bins[0]), float(histograms.energy_bins[-1])] + axis_configs = { + "core_distance": { + "x_bins": histograms.core_distance_bins, + "name": "core_scatter", + "unit": "m", + }, + "angular_distance": { + "x_bins": histograms.view_cone_bins, + "name": "viewcone", + "unit": "deg", + }, + } + + per_axis_limits = {} + + for axis_name, config in axis_configs.items(): + if bins_per_decade > 0: + low = int(np.floor(np.log10(np.min(histograms.energy_bins)))) + high = int(np.ceil(np.log10(np.max(histograms.energy_bins)))) + axis_max, curve_x, curve_y = _differential_upper_limits( + histograms.histograms[f"{axis_name}_vs_energy"]["histogram"], + config["x_bins"], + histograms.energy_bins, + np.logspace(low, high, (high - low) * bins_per_decade + 1), + allowed_losses[axis_name], + config["name"], + config["unit"], + ) + else: + axis_max = _integral_limits( + histograms.histograms[axis_name]["histogram"], + config["x_bins"], + allowed_losses[axis_name]["loss_fraction"], + allowed_losses[axis_name]["loss_min_events"], + limit_type="upper", + ) + curve_x = [axis_max, axis_max] + curve_y = energy_range + + per_axis_limits[axis_name] = { + "max": axis_max, + "curve": {"x": curve_x, "y": curve_y}, + "curve_key": f"{axis_name}_vs_energy_curve", + } + + core_max = per_axis_limits["core_distance"]["max"] + vc_max = per_axis_limits["angular_distance"]["max"] + + upper_radius_limit = core_max * u.m + upper_radius_limit = _is_close( + upper_radius_limit, histograms.file_info["core_scatter_max"].to("m") if "core_scatter_max" in histograms.file_info else None, @@ -455,19 +484,15 @@ def _compute_differential_limits( else None, "Upper viewcone limit is equal to the maximum viewcone distance of", ) - _logger.info( - f"Differential upper_radius_limit_ground (max over bins): {upper_radius_limit_ground}" - ) - _logger.info( - f"Differential upper_radius_limit_shower (max over bins): {upper_radius_limit_shower}" - ) + _logger.info(f"Differential upper_radius_limit (max over bins): {upper_radius_limit}") _logger.info(f"Differential viewcone_radius (max over bins): {viewcone_radius}") return { - "upper_radius_limit_ground": upper_radius_limit_ground, - "upper_radius_limit_shower": upper_radius_limit_shower, + "upper_radius_limit": upper_radius_limit, "viewcone_radius": viewcone_radius, - "core_vs_energy_curve": {"x": core_x, "y": core_y}, - "angular_distance_vs_energy_curve": {"x": vc_x, "y": vc_y}, + per_axis_limits["core_distance"]["curve_key"]: per_axis_limits["core_distance"]["curve"], + per_axis_limits["angular_distance"]["curve_key"]: per_axis_limits["angular_distance"][ + "curve" + ], } @@ -491,7 +516,7 @@ def _differential_upper_limits( total = float(np.sum(projected)) if total <= 0: continue - limit = _compute_limits( + limit = _integral_limits( projected, x_bins, allowed_loss["loss_fraction"], @@ -601,7 +626,7 @@ def _round_value(key, val): """Round value based on key type.""" if key == "lower_energy_limit": return np.floor(val * 1e3) / 1e3 - if key in ("upper_radius_limit_ground", "upper_radius_limit_shower"): + if key == "upper_radius_limit": return np.ceil(val / 25) * 25 if key == "viewcone_radius": return np.ceil(val / 0.25) * 0.25 @@ -613,17 +638,24 @@ def _create_table_columns(cols, columns, units): table_cols = [] for k in cols: col_data = columns[k] + col_description = COLUMN_DESCRIPTIONS.get(k) if any(isinstance(v, list | tuple) for v in col_data): - col = Column(data=col_data, name=k, unit=units.get(k), dtype=object) + col = Column( + data=col_data, + name=k, + unit=units.get(k), + dtype=object, + description=col_description, + ) else: - col = Column(data=col_data, name=k, unit=units.get(k)) + col = Column(data=col_data, name=k, unit=units.get(k), description=col_description) table_cols.append(col) return table_cols -def _compute_limits(hist, bin_edges, loss_fraction, loss_min_events=10, limit_type="lower"): +def _integral_limits(hist, bin_edges, loss_fraction, loss_min_events=10, limit_type="lower"): """ - Compute the limits based on the loss fraction and minimal required lost events. + Compute integral limits based on the loss fraction and minimal required lost events. Add or subtract one bin to be on the safe side of the limit. @@ -683,7 +715,7 @@ def compute_lower_energy_limit(histograms, loss_fraction, loss_min_events=10): Lower energy limit. """ energy_min = ( - _compute_limits( + _integral_limits( histograms.histograms["energy"]["histogram"], histograms.energy_bins, loss_fraction, @@ -709,96 +741,6 @@ def _is_close(value, reference, warning_text): return value -def _core_distance_ground_to_shower(core_distance, zenith): - """Convert core distance from ground to shower coordinates using zenith.""" - if zenith is None: - return core_distance - return core_distance * np.cos(zenith.to("rad").value) - - -def compute_upper_radius_limit(histograms, loss_fraction, loss_min_events=10): - """ - Compute the upper radial distance in both ground and shower coordinates. - - Parameters - ---------- - histograms : EventDataHistograms - Histograms. - loss_fraction : float - Fraction of events to be lost. - loss_min_events : int, optional - Minimum number of events to be lost after applying a derived limit. - - Returns - ------- - dict - Dictionary containing: - - "upper_radius_limit_ground": Upper radial distance in ground coordinates (m). - - "upper_radius_limit_shower": Upper radial distance in shower coordinates (m). - """ - radius_limit_ground = ( - _compute_limits( - histograms.histograms["core_distance"]["histogram"], - histograms.core_distance_bins, - loss_fraction, - loss_min_events, - limit_type="upper", - ) - * u.m - ) - radius_limit_shower = _core_distance_ground_to_shower( - radius_limit_ground, - histograms.file_info.get("zenith"), - ) - radius_limit_shower = _is_close( - radius_limit_shower, - histograms.file_info["core_scatter_max"].to("m") - if "core_scatter_max" in histograms.file_info - else None, - "Upper radius limit is equal to the maximum core scatter distance of", - ) - return { - "upper_radius_limit_ground": radius_limit_ground, - "upper_radius_limit_shower": radius_limit_shower, - } - - -def compute_viewcone(histograms, loss_fraction, loss_min_events=10): - """ - Compute the viewcone based on the event loss fraction. - - The shower IDs of triggered events are used to create a mask for the - azimuth and altitude of the triggered events. A mapping is created - between the triggered events and the simulated events using the shower IDs. - - Parameters - ---------- - histograms : EventDataHistograms - Histograms. - loss_fraction : float - Fraction of events to be lost. - loss_min_events : int, optional - Minimum number of events to be lost after applying a derived limit. - - Returns - ------- - astropy.units.Quantity - Viewcone radius in degrees. - """ - viewcone_limit = ( - _compute_limits( - histograms.histograms["angular_distance"]["histogram"], - histograms.view_cone_bins, - loss_fraction, - loss_min_events, - limit_type="upper", - ) - * u.deg - ) - return _is_close( - viewcone_limit, - histograms.file_info["viewcone_max"].to("deg") - if "viewcone_max" in histograms.file_info - else None, - "Upper viewcone limit is equal to the maximum viewcone distance of", - ) +def _core_distance_ground_to_shower(core_distance, _zenith): + """Return unchanged core distance; limits are handled in shower coordinates.""" + return core_distance diff --git a/src/simtools/sim_events/histograms.py b/src/simtools/sim_events/histograms.py index 4d64340bf1..0ceac42ddc 100644 --- a/src/simtools/sim_events/histograms.py +++ b/src/simtools/sim_events/histograms.py @@ -216,7 +216,7 @@ def _define_histograms(self, event_data, triggered_data, shower_data): "is_1d": False, "axis_titles": ["Core X (m)", "Core Y (m)", event_count_axis_title], }, - "core_vs_energy": { + "core_distance_vs_energy": { "event_data_column": ("core_distance_shower", "simulated_energy"), "event_data": (event_data, event_data), "bin_edges": (self.core_distance_bins, self.energy_bins), @@ -393,19 +393,14 @@ def core_distance_bins(self): """ Return bins for the core distance histogram. - CORSIKA CSCAT is defined in the shower plane, shower coordinates - are in ground coordinates. The core distance bins in ground coordinates - are therefore scaled with 1/cos(zenith). + CORSIKA CSCAT ('core_scatter_max') is defined in the shower plane. """ - zenith = self.file_info.get("zenith", 0.0 * u.deg).to("rad").value - scaling_factor = 1 / np.cos(zenith) - if "core_distance_bin_edges" in self.histograms: return self.histograms["core_distance_bin_edges"] return np.linspace( self.file_info.get("core_scatter_min", 0.0 * u.m).to("m").value, - self.file_info.get("core_scatter_max", 1.0e5 * u.m).to("m").value * scaling_factor, + self.file_info.get("core_scatter_max", 1.0e5 * u.m).to("m").value, 100, ) diff --git a/src/simtools/sim_events/reader.py b/src/simtools/sim_events/reader.py index b4f0ef9cc2..4cda2f30b8 100644 --- a/src/simtools/sim_events/reader.py +++ b/src/simtools/sim_events/reader.py @@ -157,8 +157,8 @@ def _table_to_shower_data(self, table): shower_data.x_core, shower_data.y_core, 0.0, - shower_data.shower_azimuth, - shower_data.shower_altitude, + np.deg2rad(shower_data.shower_azimuth), + np.deg2rad(shower_data.shower_altitude), ) ) shower_data.core_distance_shower = np.hypot( diff --git a/src/simtools/utils/geometry.py b/src/simtools/utils/geometry.py index ceea44c5a2..0b5e894ab6 100644 --- a/src/simtools/utils/geometry.py +++ b/src/simtools/utils/geometry.py @@ -121,11 +121,13 @@ def solid_angle(angle_max, angle_min=0 * u.rad): return 2 * np.pi * (np.cos(angle_min.to("rad")) - np.cos(angle_max.to("rad"))) * u.sr -def transform_ground_to_shower_coordinates(x_ground, y_ground, z_ground, azimuth, altitude): +def transform_ground_to_shower_coordinates( + x_ground, y_ground, z_ground, azimuth_from_north_east, altitude +): """ Transform ground to shower coordinates. - Assume ground to be of type 'North-West-Up' (NWU) coordinates. + Reverts the rotation applied in sim_telarray with the same convention as in iact.c. Parameters ---------- @@ -145,14 +147,20 @@ def transform_ground_to_shower_coordinates(x_ground, y_ground, z_ground, azimuth tuple Transformed shower coordinates (x', y', z'). """ - x, y, z, az, alt = np.broadcast_arrays(x_ground, y_ground, z_ground, azimuth, altitude) + x, y, z, az, alt = np.broadcast_arrays( + x_ground, y_ground, z_ground, azimuth_from_north_east, altitude + ) + + # Convert to CORSIKA phi convention used in iact.c: + # 0 = moves North, 90 deg = moves West + phi = np.mod(np.pi - az, 2.0 * np.pi) - ca, sa = np.cos(az), np.sin(az) - cz, sz = np.sin(alt), np.cos(alt) + ca, sa = np.cos(phi), np.sin(phi) + cos_theta = np.sin(alt) # theta = pi/2 - altitude - x_s = ca * cz * x - sa * y + ca * sz * z - y_s = sa * cz * x + ca * y + sa * sz * z - z_s = -sz * x + cz * z + x_s = (ca * ca * cos_theta + sa * sa) * x + ca * sa * (cos_theta - 1.0) * y + y_s = ca * sa * (cos_theta - 1.0) * x + (sa * sa * cos_theta + ca * ca) * y + z_s = z return x_s, y_s, z_s diff --git a/src/simtools/visualization/plot_simtel_event_histograms.py b/src/simtools/visualization/plot_simtel_event_histograms.py index 74a477d214..9c9218d1d0 100644 --- a/src/simtools/visualization/plot_simtel_event_histograms.py +++ b/src/simtools/visualization/plot_simtel_event_histograms.py @@ -21,11 +21,9 @@ def plot(histograms, output_path=None, limits=None, array_name=None): Directory to save plots. If None, plots will be displayed. limits: dict, optional Dictionary containing limits for plotting. Keys can include: - - "upper_radius_limit_ground": Upper limit for core distance (ground coordinates) - - "upper_radius_limit_shower": Upper limit for core distance (shower coordinates) + - "upper_radius_limit": Upper limit for core distance - "lower_energy_limit": Lower limit for energy - "viewcone_radius": Radius for the viewcone - Note: Plotting always uses ground coordinates (upper_radius_limit_ground). array_name: str, optional Name of the telescope array configuration. """ @@ -40,7 +38,6 @@ def _get_limits(name, limits): Extract limits from the provided dictionary for plotting. Fine tuned to expected histograms to be plotted. - Always uses ground coordinates for upper_radius_limit. """ def _safe_value(limits, key): @@ -49,17 +46,17 @@ def _safe_value(limits, key): mapping = { "energy": {"x": _safe_value(limits, "lower_energy_limit")}, - "core_distance": {"x": _safe_value(limits, "upper_radius_limit_ground")}, + "core_distance": {"x": _safe_value(limits, "upper_radius_limit")}, "angular_distance": {"x": _safe_value(limits, "viewcone_radius")}, - "core_vs_energy": { - "x": _safe_value(limits, "upper_radius_limit_ground"), + "core_distance_vs_energy": { + "x": _safe_value(limits, "upper_radius_limit"), "y": _safe_value(limits, "lower_energy_limit"), - "curve": limits.get("core_vs_energy_curve"), + "curve": limits.get("core_distance_vs_energy_curve"), }, - "core_vs_energy_cumulative": { - "x": _safe_value(limits, "upper_radius_limit_ground"), + "core_distance_vs_energy_cumulative": { + "x": _safe_value(limits, "upper_radius_limit"), "y": _safe_value(limits, "lower_energy_limit"), - "curve": limits.get("core_vs_energy_curve"), + "curve": limits.get("core_distance_vs_energy_curve"), }, "angular_distance_vs_energy": { "x": _safe_value(limits, "viewcone_radius"), @@ -71,7 +68,7 @@ def _safe_value(limits, key): "y": _safe_value(limits, "lower_energy_limit"), "curve": limits.get("angular_distance_vs_energy_curve"), }, - "x_core_shower_vs_y_core_shower": {"r": _safe_value(limits, "upper_radius_limit_ground")}, + "x_core_shower_vs_y_core_shower": {"r": _safe_value(limits, "upper_radius_limit")}, } return mapping.get(name) diff --git a/tests/integration_tests/config/plot_simulated_event_distributions.yml b/tests/integration_tests/config/plot_simulated_event_distributions.yml index eed7bbea47..71ede0d103 100644 --- a/tests/integration_tests/config/plot_simulated_event_distributions.yml +++ b/tests/integration_tests/config/plot_simulated_event_distributions.yml @@ -5,7 +5,7 @@ applications: event_data_file: tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_reduced_event_data.hdf5 output_path: simtools-output integration_tests: - - output_file: core_vs_energy.png + - output_file: core_distance_vs_energy.png - output_file: energy_mc.png test_name: sim_telarray_input schema_name: application_workflow.metaschema diff --git a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py index cfc5a8deae..235fe0c825 100644 --- a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py @@ -103,8 +103,18 @@ def test_create_results_table(mock_results): assert isinstance(table, Table) assert len(table) == 1 + assert "telescope_ids" not in table.colnames assert "zenith" in table.colnames assert table["zenith"].unit == u.deg + assert "br_energy_min" in table.colnames + assert "br_energy_max" in table.colnames + assert "br_core_scatter_max" in table.colnames + assert "br_viewcone_max" in table.colnames + assert table["br_energy_min"].unit == u.TeV + assert table["br_energy_max"].unit == u.TeV + assert table["br_core_scatter_max"].unit == u.m + assert table["br_viewcone_max"].unit == u.deg + assert table["br_viewcone_max"].description == "Viewcone max from broad-range simulations." assert table.meta["loss_fraction_energy"] == pytest.approx(0.2) assert table.meta["loss_min_events_energy"] == 10 assert table.meta["loss_fraction_core_distance"] == pytest.approx(0.2) @@ -483,12 +493,12 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): ) mock_compute_viewcone = mocker.patch(COMPUTE_VIEWCONE_PATH) mock_differential = mocker.patch( - "simtools.production_configuration.derive_corsika_limits._compute_differential_limits", + "simtools.production_configuration.derive_corsika_limits._compute_limits", return_value={ "upper_radius_limit_ground": 120.0 * u.m, "upper_radius_limit_shower": 114.0 * u.m, "viewcone_radius": 3.0 * u.deg, - "core_vs_energy_curve": {"x": [100.0, 120.0], "y": [0.1, 1.0]}, + "core_distance_vs_energy_curve": {"x": [100.0, 120.0], "y": [0.1, 1.0]}, "angular_distance_vs_energy_curve": {"x": [2.5, 3.0], "y": [0.1, 1.0]}, }, ) @@ -506,7 +516,7 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): assert result["upper_radius_limit_ground"].value == pytest.approx(120.0) assert result["upper_radius_limit_shower"].value == pytest.approx(114.0) assert result["viewcone_radius"].value == pytest.approx(3.0) - assert result["core_vs_energy_curve"] == {"x": [100.0, 120.0], "y": [0.1, 1.0]} + assert result["core_distance_vs_energy_curve"] == {"x": [100.0, 120.0], "y": [0.1, 1.0]} assert result["angular_distance_vs_energy_curve"] == {"x": [2.5, 3.0], "y": [0.1, 1.0]} mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2, 10) @@ -526,16 +536,14 @@ def test_process_file_with_differential_loss_per_energy_bin(mocker): ({}, None, None), ], ) -def test_compute_differential_limits( - mocker, file_info, expected_core_scatter_max, expected_viewcone_max -): - """Test _compute_differential_limits forwards slices and preserves units.""" +def test_compute_limits(mocker, file_info, expected_core_scatter_max, expected_viewcone_max): + """Test _compute_limits forwards slices and preserves units.""" histograms = mocker.MagicMock() histograms.energy_bins = np.array([1.0, 10.0]) histograms.core_distance_bins = np.array([0.0, 100.0]) histograms.view_cone_bins = np.array([0.0, 5.0]) histograms.histograms = { - "core_vs_energy": {"histogram": "core-hist"}, + "core_distance_vs_energy": {"histogram": "core-hist"}, "angular_distance_vs_energy": {"histogram": "viewcone-hist"}, } histograms.file_info = file_info @@ -552,7 +560,7 @@ def test_compute_differential_limits( side_effect=[125.0 * u.m, 3.25 * u.deg], ) - derive_corsika_limits._compute_differential_limits(histograms, DEFAULT_ALLOWED_LOSSES, 2) + derive_corsika_limits._compute_limits(histograms, DEFAULT_ALLOWED_LOSSES, 2) expected_diff_bins = np.logspace(0, 1, 3) np.testing.assert_allclose(mock_diff_limits.call_args_list[0].args[3], expected_diff_bins) @@ -570,6 +578,42 @@ def test_compute_differential_limits( assert mock_is_close.call_args_list[1].args[1] == expected_viewcone_max +def test_compute_limits_with_integral_fallback_curves(mocker): + """Test _compute_limits returns energy-independent curves for integral limits.""" + histograms = mocker.MagicMock() + histograms.energy_bins = np.array([1.0, 10.0]) + histograms.core_distance_bins = np.array([0.0, 100.0]) + histograms.view_cone_bins = np.array([0.0, 5.0]) + histograms.histograms = { + "core_distance": {"histogram": np.array([1.0, 2.0])}, + "angular_distance": {"histogram": np.array([3.0, 4.0])}, + } + histograms.file_info = {} + + mock_integral_limits = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._integral_limits", + side_effect=[120.0, 3.0], + ) + mock_diff_limits = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._differential_upper_limits" + ) + mocker.patch( + "simtools.production_configuration.derive_corsika_limits._is_close", + side_effect=lambda value, *_: value, + ) + + result = derive_corsika_limits._compute_limits( + histograms, + DEFAULT_ALLOWED_LOSSES, + bins_per_decade=0, + ) + + assert mock_integral_limits.call_count == 2 + mock_diff_limits.assert_not_called() + assert result["core_distance_vs_energy_curve"] == {"x": [120.0, 120.0], "y": [1.0, 10.0]} + assert result["angular_distance_vs_energy_curve"] == {"x": [3.0, 3.0], "y": [1.0, 10.0]} + + def test_process_file_passes_energy_bins_per_decade_to_histograms(mocker): """Test differential binning resolution is forwarded to EventDataHistograms.""" mock_histograms = mocker.MagicMock() @@ -580,12 +624,12 @@ def test_process_file_passes_energy_bins_per_decade_to_histograms(mocker): ) mocker.patch(COMPUTE_LOWER_ENERGY_LIMIT_PATH, return_value=1.0 * u.TeV) mocker.patch( - "simtools.production_configuration.derive_corsika_limits._compute_differential_limits", + "simtools.production_configuration.derive_corsika_limits._compute_limits", return_value={ "upper_radius_limit_ground": 120.0 * u.m, "upper_radius_limit_shower": 114.0 * u.m, "viewcone_radius": 3.0 * u.deg, - "core_vs_energy_curve": {"x": [100.0], "y": [1.0]}, + "core_distance_vs_energy_curve": {"x": [100.0], "y": [1.0]}, "angular_distance_vs_energy_curve": {"x": [3.0], "y": [1.0]}, }, ) @@ -862,11 +906,11 @@ def test_build_production_subdirectories_creates_dirs(tmp_test_directory): assert output_subdir.isdir() -def test_core_distance_ground_to_shower_converts_with_zenith(): - """Test _core_distance_ground_to_shower scales by cos(zenith).""" +def test_core_distance_ground_to_shower_keeps_value_with_zenith(): + """Test _core_distance_ground_to_shower keeps value unchanged.""" result = derive_corsika_limits._core_distance_ground_to_shower(100.0 * u.m, 60.0 * u.deg) assert result.unit == u.m - assert result.value == pytest.approx(50.0) + assert result.value == pytest.approx(100.0) def test_core_distance_ground_to_shower_no_zenith_returns_input(): @@ -1146,7 +1190,6 @@ def mock_results(): { "primary_particle": "gamma", "array_name": "LST", - "telescope_ids": ["LSTN-01"], "zenith": 20.0 * u.deg, "azimuth": 180.0 * u.deg, "nsb_level": 1.0, @@ -1154,5 +1197,9 @@ def mock_results(): "upper_radius_limit_ground": 400.0 * u.m, "upper_radius_limit_shower": 380.0 * u.m, "viewcone_radius": 5.0 * u.deg, + "br_energy_min": 0.03 * u.TeV, + "br_energy_max": 300.0 * u.TeV, + "br_core_scatter_max": 800.0 * u.m, + "br_viewcone_max": 10.0 * u.deg, } ] diff --git a/tests/unit_tests/sim_events/test_histograms.py b/tests/unit_tests/sim_events/test_histograms.py index 842450ae49..bf26646979 100644 --- a/tests/unit_tests/sim_events/test_histograms.py +++ b/tests/unit_tests/sim_events/test_histograms.py @@ -703,7 +703,10 @@ def test_calculate_cumulative_data(mock_reader, hdf5_file_name): "energy": {"histogram": np.array([10, 20, 30, 40]), "axis_titles": ["E", ""]}, "core_distance": {"histogram": np.array([5, 15, 25, 35]), "axis_titles": ["r", ""]}, "angular_distance": {"histogram": np.array([2, 4, 6, 8]), "axis_titles": ["theta", ""]}, - "core_vs_energy": {"histogram": np.array([[1, 2], [3, 4]]), "axis_titles": ["r", "E"]}, + "core_distance_vs_energy": { + "histogram": np.array([[1, 2], [3, 4]]), + "axis_titles": ["r", "E"], + }, "angular_distance_vs_energy": { "histogram": np.array([[2, 3], [4, 5]]), "axis_titles": ["theta", "E"], @@ -715,13 +718,13 @@ def test_calculate_cumulative_data(mock_reader, hdf5_file_name): expected_cumulative_core_distance = np.array([5, 20, 45, 80]) expected_cumulative_angular_distance = np.array([2, 6, 12, 20]) # Expected normalized cumulative for 2D histograms along axis=0 (column-wise) - expected_norm_core_vs_energy = np.array([[1 / 4, 2 / 6], [1.0, 1.0]]) + expected_norm_core_distance_vs_energy = np.array([[1 / 4, 2 / 6], [1.0, 1.0]]) expected_norm_ang_vs_energy = np.array([[2 / 6, 3 / 8], [1.0, 1.0]]) assert set(cumulative_data.keys()) == { "core_distance_cumulative", "angular_distance_cumulative", "angular_distance_vs_energy_cumulative", - "core_vs_energy_cumulative", + "core_distance_vs_energy_cumulative", "energy_cumulative", } np.testing.assert_array_equal( @@ -735,7 +738,8 @@ def test_calculate_cumulative_data(mock_reader, hdf5_file_name): expected_cumulative_angular_distance, ) np.testing.assert_allclose( - cumulative_data["core_vs_energy_cumulative"]["histogram"], expected_norm_core_vs_energy + cumulative_data["core_distance_vs_energy_cumulative"]["histogram"], + expected_norm_core_distance_vs_energy, ) np.testing.assert_allclose( cumulative_data["angular_distance_vs_energy_cumulative"]["histogram"], diff --git a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py index 8fdf6ba6f9..458f96bd51 100644 --- a/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py +++ b/tests/unit_tests/visualization/test_plot_simtel_event_histograms.py @@ -423,11 +423,11 @@ def test_create_2d_plot_config(): "histogram": np.array([[1, 2], [3, 4]]), "bin_edges": [np.array([0, 1, 2]), np.array([0, 1, 2])], "plot_scales": {"y": "log"}, - "title": "Triggered events: core vs energy: core vs energy", + "title": "Triggered events: core distance vs energy", "axis_titles": [CORE_DIST_LABEL, ENERGY_LABEL, EVENT_COUNT], } config = { - "base_key": "core_vs_energy", + "base_key": "core_distance_vs_energy", "x_label": CORE_DIST_LABEL, "y_label": ENERGY_LABEL, "plot_params": {"norm": "log", "cmap": "viridis"}, @@ -443,7 +443,7 @@ def test_create_2d_plot_config(): } result = plot_simtel_event_histograms._create_2d_plot_config( histograms, - "core_vs_energy", + "core_distance_vs_energy", config, limits, ) @@ -456,16 +456,16 @@ def test_create_2d_plot_config(): assert result["labels"]["y"] == config["y_label"] assert ( result["labels"]["title"] - == "Triggered events: core vs energy: core vs energy: core vs energy" + == "Triggered events: core distance vs energy: core distance vs energy" ) # Accept lines from limits dict, not config assert ( result["labels"]["title"] - == "Triggered events: core vs energy: core vs energy: core vs energy" + == "Triggered events: core distance vs energy: core distance vs energy" ) assert result["colorbar_label"] in (config["colorbar_label"], None) - assert result["filename"] == "core_vs_energy" + assert result["filename"] == "core_distance_vs_energy" def test_create_2d_plot_config_core_xy(): @@ -607,18 +607,18 @@ def test_get_limits(): result = _get_limits("angular_distance", limits) assert result == {"x": 5} - limits["core_vs_energy_curve"] = {"x": [10, 20], "y": [0.1, 1.0]} + limits["core_distance_vs_energy_curve"] = {"x": [10, 20], "y": [0.1, 1.0]} limits["angular_distance_vs_energy_curve"] = {"x": [2.5, 3.0], "y": [0.1, 1.0]} - result = _get_limits("core_vs_energy", limits) + result = _get_limits("core_distance_vs_energy", limits) assert result["x"] == 100 assert result["y"] == pytest.approx(0.1) - assert result["curve"] == limits["core_vs_energy_curve"] + assert result["curve"] == limits["core_distance_vs_energy_curve"] - result = _get_limits("core_vs_energy_cumulative", limits) + result = _get_limits("core_distance_vs_energy_cumulative", limits) assert result["x"] == 100 assert result["y"] == pytest.approx(0.1) - assert result["curve"] == limits["core_vs_energy_curve"] + assert result["curve"] == limits["core_distance_vs_energy_curve"] result = _get_limits("angular_distance_vs_energy", limits) assert result["x"] == 5