Add production limits for list of telescopes

docs/changes/2128.feature.md

@@ -0,0 +1 @@
+Allow custom list of telescopes for corsika limits derivation.

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

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():

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)

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.

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):
         """

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])

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)

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)

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