diff --git a/docs/changes/2128.feature.md b/docs/changes/2128.feature.md new file mode 100644 index 0000000000..e54c79913f --- /dev/null +++ b/docs/changes/2128.feature.md @@ -0,0 +1 @@ +Allow custom list of telescopes for corsika limits derivation. diff --git a/src/simtools/applications/production_derive_corsika_limits.py b/src/simtools/applications/production_derive_corsika_limits.py index b90ff829c3..e8960f1885 100644 --- a/src/simtools/applications/production_derive_corsika_limits.py +++ b/src/simtools/applications/production_derive_corsika_limits.py @@ -80,6 +80,17 @@ --loss_fraction 1e-6 \\ --plot_histograms \\ --output_file corsika_simulation_limits.ecsv + +Derive limits for an inline list of array elements: + +.. code-block:: console + + simtools-production-derive-corsika-limits \\ + --event_data_file event_dat_file.hdf5 \\ + --array_element_list LSTN-01 LSTN-02 MSTN-03 \\ + --loss_fraction 1e-6 \\ + --plot_histograms \\ + --output_file corsika_simulation_limits.ecsv """ from simtools.application_control import build_application diff --git a/src/simtools/production_configuration/derive_corsika_limits.py b/src/simtools/production_configuration/derive_corsika_limits.py index d5f0c9eb70..df15722484 100644 --- a/src/simtools/production_configuration/derive_corsika_limits.py +++ b/src/simtools/production_configuration/derive_corsika_limits.py @@ -32,10 +32,17 @@ def generate_corsika_limits_grid(args_dict): args_dict.get("site"), args_dict.get("model_version"), ) - else: + elif args_dict.get("array_element_list"): + telescope_configs = {"array_element_list": args_dict["array_element_list"]} + elif args_dict.get("telescope_ids"): telescope_configs = ascii_handler.collect_data_from_file(args_dict["telescope_ids"])[ "telescope_configs" ] + else: + raise ValueError( + "No telescope configuration provided. Use one of --array_layout_name, " + "--array_element_list, or --telescope_ids." + ) results = [] for array_name, telescope_ids in telescope_configs.items(): diff --git a/src/simtools/sim_events/histograms.py b/src/simtools/sim_events/histograms.py index a4beb47d6b..f84e3e16b6 100644 --- a/src/simtools/sim_events/histograms.py +++ b/src/simtools/sim_events/histograms.py @@ -61,7 +61,13 @@ def fill(self): "scatter_area": _file_info_table["scatter_area"].to("cm2"), } - self.histograms = self._define_histograms(event_data, triggered_data, shower_data) + current_histograms = self._define_histograms(event_data, triggered_data, shower_data) + for name, hist in current_histograms.items(): + previous = self.histograms.get(name) + if previous is not None: + hist["histogram"] = previous["histogram"] + + self.histograms = current_histograms for data in self.histograms.values(): self._fill_histogram_and_bin_edges(data) diff --git a/src/simtools/sim_events/reader.py b/src/simtools/sim_events/reader.py index 640907b51f..78a7045d77 100644 --- a/src/simtools/sim_events/reader.py +++ b/src/simtools/sim_events/reader.py @@ -345,7 +345,13 @@ def get_reduced_simulation_file_info(self, simulation_file_info): dict Dictionary containing the reduced simulation file info. """ - particle_id = np.unique(simulation_file_info["particle_id"].data) + particle_id = np.unique( + self._convert_numeric_column( + simulation_file_info["particle_id"].data, + key="particle_id", + dtype=np.int64, + ) + ) keys = [ "zenith", "azimuth", @@ -359,10 +365,15 @@ def get_reduced_simulation_file_info(self, simulation_file_info): ] float_arrays = {} for key in keys: + column_values = self._convert_numeric_column( + simulation_file_info[key].data, + key=key, + dtype=np.float64, + ) if key == "energy_min": - float_arrays[key] = np.unique(np.round(simulation_file_info[key].data, decimals=3)) + float_arrays[key] = np.unique(np.round(column_values, decimals=3)) else: - float_arrays[key] = np.unique(np.round(simulation_file_info[key].data, decimals=2)) + float_arrays[key] = np.unique(np.round(column_values, decimals=2)) if any(len(arr) > 1 for arr in (particle_id, *(float_arrays[key] for key in keys))): self._logger.warning("Simulation file info has non-unique values.") @@ -391,6 +402,27 @@ def get_reduced_simulation_file_info(self, simulation_file_info): return reduced_info + def _convert_numeric_column(self, values, key, dtype): + """Convert table column values to a numeric numpy array.""" + array_values = np.asarray(values) + + if array_values.dtype.kind == "O": + array_values = np.array( + [ + value.decode("utf-8") if isinstance(value, bytes | bytearray) else value + for value in array_values + ] + ) + + try: + if array_values.dtype.kind in ("S", "U"): + return np.char.strip(array_values.astype("U")).astype(dtype) + return array_values.astype(dtype) + except (TypeError, ValueError) as exc: + raise ValueError( + f"Unable to convert FILE_INFO column '{key}' to numeric values." + ) from exc + def scatter_area(self, core_scatter_min, core_scatter_max): """ Calculate the scatter area of the core. diff --git a/src/simtools/sim_events/writer.py b/src/simtools/sim_events/writer.py index 27143e2ad7..5219e559a0 100644 --- a/src/simtools/sim_events/writer.py +++ b/src/simtools/sim_events/writer.py @@ -378,7 +378,15 @@ def get_nsb_level_from_sim_telarray_metadata(self, file): """ metadata, _ = read_sim_telarray_metadata(file) nsb_integrated_flux = metadata.get("nsb_integrated_flux") - return nsb_integrated_flux or self._get_nsb_level_from_file_name(str(file)) + if nsb_integrated_flux is not None: + try: + return float(nsb_integrated_flux) + except (TypeError, ValueError): + self._logger.warning( + f"Invalid nsb_integrated_flux value '{nsb_integrated_flux}' for {file}" + ) + + return self._get_nsb_level_from_file_name(str(file)) def _get_nsb_level_from_file_name(self, file): """ 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 597f936fc6..a82b3e72be 100644 --- a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py @@ -228,6 +228,41 @@ def test_generate_corsika_limits_grid_with_db_layouts(mocker, mock_args_dict): assert mock_write.call_count == 1 +def test_generate_corsika_limits_grid_with_array_element_list(mocker, mock_args_dict): + """Test generate_corsika_limits_grid using inline array_element_list.""" + args = mock_args_dict.copy() + args["array_element_list"] = ["LSTN-01", "LSTN-02", "MSTN-03"] + args["telescope_ids"] = None + + mock_collect = mocker.patch("simtools.io.ascii_handler.collect_data_from_file") + mock_process = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._process_file" + ) + mock_write = mocker.patch( + "simtools.production_configuration.derive_corsika_limits.write_results" + ) + + derive_corsika_limits.generate_corsika_limits_grid(args) + + mock_collect.assert_not_called() + mock_process.assert_called_once() + call_args = mock_process.call_args[0] + assert call_args[1] == "array_element_list" + assert call_args[2] == ["LSTN-01", "LSTN-02", "MSTN-03"] + assert mock_write.call_count == 1 + + +def test_generate_corsika_limits_grid_without_telescope_configuration(mock_args_dict): + """Test generate_corsika_limits_grid raises if no telescope input is provided.""" + args = mock_args_dict.copy() + args["array_layout_name"] = None + args["array_element_list"] = None + args["telescope_ids"] = None + + with pytest.raises(ValueError, match="No telescope configuration provided"): + derive_corsika_limits.generate_corsika_limits_grid(args) + + def test_compute_limits_lower(): hist = np.array([1, 2, 3, 4, 5]) bin_edges = np.array([0, 1, 2, 3, 4, 5]) diff --git a/tests/unit_tests/sim_events/test_histograms.py b/tests/unit_tests/sim_events/test_histograms.py index 018779822d..4fa53db5d0 100644 --- a/tests/unit_tests/sim_events/test_histograms.py +++ b/tests/unit_tests/sim_events/test_histograms.py @@ -206,6 +206,61 @@ def fake_fill_histogram_and_bin_edges(data): ) +def test_fill_accumulates_histograms_across_data_sets(mock_reader, hdf5_file_name, mocker): + """Test fill accumulates histogram counts across all indexed datasets.""" + histograms = EventDataHistograms(hdf5_file_name) + + mock_reader.return_value.data_sets = ["dataset_0", "dataset_1"] + mock_reader.return_value.read_event_data.return_value = ( + mocker.Mock(), + mocker.Mock(), + mocker.Mock(), + mocker.Mock(), + ) + mock_reader.return_value.get_reduced_simulation_file_info.return_value = { + "energy_min": 0.1 * u.TeV, + "core_scatter_max": 100.0 * u.m, + "viewcone_max": 2.0 * u.deg, + "solid_angle": 1.0 * u.sr, + "scatter_area": 1.0 * u.cm**2, + } + + dataset_0 = mocker.Mock() + dataset_0.simulated_energy = np.array([0.2, 0.4]) + dataset_1 = mocker.Mock() + dataset_1.simulated_energy = np.array([1.4]) + + mocker.patch.object( + histograms, + "_define_histograms", + side_effect=[ + { + "energy": { + "1d": True, + "event_data": dataset_0, + "event_data_column": "simulated_energy", + "bin_edges": np.array([0.0, 1.0, 2.0]), + "histogram": None, + } + }, + { + "energy": { + "1d": True, + "event_data": dataset_1, + "event_data_column": "simulated_energy", + "bin_edges": np.array([0.0, 1.0, 2.0]), + "histogram": None, + } + }, + ], + ) + + histograms.fill() + + np.testing.assert_array_equal(histograms.histograms["energy"]["histogram"], np.array([2, 1])) + assert mock_reader.return_value.read_event_data.call_count == 2 + + def test_calculate_cumulative_histogram(mock_reader, hdf5_file_name): """Test calculation of cumulative histogram.""" histograms = EventDataHistograms(hdf5_file_name) diff --git a/tests/unit_tests/sim_events/test_reader.py b/tests/unit_tests/sim_events/test_reader.py index fbd84e105f..405f19a196 100644 --- a/tests/unit_tests/sim_events/test_reader.py +++ b/tests/unit_tests/sim_events/test_reader.py @@ -2,6 +2,7 @@ from unittest.mock import patch import astropy.units as u +import numpy as np import pytest from astropy.table import Table from astropy.tests.helper import assert_quantity_allclose @@ -160,6 +161,54 @@ def test_get_reduced_simulation_info_with_warning(mock_primary_particle, mock_fi assert info["zenith"] == pytest.approx(20.0) # Should use first value +def test_get_reduced_simulation_info_with_string_encoded_numeric_values(mock_fits_file): + """Test conversion of numeric FILE_INFO columns stored as byte strings.""" + reader = EventDataReader(mock_fits_file) + + file_info = Table() + file_info.meta["EXTNAME"] = "FILE_INFO" + file_info["file_name"] = ["test.fits"] + file_info["file_id"] = [0] + file_info["particle_id"] = [b"1"] + file_info["zenith"] = [20.0] * u.deg + file_info["azimuth"] = [0.0] * u.deg + file_info["nsb_level"] = [b"0.24053832149999996"] + file_info["energy_min"] = [0.003] * u.TeV + file_info["energy_max"] = [330.0] * u.TeV + file_info["viewcone_min"] = [0.0] * u.deg + file_info["viewcone_max"] = [10.0] * u.deg + file_info["core_scatter_min"] = [0.0] * u.m + file_info["core_scatter_max"] = [1900.0] * u.m + + info = reader.get_reduced_simulation_file_info(file_info) + + assert info["primary_particle"] == "gamma" + assert info["nsb_level"] == pytest.approx(0.24) + + +def test_convert_numeric_column_decodes_object_bytes(mock_fits_file): + """Test object arrays with byte entries are decoded and converted.""" + reader = EventDataReader(mock_fits_file) + + values = np.array([b" 1.5 ", b"2.0"], dtype=object) + result = reader._convert_numeric_column(values, key="nsb_level", dtype=np.float64) + + np.testing.assert_allclose(result, np.array([1.5, 2.0])) + + +def test_convert_numeric_column_raises_on_invalid_values(mock_fits_file): + """Test conversion helper raises clear error message for invalid numeric data.""" + reader = EventDataReader(mock_fits_file) + + values = np.array(["not-a-number"], dtype=np.str_) + + with pytest.raises( + ValueError, + match=r"Unable to convert FILE_INFO column 'nsb_level' to numeric values\.", + ): + reader._convert_numeric_column(values, key="nsb_level", dtype=np.float64) + + def test_get_triggered_shower_data_single_match(mock_fits_file): """Test _get_triggered_shower_data with single matches.""" reader = EventDataReader(mock_fits_file) diff --git a/tests/unit_tests/sim_events/test_writer.py b/tests/unit_tests/sim_events/test_writer.py index 52308676fb..a6fc10035c 100644 --- a/tests/unit_tests/sim_events/test_writer.py +++ b/tests/unit_tests/sim_events/test_writer.py @@ -303,6 +303,42 @@ def test_get_nsb_level_from_file_name_invalid_input(lookup_table_generator): lookup_table_generator._get_nsb_level_from_file_name(None) +def test_get_nsb_level_from_metadata_casts_string_value(mocker, lookup_table_generator): + """Test NSB metadata string values are converted to float.""" + mocker.patch( + "simtools.sim_events.writer.read_sim_telarray_metadata", + return_value=({"nsb_integrated_flux": "0.24053832149999996"}, {}), + ) + + nsb = lookup_table_generator.get_nsb_level_from_sim_telarray_metadata("dummy.simtel.zst") + + assert nsb == pytest.approx(0.24053832149999996) + + +def test_get_nsb_level_from_metadata_invalid_value_falls_back( + mocker, lookup_table_generator, caplog +): + """Test invalid NSB metadata values trigger fallback to filename parsing.""" + mocker.patch( + "simtools.sim_events.writer.read_sim_telarray_metadata", + return_value=({"nsb_integrated_flux": "not-a-number"}, {}), + ) + fallback = mocker.patch.object( + lookup_table_generator, + "_get_nsb_level_from_file_name", + return_value=0.24, + ) + + with caplog.at_level("WARNING"): + nsb = lookup_table_generator.get_nsb_level_from_sim_telarray_metadata( + "dummy_dark.simtel.zst" + ) + + fallback.assert_called_once_with("dummy_dark.simtel.zst") + assert "Invalid nsb_integrated_flux value 'not-a-number'" in caplog.text + assert nsb == pytest.approx(0.24) + + def test_process_mc_event(lookup_table_generator): """Test processing of MC events.""" lookup_table_generator.n_use = 2