diff --git a/docs/changes/2185.feature.md b/docs/changes/2185.feature.md new file mode 100644 index 0000000000..bc503ffd63 --- /dev/null +++ b/docs/changes/2185.feature.md @@ -0,0 +1 @@ +Improve derivation of CORSIKA limits and allow to define a loss fraction per energy bin ('differential loss'). diff --git a/docs/changes/2196.feature.md b/docs/changes/2196.feature.md new file mode 100644 index 0000000000..7391b8458d --- /dev/null +++ b/docs/changes/2196.feature.md @@ -0,0 +1 @@ +Add application to plot CORSIKA limits. Rearrange limits plotting for `production_merge_corsika_limits.py`. diff --git a/docs/source/api-reference/production_configuration.md b/docs/source/api-reference/production_configuration.md index a1f7aef886..7f1744570f 100644 --- a/docs/source/api-reference/production_configuration.md +++ b/docs/source/api-reference/production_configuration.md @@ -61,17 +61,6 @@ the calculation of the number of events to be simulated given a pre-determined m :members: ``` -(plot-production-grid)= - -## plot_production_grid - -```{eval-rst} -.. automodule:: production_configuration.plot_production_grid - :members: -``` - -(merge-corsika-limits)= - ## merge_corsika_limits ```{eval-rst} diff --git a/docs/source/api-reference/visualization.md b/docs/source/api-reference/visualization.md index 45caf71a98..a4cb14f40b 100644 --- a/docs/source/api-reference/visualization.md +++ b/docs/source/api-reference/visualization.md @@ -121,3 +121,24 @@ the visualization module. .. automodule:: visualization.plot_event_level_production_comparison :members: ``` + +(plot-production-grid)= + +## plot_production_grid + +```{eval-rst} +.. automodule:: visualization.plot_production_grid + :members: +``` + +(merge-corsika-limits)= + + +(plot-corsika-limits)= + +## plot_corsika_limits + +```{eval-rst} +.. automodule:: visualization.plot_corsika_limits + :members: +``` diff --git a/docs/source/user-guide/applications.md b/docs/source/user-guide/applications.md index b8d38cacba..0f70929518 100644 --- a/docs/source/user-guide/applications.md +++ b/docs/source/user-guide/applications.md @@ -94,6 +94,7 @@ simtools-production-derive-corsika-limits simtools-production-generate-grid simtools-production-merge-corsika-limits +simtools-production-plot-corsika-limits simtools-run-application simtools-simulate-flasher simtools-simulate-illuminator diff --git a/docs/source/user-guide/applications/simtools-production-plot-corsika-limits.rst b/docs/source/user-guide/applications/simtools-production-plot-corsika-limits.rst new file mode 100644 index 0000000000..6a8f1ce125 --- /dev/null +++ b/docs/source/user-guide/applications/simtools-production-plot-corsika-limits.rst @@ -0,0 +1,5 @@ +simtools-production-plot-corsika-limits +======================================= + +.. automodule:: production_plot_corsika_limits + :members: diff --git a/pyproject.toml b/pyproject.toml index d253530718..fc76d05bae 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -118,6 +118,7 @@ scripts.simtools-production-derive-corsika-limits = "simtools.applications.produ scripts.simtools-production-derive-statistics = "simtools.applications.production_derive_statistics:main" scripts.simtools-production-generate-grid = "simtools.applications.production_generate_grid:main" scripts.simtools-production-merge-corsika-limits = "simtools.applications.production_merge_corsika_limits:main" +scripts.simtools-production-plot-corsika-limits = "simtools.applications.production_plot_corsika_limits:main" scripts.simtools-run-application = "simtools.applications.run_application:main" scripts.simtools-simulate-flasher = "simtools.applications.simulate_flasher:main" scripts.simtools-simulate-illuminator = "simtools.applications.simulate_illuminator:main" diff --git a/src/simtools/applications/plot_production_grid.py b/src/simtools/applications/plot_production_grid.py index 7a067f957a..f1500dcae6 100644 --- a/src/simtools/applications/plot_production_grid.py +++ b/src/simtools/applications/plot_production_grid.py @@ -52,7 +52,7 @@ from simtools.application_control import build_application from simtools.model.site_model import SiteModel -from simtools.production_configuration.plot_production_grid import ProductionGridPlotter +from simtools.visualization.plot_production_grid import ProductionGridPlotter logger = logging.getLogger(__name__) diff --git a/src/simtools/applications/production_derive_corsika_limits.py b/src/simtools/applications/production_derive_corsika_limits.py index 3fe9075834..2936c0a856 100644 --- a/src/simtools/applications/production_derive_corsika_limits.py +++ b/src/simtools/applications/production_derive_corsika_limits.py @@ -3,41 +3,74 @@ r""" Derive CORSIKA configuration limits for energy, core distance, and viewcone radius. -This tool determines configuration limits based on triggered events from broad-range +This tool derives configuration limits from reduced event data generated by broad-range simulations. It supports the derivation of the following CORSIKA configuration parameters: - **ERANGE**: lower energy limit - **CSCAT**: upper core distance - **VIEWCONE**: viewcone radius -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. +Broad-range simulations in this context are simulation sets generated with wide ranges +for these parameters. + +How limits are derived +---------------------- +Lower energy limit (**ERANGE**) + Derived from the 1D triggered-energy histogram using a peak-relative threshold: + + 1. Find the absolute maximum bin in the triggered-energy histogram. + 2. Compute a stable peak by averaging the maximum bin and available immediate neighbors. + 3. Define a threshold as fraction of this stable peak. + 4. Walk toward lower energies from the peak to the first bin below threshold. + 5. Use that bin's lower edge as the energy limit. + +Upper core distance (**CSCAT**) and viewcone radius (**VIEWCONE**) + Derived from per-axis allowed-loss settings (``axis,fraction,min_events``), either as + integrated limits or per-energy-bin + differential limits (depending on ``--differential_loss_bins_per_decade``). + +energy_threshold_fraction (float, optional) + Fraction of the stable energy-peak count used to derive ERANGE. + +``--allowed_losses`` + Required for core_distance and angular_distance, or ``all`` plus overrides. + 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 @@ -57,14 +90,23 @@ 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. +allowed_losses (str, required, repeatable) + Per-axis allowed-loss tuple in the form + ``axis,fraction,min_events``. + Use once per axis (core_distance, angular_distance), or use ``all`` + to set both axes and optionally override selected axes. + Core distance and angular distance limits use these settings directly. +energy_threshold_fraction (float, optional) + Fraction of the stable energy-peak count used to derive ERANGE (default: 0.01). plot_histograms (bool, optional) Plot histograms of the event data. output_file (str, optional) Path to the output file for the derived limits. n_workers (int, optional) Number of worker processes to use for execution. Default is 1. +differential_loss_bins_per_decade (int, optional) + Number of differential energy bins per decade for per-bin limit computation. + Set to 0 (default) to use integrated limits. Example ------- @@ -76,18 +118,21 @@ simtools-production-derive-corsika-limits \\ --event_data_file event_dat_file.hdf5 \\ --array_layout_name alpha,beta \\ - --loss_fraction 1e-6 \\ + --allowed_losses core_distance,1e-6,10 \ + --allowed_losses angular_distance,1e-6,10 \ + --energy_threshold_fraction 0.01 \ --plot_histograms \\ --output_file corsika_simulation_limits.ecsv -Derive limits for a single production with a given file for custom defined array layouts: +Derive limits for a single production with a custom telescope-config file: .. code-block:: console simtools-production-derive-corsika-limits \\ --event_data_file event_dat_file.hdf5 \\ --telescope_ids path/to/telescope_configs.yaml \\ - --loss_fraction 1e-6 \\ + --allowed_losses all,1e-6,10 \ + --energy_threshold_fraction 0.01 \ --plot_histograms \\ --output_file corsika_simulation_limits.ecsv @@ -99,7 +144,8 @@ --event_data_file pattern_1_*.hdf5 \\ --event_data_file pattern_2_*.hdf5 \\ --array_layout_name alpha \\ - --loss_fraction 1e-6 \\ + --allowed_losses all,1e-6,10 \ + --energy_threshold_fraction 0.01 \ --plot_histograms \\ --n_workers 4 \\ --output_file corsika_simulation_limits.ecsv @@ -130,10 +176,25 @@ 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.", + 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 core_distance, angular_distance, " + "or use all to set both. " + "Example: --allowed_losses core_distance,1e-6,10" + ), + ) + parser.add_argument( + "--energy_threshold_fraction", + help="Fraction of the stable energy-peak count used to derive ERANGE ", + type=parser.efficiency_interval, + required=False, + default=0.01, ) parser.add_argument( "--plot_histograms", @@ -151,6 +212,16 @@ def _add_arguments(parser): required=False, default=1, ) + parser.add_argument( + "--differential_loss_bins_per_decade", + help=( + "Number of differential energy bins per decade for per-bin limit computation. " + "Set to 0 (default) to use integrated limits." + ), + type=int, + required=False, + default=0, + ) def main(): diff --git a/src/simtools/applications/production_merge_corsika_limits.py b/src/simtools/applications/production_merge_corsika_limits.py index a6ce1a1454..dd50d7dbfb 100644 --- a/src/simtools/applications/production_merge_corsika_limits.py +++ b/src/simtools/applications/production_merge_corsika_limits.py @@ -13,10 +13,9 @@ This tool supports three main use cases: 1. Merge multiple CORSIKA limit tables into a single file and optionally generate - plots of the derived limits. + a merged limits table. 2. Merge tables and also check for grid completeness against a provided grid - definition file. This requires the --grid_definition parameter. Coverage plots - can also be generated. + definition file. This requires the --grid_definition parameter. 3. Check grid completeness of an already merged table file. This requires both the --merged_table and --grid_definition parameters. @@ -36,10 +35,6 @@ output_file (str, optional) Name of the output file for the merged limits table. Default is "merged_corsika_limits.ecsv". -plot_grid_coverage (bool, optional) - Flag to generate plots showing grid coverage. Requires --grid_definition. -plot_limits (bool, optional) - Flag to generate plots showing the derived limits. Examples -------- @@ -49,7 +44,7 @@ simtools-production-merge-corsika-limits \\ --input_files "simtools-output/corsika_limits/" \\ - --output_file merged_limits.ecsv --plot_limits + --output_file merged_limits.ecsv 2. Merge tables and check grid completeness: @@ -58,7 +53,7 @@ simtools-production-merge-corsika-limits \\ --input_files "simtools-output/corsika_limits/" \\ --grid_definition grid_definition.yaml \\ - --output_file merged_limits.ecsv --plot_grid_coverage + --output_file merged_limits.ecsv 3. Check grid completeness of an existing merged table: @@ -66,7 +61,7 @@ simtools-production-merge-corsika-limits \\ --merged_table merged_limits.ecsv \\ - --grid_definition grid_definition.yaml --plot_grid_coverage + --grid_definition grid_definition.yaml 4. Merge tables using a list of files from a text file: @@ -114,18 +109,6 @@ def _add_arguments(parser): default=None, help="Path to YAML file defining the expected grid points.", ) - parser.add_argument( - "--plot_grid_coverage", - help="Generate plots showing grid coverage.", - action="store_true", - default=False, - ) - parser.add_argument( - "--plot_limits", - help="Generate plots showing the derived limits.", - action="store_true", - default=False, - ) def main(): diff --git a/src/simtools/applications/production_plot_corsika_limits.py b/src/simtools/applications/production_plot_corsika_limits.py new file mode 100644 index 0000000000..c4bab97585 --- /dev/null +++ b/src/simtools/applications/production_plot_corsika_limits.py @@ -0,0 +1,51 @@ +#!/usr/bin/python3 + +r""" +Plot CORSIKA limits from a ECSV table. + +This application reads a CORSIKA limits table and plots the limits +as function of zenith angle. + + +Command line arguments +---------------------- +input (str, required) + Path to a CORSIKA limits table in ECSV format. + +Example +------- + +.. code-block:: console + + simtools-production-plot-corsika-limits \ + --input simtools-output/merged_corsika_limits.ecsv \ + --output_path simtools-output +""" + +from simtools.application_control import build_application +from simtools.data_model import data_reader +from simtools.visualization.plot_corsika_limits import plot_limits + + +def _add_arguments(parser): + """Register application-specific command line arguments.""" + parser.add_argument( + "--input", + type=str, + required=True, + help="Path to a merged CORSIKA limits table in ECSV format.", + ) + + +def main(): + """Run CORSIKA limits plotting.""" + app_context = build_application(initialization_kwargs={"output": True}) + + plot_limits( + data_reader.read_table_from_file(app_context.args["input"]), + app_context.io_handler.get_output_directory(), + ) + + +if __name__ == "__main__": + main() diff --git a/src/simtools/production_configuration/derive_corsika_limits.py b/src/simtools/production_configuration/derive_corsika_limits.py index d3e5fe234e..70706ecf67 100644 --- a/src/simtools/production_configuration/derive_corsika_limits.py +++ b/src/simtools/production_configuration/derive_corsika_limits.py @@ -21,18 +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 = ("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", "viewcone_radius", + "br_energy_min", + "br_energy_max", + "br_core_scatter_max", + "br_viewcone_max", ] @@ -122,9 +138,11 @@ 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"] + allowed_losses = job_spec["allowed_losses"] + energy_threshold_fraction = job_spec["energy_threshold_fraction"] 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) _logger.info( f"Processing production {production_index}: pattern={production_pattern}, " @@ -135,9 +153,11 @@ def _execute_production_job(job_spec): production_pattern, array_name, telescope_ids, - loss_fraction, + allowed_losses, + energy_threshold_fraction, plot_histograms, output_subdir=output_subdir, + differential_loss_bins_per_decade=differential_loss_bins_per_decade, ) result.update( @@ -171,6 +191,71 @@ def _resolve_telescope_configs(args_dict): ) +def _parse_allowed_losses(allowed_losses_args): + """ + Parse repeatable --allowed_losses values for core/viewcone axes. + + 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( + f"Invalid axis for --allowed_losses. Allowed axes: {', '.join(LOSS_AXES)}, 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: @@ -201,6 +286,9 @@ 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")) + energy_threshold_fraction = float(args_dict.get("energy_threshold_fraction", 0.01)) + 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 @@ -231,9 +319,11 @@ 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"], + "allowed_losses": allowed_losses, + "energy_threshold_fraction": energy_threshold_fraction, "plot_histograms": args_dict["plot_histograms"], "output_subdir": output_subdir, + "differential_loss_bins_per_decade": differential_loss_bins_per_decade, } job_specs.append(job_spec) @@ -245,16 +335,18 @@ def generate_corsika_limits_grid(args_dict): max_workers=n_workers, ) - write_results(results, args_dict) + write_results(results, args_dict, allowed_losses, energy_threshold_fraction) def _process_file( file_path, array_name, telescope_ids, - loss_fraction, - plot_histograms, + allowed_losses, + energy_threshold_fraction=0.01, + plot_histograms=False, output_subdir=None, + differential_loss_bins_per_decade=0, ): """ Compute limits for a given event data file and telescope configuration. @@ -269,12 +361,17 @@ 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. + allowed_losses : dict + Per-axis loss settings for core_distance/angular_distance. + energy_threshold_fraction : float, optional + Fraction of the stable energy-peak count used to derive ERANGE. plot_histograms : bool Whether to plot histograms. output_subdir : Path or None, optional Output subdirectory for plots. If None, uses default output directory. + differential_loss_bins_per_decade : int, optional + Number of energy bins per decade for differential per-bin limits. + Set to 0 (default) to use integrated limits. Returns ------- @@ -283,17 +380,32 @@ def _process_file( """ histograms = EventDataHistograms( file_path, - array_name=array_name, - telescope_list=telescope_ids, + array_name, + telescope_ids, + differential_loss_bins_per_decade or 10, ) - histograms.fill() + histograms.fill(fill_efficiency_histogram=False) limits = { - "lower_energy_limit": compute_lower_energy_limit(histograms, loss_fraction), - "upper_radius_limit": compute_upper_radius_limit(histograms, loss_fraction), - "viewcone_radius": compute_viewcone(histograms, loss_fraction), + "lower_energy_limit": compute_lower_energy_limit( + histograms, + energy_threshold_fraction, + ), } + limits.update( + _compute_limits( + histograms, + 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() @@ -307,7 +419,132 @@ def _process_file( return limits -def write_results(results, args_dict): +def _compute_limits(histograms, allowed_losses, bins_per_decade): + """ + Compute core and viewcone limits per energy bin and return max limits. + + 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 = {} + differential_energy_bins = None + 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)))) + differential_energy_bins = np.logspace(low, high, (high - low) * bins_per_decade + 1) + + for axis_name, config in axis_configs.items(): + if bins_per_decade > 0: + axis_max, curve_x, curve_y = _differential_upper_limits( + histograms.histograms[f"{axis_name}_vs_energy"]["histogram"], + config["x_bins"], + histograms.energy_bins, + differential_energy_bins, + 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, + "Upper radius limit is equal to the maximum core scatter distance of", + ) + viewcone_radius = _is_close( + vc_max * u.deg, + 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", + ) + _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": upper_radius_limit, + "viewcone_radius": viewcone_radius, + 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" + ], + } + + +def _differential_upper_limits( + histogram2d, + x_bins, + y_bins, + diff_e_bins, + allowed_loss, + name, + unit, +): + """Compute upper limits per energy slice of a 2D (x, energy) histogram.""" + y_centers = 0.5 * (y_bins[:-1] + y_bins[1:]) + limits, energy_centers = [], [] + n = len(diff_e_bins) - 1 + for i in range(n): + e_low, e_high = diff_e_bins[i], diff_e_bins[i + 1] + hi_op = np.less_equal if i == n - 1 else np.less + projected = np.sum(histogram2d[:, (y_centers >= e_low) & hi_op(y_centers, e_high)], axis=1) + total = float(np.sum(projected)) + if total <= 0: + continue + limit = _integral_limits( + projected, + x_bins, + allowed_loss["loss_fraction"], + allowed_loss["loss_min_events"], + limit_type="upper", + ) + keep = np.searchsorted(x_bins, limit, side="left") + loss = (total - float(np.sum(projected[:keep]))) / total + limits.append(limit) + energy_centers.append(float(np.sqrt(e_low * e_high))) + _logger.info( + f"Differential {name}: E=[{e_low:.4g}, {e_high:.4g}] TeV, " + f"N={int(total)}, limit={limit:.4g} {unit}, loss={loss:.5f}" + ) + return ( + (float(np.max(limits)), limits, energy_centers) if limits else (float(x_bins[-1]), [], []) + ) + + +def write_results(results, args_dict, allowed_losses, energy_threshold_fraction): """ Write the computed limits as astropy table to file. @@ -317,8 +554,16 @@ def write_results(results, args_dict): List of computed limits. args_dict : dict Dictionary containing command line arguments. + allowed_losses : dict + Per-axis loss settings for core_distance/angular_distance. + energy_threshold_fraction : float + Fraction used for deriving the lower energy threshold. """ - table = _create_results_table(results, args_dict["loss_fraction"]) + table = _create_results_table( + results, + allowed_losses, + energy_threshold_fraction, + ) output_dir = io_handler.IOHandler().get_output_directory() output_file = output_dir / args_dict["output_file"] @@ -329,7 +574,7 @@ def write_results(results, args_dict): MetadataCollector.dump(args_dict, output_file) -def _create_results_table(results, loss_fraction): +def _create_results_table(results, allowed_losses, energy_threshold_fraction): """ Convert list of simulation results to an astropy Table with metadata. @@ -339,8 +584,10 @@ def _create_results_table(results, loss_fraction): ---------- results : list[dict] Computed limits per file and telescope configuration. - loss_fraction : float - Fraction of lost events (added to metadata). + allowed_losses : dict + Per-axis loss settings added to metadata. + energy_threshold_fraction : float + Fraction used for deriving the lower energy threshold. Returns ------- @@ -362,9 +609,14 @@ def _create_results_table(results, loss_fraction): { "created": datetime.datetime.now().isoformat(), "description": "Lookup table for CORSIKA limits computed from simulations.", - "loss_fraction": loss_fraction, } ) + 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"] + ) + table.meta["energy_threshold_fraction"] = energy_threshold_fraction return table @@ -402,17 +654,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, 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. + 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. @@ -424,6 +683,8 @@ def _compute_limits(hist, bin_edges, loss_fraction, limit_type="lower"): Array of bin edges. loss_fraction : float Fraction of events to be lost. + loss_min_events : int, optional + Minimum number of events to be lost after applying a limit. limit_type : str, optional Type of limit ('lower' or 'upper'). Default is 'lower'. @@ -433,28 +694,84 @@ def _compute_limits(hist, bin_edges, loss_fraction, limit_type="lower"): Bin edge value corresponding to the threshold. """ total_events = np.sum(hist) - threshold = (1 - loss_fraction) * total_events + # Keep-threshold corresponding to a strictly greater-than requested absolute event loss. + max_kept_for_min_loss = np.nextafter(total_events - float(loss_min_events), -np.inf) + threshold = min( + (1 - loss_fraction) * total_events, + max_kept_for_min_loss, + ) + threshold = np.clip(threshold, 0.0, total_events) if limit_type == "upper": cum = np.cumsum(hist) idx = np.searchsorted(cum, threshold) + 1 return bin_edges[min(idx, len(bin_edges) - 1)] if limit_type == "lower": cum = np.cumsum(hist[::-1]) - idx = np.searchsorted(cum, threshold) + 1 + idx = np.searchsorted(cum, threshold) return bin_edges[max(len(bin_edges) - 1 - idx, 0)] raise ValueError("limit_type must be 'lower' or 'upper'") -def compute_lower_energy_limit(histograms, loss_fraction): +def _find_low_energy_threshold_from_histogram(counts, bin_edges, threshold_fraction=0.1): + """Find low-energy threshold from a 1D histogram using a peak-relative criterion. + + The threshold is defined as the first bin (walking to lower energies from the + histogram maximum) where the count drops below ``threshold_fraction`` times a + stable peak estimate. The stable peak estimate is the mean of the maximum bin + and its immediate neighbors (neighbors included only when available). + + Parameters + ---------- + counts : np.ndarray + Histogram bin counts. + bin_edges : np.ndarray + Histogram bin edges (length must be ``len(counts) + 1``). + threshold_fraction : float, optional + Fraction of the stable peak used as threshold. + + Returns + ------- + float + Derived energy threshold. + """ + counts = np.asarray(counts, dtype=float) + bin_edges = np.asarray(bin_edges, dtype=float) + + if counts.ndim != 1 or bin_edges.ndim != 1: + raise ValueError("counts and bin_edges must be one-dimensional arrays") + if counts.size == 0: + raise ValueError("counts must not be empty") + if not np.any(counts > 0): + raise ValueError("counts must contain at least one positive entry") + if bin_edges.size != counts.size + 1: + raise ValueError("bin_edges length must be len(counts) + 1") + if not 0.0 < threshold_fraction <= 1.0: + raise ValueError("threshold_fraction must be in the interval (0, 1]") + + peak_idx = int(np.argmax(counts)) + + left = max(peak_idx - 1, 0) + right = min(peak_idx + 1, counts.size - 1) + n_peak = float(np.mean(counts[left : right + 1])) + threshold = threshold_fraction * n_peak + + for idx in range(peak_idx, -1, -1): + if counts[idx] < threshold: + return float(bin_edges[idx]) + + return float(bin_edges[0]) + + +def compute_lower_energy_limit(histograms, threshold_fraction): """ - Compute the lower energy limit in TeV based on the event loss fraction. + Compute the lower energy limit in TeV based on the threshold fraction. Parameters ---------- histograms : EventDataHistograms Histograms. - loss_fraction : float - Fraction of events to be lost. + threshold_fraction : float + Fraction of the stable peak used as threshold. Returns ------- @@ -462,11 +779,10 @@ def compute_lower_energy_limit(histograms, loss_fraction): Lower energy limit. """ energy_min = ( - _compute_limits( + _find_low_energy_threshold_from_histogram( histograms.histograms["energy"]["histogram"], histograms.energy_bins, - loss_fraction, - limit_type="lower", + threshold_fraction=threshold_fraction, ) * u.TeV ) @@ -485,75 +801,3 @@ def _is_close(value, reference, warning_text): if reference is not None and np.isclose(value.value, reference.value, rtol=1.0e-2): _logger.warning(f"{warning_text} {value}.") return value - - -def compute_upper_radius_limit(histograms, loss_fraction): - """ - Compute the upper radial distance based on the event loss fraction. - - Parameters - ---------- - histograms : EventDataHistograms - Histograms. - loss_fraction : float - Fraction of events to be lost. - - Returns - ------- - astropy.units.Quantity - Upper radial distance in m. - """ - radius_limit = ( - _compute_limits( - histograms.histograms["core_distance"]["histogram"], - histograms.core_distance_bins, - loss_fraction, - limit_type="upper", - ) - * u.m - ) - return _is_close( - radius_limit, - 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", - ) - - -def compute_viewcone(histograms, loss_fraction): - """ - 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. - - Returns - ------- - astropy.units.Quantity - Viewcone radius in degrees. - """ - viewcone_limit = ( - _compute_limits( - histograms.histograms["angular_distance"]["histogram"], - histograms.view_cone_bins, - loss_fraction, - 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", - ) diff --git a/src/simtools/production_configuration/merge_corsika_limits.py b/src/simtools/production_configuration/merge_corsika_limits.py index 6afc2aeabe..ddc9ad10df 100644 --- a/src/simtools/production_configuration/merge_corsika_limits.py +++ b/src/simtools/production_configuration/merge_corsika_limits.py @@ -4,7 +4,6 @@ from itertools import product from pathlib import Path -import matplotlib.pyplot as plt import numpy as np from astropy.table import unique, vstack @@ -15,8 +14,6 @@ _logger = logging.getLogger(__name__) -ZENITH_LABEL = "Zenith [deg]" - class CorsikaMergeLimits: """Class for merging CORSIKA limit tables and checking grid completeness.""" @@ -310,179 +307,6 @@ def check_grid_completeness(self, merged_table, grid_definition): "expected_str": expected_combinations_str, } - def _plot_single_grid_coverage( - self, ax, zeniths, azimuths, nsb, array_name, found_combinations_str - ): - """Plot grid coverage for a single NSB and array_name.""" - z_grid = np.zeros((len(zeniths), len(azimuths))) - for i, zenith in enumerate(zeniths): - for j, azimuth in enumerate(azimuths): - point_str = (str(zenith), str(azimuth), str(nsb), str(array_name)) - if point_str in found_combinations_str: - z_grid[i, j] = 1 - - az_vals = azimuths.value if hasattr(azimuths, "value") else azimuths - zen_vals = zeniths.value if hasattr(zeniths, "value") else zeniths - extent = [ - min(az_vals) - 0.5, - max(az_vals) + 0.5, - max(zen_vals) + 0.5, - min(zen_vals) - 0.5, - ] - colors = ["red", "green"] - cmap = plt.matplotlib.colors.ListedColormap(colors) - im = ax.imshow(z_grid, cmap=cmap, vmin=0, vmax=1, extent=extent) - - cbar = plt.colorbar( - im, - ax=ax, - ticks=[0, 1], - label="Coverage", - shrink=0.25, - pad=0.02, - ) - cbar.set_ticklabels(["Missing", "Present"]) - ax.set_title(f"Grid Coverage: NSB={nsb}, Array Name={array_name}") - ax.set_xlabel("Azimuth [deg]") - ax.set_ylabel(ZENITH_LABEL) - ax.set_xticks(az_vals) - ax.set_yticks(zen_vals) - ax.grid(which="major", linestyle="-", linewidth="0.5", color="black", alpha=0.3) - - def plot_grid_coverage(self, merged_table, grid_definition): - """Generate plots showing grid coverage for each combination of NSB level and array name. - - Creates a series of heatmap plots showing which grid points (combinations of zenith and - azimuth angles) are present or missing in the merged table, for each combination of - NSB level and array name. - - Parameters - ---------- - merged_table : astropy.table.Table - The merged table containing CORSIKA limit data. - grid_definition : dict - Dictionary defining the grid dimensions with keys: - 'zenith': list of zenith angles, - 'azimuth': list of azimuth angles, - 'nsb_level': list of NSB levels, - 'array_name': list of array names - - Returns - ------- - list - List of Path objects pointing to the saved plot files. - """ - if not grid_definition: - _logger.info("No grid definition provided, skipping grid coverage plots.") - return [] - - _logger.info("Generating grid coverage plots") - output_files = [] - - _, completeness_info = self.check_grid_completeness(merged_table, grid_definition) - found_combinations_str = completeness_info.get("found_str", set()) - - unique_values = { - "zeniths": np.array(grid_definition.get("zenith", [])), - "azimuths": np.array(grid_definition.get("azimuth", [])), - "nsb_levels": np.array(grid_definition.get("nsb_level", [])), - "array_names": np.array(grid_definition.get("array_name", [])), - } - - for nsb, array_name in product(unique_values["nsb_levels"], unique_values["array_names"]): - _, ax = plt.subplots(figsize=(10, 8)) - self._plot_single_grid_coverage( - ax, - unique_values["zeniths"], - unique_values["azimuths"], - nsb, - array_name, - found_combinations_str, - ) - output_file = self.output_dir / f"grid_coverage_{nsb}_{array_name}.png" - plt.tight_layout() - plt.savefig(output_file, bbox_inches="tight") - plt.close() - output_files.append(output_file) - return output_files - - def plot_limits(self, merged_table): - """Create plots showing the derived limits for each combination of array_name and azimuth. - - Creates plots showing the lower energy limit, upper radius limit, and viewcone radius - versus zenith angle for each combination of array_name and azimuth angle. Each plot has - lines for different NSB levels. - - Parameters - ---------- - merged_table : astropy.table.Table - The merged table containing CORSIKA limit data. - - Returns - ------- - list - List of Path objects pointing to the saved plot files. - """ - _logger.info("Generating limit plots") - output_files = [] - - grouped_by_layout_az = merged_table.group_by(["array_name", "azimuth"]) - - for group in grouped_by_layout_az.groups: - array_name = group["array_name"][0] - azimuth = group["azimuth"][0] - azimuth_value = azimuth.value if hasattr(azimuth, "value") else azimuth - - fig, axes = plt.subplots(1, 3, figsize=(18, 6)) - legend_handles, legend_labels = [], [] - - grouped_by_nsb = group.group_by("nsb_level") - colors = plt.get_cmap("viridis")(np.linspace(0, 1, len(grouped_by_nsb.groups))) - - for i, nsb_group in enumerate(grouped_by_nsb.groups): - nsb_level = nsb_group["nsb_level"][0] - plot_columns = [ - "zenith", - "lower_energy_limit", - "upper_radius_limit", - "viewcone_radius", - ] - agg_data = nsb_group[plot_columns].group_by("zenith").groups.aggregate(np.mean) - agg_data.sort("zenith") - zeniths = agg_data["zenith"].value - - (line,) = axes[0].plot( - zeniths, agg_data["lower_energy_limit"], "o-", color=colors[i] - ) - axes[1].plot(zeniths, agg_data["upper_radius_limit"], "o-", color=colors[i]) - axes[2].plot(zeniths, agg_data["viewcone_radius"], "o-", color=colors[i]) - legend_handles.append(line) - legend_labels.append(f"NSB={nsb_level}") - - axes[0].set_title("Lower Energy Limit vs Zenith") - axes[0].set_xlabel(ZENITH_LABEL) - axes[0].set_ylabel("Lower Energy Limit [TeV]") - axes[0].grid(True) - axes[1].set_title("Upper Radius Limit vs Zenith") - axes[1].set_xlabel(ZENITH_LABEL) - axes[1].set_ylabel("Upper Radius Limit [m]") - axes[1].grid(True) - axes[2].set_title("Viewcone Radius vs Zenith") - axes[2].set_xlabel(ZENITH_LABEL) - axes[2].set_ylabel("Viewcone Radius [deg]") - axes[2].grid(True) - - fig.legend(legend_handles, legend_labels, loc="lower center", ncol=len(legend_labels)) - plt.suptitle(f"CORSIKA Limits: Array Name={array_name}, Azimuth={azimuth_value} deg") - plt.tight_layout() - plt.subplots_adjust(bottom=0.15) - - output_file = self.output_dir / f"limits_{array_name}_azimuth{azimuth_value}.png" - plt.savefig(output_file) - plt.close(fig) - output_files.append(output_file) - return output_files - def write_merged_table(self, merged_table, output_file, input_files, grid_completeness): """Write the merged table to file and save metadata. @@ -567,7 +391,7 @@ def resolve_input_files_and_table(args_dict, merger): def merge_corsika_limits(args_dict, merger=None): """ - Run table merge, completeness checks, optional plotting, and optional write-out. + Run table merge, completeness checks, and optional write-out. Parameters ---------- @@ -584,12 +408,6 @@ def merge_corsika_limits(args_dict, merger=None): is_complete, grid_completeness = merger.check_grid_completeness(merged_table, grid_definition) - if args_dict.get("plot_grid_coverage"): - merger.plot_grid_coverage(merged_table, grid_definition) - - if args_dict.get("plot_limits"): - merger.plot_limits(merged_table) - if not from_merged_table: output_file = merger.output_dir / args_dict["output_file"] merger.write_merged_table( diff --git a/src/simtools/sim_events/histograms.py b/src/simtools/sim_events/histograms.py index 5084104034..87ddc8ec95 100644 --- a/src/simtools/sim_events/histograms.py +++ b/src/simtools/sim_events/histograms.py @@ -33,14 +33,23 @@ class EventDataHistograms: Name of the telescope array configuration (default is None). telescope_list : list, optional List of telescope IDs to filter the events (default is None). + energy_bins_per_decade : int, optional + Number of energy bins per decade for logarithmic energy histograms. """ - def __init__(self, event_data_file, array_name=None, telescope_list=None): + def __init__( + self, + event_data_file, + array_name=None, + telescope_list=None, + energy_bins_per_decade=10, + ): """Initialize.""" self._logger = logging.getLogger(__name__) self.event_data_file = event_data_file self.event_data_files = self._normalize_event_data_files(event_data_file) self.array_name = array_name + self.energy_bins_per_decade = max(int(energy_bins_per_decade), 1) self.histograms = {} self.file_info = {} @@ -96,6 +105,7 @@ def _update_file_info(self, file_info_table): "azimuth": self._get_file_info_value(file_info_table, "azimuth", "deg"), "nsb_level": self._get_file_info_value(file_info_table, "nsb_level"), "energy_min": self._get_file_info_value(file_info_table, "energy_min", "TeV"), + "energy_max": self._get_file_info_value(file_info_table, "energy_max", "TeV"), "core_scatter_max": self._get_file_info_value(file_info_table, "core_scatter_max", "m"), "viewcone_max": self._get_file_info_value(file_info_table, "viewcone_max", "deg"), "solid_angle": self._get_file_info_value(file_info_table, "solid_angle", "sr"), @@ -115,7 +125,7 @@ def _fill_current_histograms(self): for data in self.histograms.values(): self._fill_histogram_and_bin_edges(data) - def fill(self): + def fill(self, fill_efficiency_histogram=True): """ Fill histograms with event data. @@ -124,6 +134,11 @@ def fill(self): Assume that all event data files are generated with similar configurations (self.file_info contains the file info of the last file). + + Parameters + ---------- + fill_efficiency_histogram : bool, optional + Whether to calculate and fill the efficiency histograms. """ total_files = len(self.event_data_files) for file_index, (event_data_file, reader) in enumerate(self._iter_readers(), start=1): @@ -143,7 +158,8 @@ def fill(self): self._fill_current_histograms() self.print_summary() - self.calculate_efficiency_data() + if fill_efficiency_histogram: + self.calculate_efficiency_data() self.calculate_cumulative_data() def _define_histograms(self, event_data, triggered_data, shower_data): @@ -190,12 +206,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": {"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": {"y": "log"}, }, "x_core_shower_vs_y_core_shower": { "event_data_column": ("x_core_shower", "y_core_shower"), @@ -204,7 +222,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), @@ -264,6 +282,7 @@ def get_histogram_definition( "1d": is_1d, "bin_edges": bin_edges, "title": title, + "title_fontsize": "xx-small", "axis_titles": axis_titles, "suffix": suffix, "plot_scales": plot_scales, @@ -349,20 +368,42 @@ def calculate_efficiency(trig_hist, mc_hist): @property def energy_bins(self): - """Return bins for the energy histogram.""" + """ + Return bins for the energy histogram. + + Align bins to full decades of energy, using the configured bins per decade, + and ensure that the range covers the energy range of the events. + + Returns + ------- + np.ndarray Array of energy bin edges in TeV. + """ if "energy_bin_edges" in self.histograms: return self.histograms["energy_bin_edges"] - return np.logspace( - np.log10(self.file_info.get("energy_min", 1.0e-3 * u.TeV).to("TeV").value), - np.log10(self.file_info.get("energy_max", 1.0e3 * u.TeV).to("TeV").value), - 100, - ) + + energy_min = self.file_info.get("energy_min", 1.0e-3 * u.TeV).to("TeV").value + energy_max = self.file_info.get("energy_max", 1.0e3 * u.TeV).to("TeV").value + energy_min = max(energy_min, 1e-3) + energy_max = max(energy_max, 10 * energy_min) + + lower_decade = np.floor(np.log10(energy_min)) + upper_decade = np.ceil(np.log10(energy_max)) + if upper_decade <= lower_decade: + upper_decade = lower_decade + 1 + + n_bins = int((upper_decade - lower_decade) * self.energy_bins_per_decade) + return np.logspace(lower_decade, upper_decade, n_bins + 1) @property def core_distance_bins(self): - """Return bins for the core distance histogram.""" + """ + Return bins for the core distance histogram. + + CORSIKA CSCAT ('core_scatter_max') is defined in the shower plane. + """ 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, diff --git a/src/simtools/visualization/plot_corsika_limits.py b/src/simtools/visualization/plot_corsika_limits.py new file mode 100644 index 0000000000..f58c803fcc --- /dev/null +++ b/src/simtools/visualization/plot_corsika_limits.py @@ -0,0 +1,265 @@ +"""Plotting utilities for CORSIKA limits tables.""" + +import logging +from itertools import product +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.colors import ListedColormap + +_logger = logging.getLogger(__name__) + +ZENITH_LABEL = "Zenith [deg]" +BROAD_RANGE_COLUMN_ALIASES = { + "lower_energy_limit": ["br_energy_min", "br_lower_energy_limit"], + "upper_radius_limit": ["br_core_scatter_max", "br_upper_radius_limit"], + "viewcone_radius": ["br_viewcone_max", "br_viewcone_radius"], +} + + +def _get_primary_particle_label(table): + """Return a primary particle label derived from table content.""" + if "primary_particle" not in table.colnames: + return "unknown" + + unique_particles = np.unique(np.array(table["primary_particle"], dtype=str)) + if len(unique_particles) == 1: + return unique_particles[0] + + return "/".join(unique_particles) + + +def _resolve_broad_range_columns(limits_table): + """Resolve broad-range column names from supported aliases.""" + resolved_columns = {} + for column_key, aliases in BROAD_RANGE_COLUMN_ALIASES.items(): + for alias in aliases: + if alias in limits_table.colnames: + resolved_columns[column_key] = alias + break + + if len(resolved_columns) != len(BROAD_RANGE_COLUMN_ALIASES): + _logger.warning( + "Not all broad-range columns found. Expected aliases: %s. Found: %s", + BROAD_RANGE_COLUMN_ALIASES, + resolved_columns, + ) + return None + + return resolved_columns + + +def _plot_single_grid_coverage( + ax, zeniths, azimuths, nsb, array_name, found_combinations_str, primary_particle +): + """Plot grid coverage for a single NSB and array name.""" + z_grid = np.zeros((len(zeniths), len(azimuths))) + for i, zenith in enumerate(zeniths): + for j, azimuth in enumerate(azimuths): + point_str = (str(zenith), str(azimuth), str(nsb), str(array_name)) + if point_str in found_combinations_str: + z_grid[i, j] = 1 + + az_vals = azimuths.value if hasattr(azimuths, "value") else azimuths + zen_vals = zeniths.value if hasattr(zeniths, "value") else zeniths + extent = [ + min(az_vals) - 0.5, + max(az_vals) + 0.5, + max(zen_vals) + 0.5, + min(zen_vals) - 0.5, + ] + + im = ax.imshow(z_grid, cmap=ListedColormap(["red", "green"]), vmin=0, vmax=1, extent=extent) + cbar = plt.colorbar(im, ax=ax, ticks=[0, 1], label="Coverage", shrink=0.25, pad=0.02) + cbar.set_ticklabels(["Missing", "Present"]) + + ax.set_title( + f"Grid Coverage: NSB={nsb}, Array Name={array_name}, Primary Particle={primary_particle}" + ) + ax.set_xlabel("Azimuth [deg]") + ax.set_ylabel(ZENITH_LABEL) + ax.set_xticks(az_vals) + ax.set_yticks(zen_vals) + ax.grid(which="major", linestyle="-", linewidth="0.5", color="black", alpha=0.3) + + +def plot_grid_coverage(limits_table, grid_definition, output_dir): + """ + Generate grid coverage plots for each NSB level and array name combination. + + Parameters + ---------- + limits_table : Table + An astropy Table containing the CORSIKA limits data. + grid_definition : dict or None + A dictionary defining the expected grid points for zenith, + azimuth, NSB level, and array name. + output_dir : str or Path + Directory where the generated grid coverage plots will be saved. + + Returns + ------- + list of Path + List of file paths to the generated grid coverage plots. + """ + if not grid_definition: + _logger.info("No grid definition provided, skipping grid coverage plots.") + return [] + + _logger.info("Generating grid coverage plots") + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + output_files = [] + primary_particle = _get_primary_particle_label(limits_table) + + found_combinations_str = set( + zip( + np.array(limits_table["zenith"].value, dtype=str), + np.array(limits_table["azimuth"].value, dtype=str), + np.array(limits_table["nsb_level"], dtype=str), + np.array(limits_table["array_name"], dtype=str), + ) + ) + + unique_values = { + "zeniths": np.array(grid_definition.get("zenith", [])), + "azimuths": np.array(grid_definition.get("azimuth", [])), + "nsb_levels": np.array(grid_definition.get("nsb_level", [])), + "array_names": np.array(grid_definition.get("array_name", [])), + } + + for nsb, array_name in product(unique_values["nsb_levels"], unique_values["array_names"]): + _, ax = plt.subplots(figsize=(10, 8)) + _plot_single_grid_coverage( + ax, + unique_values["zeniths"], + unique_values["azimuths"], + nsb, + array_name, + found_combinations_str, + primary_particle, + ) + output_file = output_dir / f"grid_coverage_{nsb}_{array_name}.png" + plt.tight_layout() + plt.savefig(output_file, bbox_inches="tight") + plt.close() + output_files.append(output_file) + + return output_files + + +def plot_limits(limits_table, output_dir): + """ + Create plots of derived CORSIKA limits for each array name and azimuth. + + Parameters + ---------- + limits_table (Table) + An astropy Table containing the CORSIKA limits data. + output_dir (str or Path) + Directory where the generated plots will be saved. + + Returns + ------- + list of Path + List of file paths to the generated plots. + """ + _logger.info("Generating limit plots") + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + output_files = [] + + grouped_by_layout_az = limits_table.group_by(["array_name", "azimuth"]) + broad_range_columns = _resolve_broad_range_columns(limits_table) + + for group in grouped_by_layout_az.groups: + array_name = group["array_name"][0] + azimuth = group["azimuth"][0] + azimuth_value = azimuth.value if hasattr(azimuth, "value") else azimuth + primary_particle = _get_primary_particle_label(group) + + fig, axes = plt.subplots(1, 3, figsize=(18, 6)) + legend_handles, legend_labels = [], [] + + grouped_by_nsb = group.group_by("nsb_level") + colors = plt.get_cmap("viridis")(np.linspace(0, 1, len(grouped_by_nsb.groups))) + + for i, nsb_group in enumerate(grouped_by_nsb.groups): + nsb_level = nsb_group["nsb_level"][0] + plot_columns = [ + "zenith", + "lower_energy_limit", + "upper_radius_limit", + "viewcone_radius", + ] + agg_data = nsb_group[plot_columns].group_by("zenith").groups.aggregate(np.mean) + agg_data.sort("zenith") + zeniths = agg_data["zenith"].value + + (line,) = axes[0].plot(zeniths, agg_data["lower_energy_limit"], "o-", color=colors[i]) + axes[1].plot(zeniths, agg_data["upper_radius_limit"], "o-", color=colors[i]) + axes[2].plot(zeniths, agg_data["viewcone_radius"], "o-", color=colors[i]) + legend_handles.append(line) + legend_labels.append(f"NSB={nsb_level}") + + if broad_range_columns: + broad_columns = [ + "zenith", + broad_range_columns["lower_energy_limit"], + broad_range_columns["upper_radius_limit"], + broad_range_columns["viewcone_radius"], + ] + broad_data = nsb_group[broad_columns].group_by("zenith").groups.aggregate(np.mean) + broad_data.sort("zenith") + + axes[0].plot( + zeniths, + broad_data[broad_range_columns["lower_energy_limit"]], + linestyle="--", + color="gray", + linewidth=1.5, + ) + axes[1].plot( + zeniths, + broad_data[broad_range_columns["upper_radius_limit"]], + linestyle="--", + color="gray", + linewidth=1.5, + ) + axes[2].plot( + zeniths, + broad_data[broad_range_columns["viewcone_radius"]], + linestyle="--", + color="gray", + linewidth=1.5, + ) + + axes[0].set_title("Lower Energy Limit vs Zenith") + axes[0].set_xlabel(ZENITH_LABEL) + axes[0].set_ylabel("Lower Energy Limit [TeV]") + axes[0].grid(True) + axes[1].set_title("Upper Radius Limit vs Zenith") + axes[1].set_xlabel(ZENITH_LABEL) + axes[1].set_ylabel("Upper Radius Limit [m]") + axes[1].grid(True) + axes[2].set_title("Viewcone Radius vs Zenith") + axes[2].set_xlabel(ZENITH_LABEL) + axes[2].set_ylabel("Viewcone Radius [deg]") + axes[2].grid(True) + + fig.legend(legend_handles, legend_labels, loc="lower center", ncol=len(legend_labels)) + plt.suptitle( + "CORSIKA Limits: " + f"Array Name={array_name}, Azimuth={azimuth_value} deg, " + f"Primary Particle={primary_particle}" + ) + plt.tight_layout() + plt.subplots_adjust(bottom=0.15) + + output_file = output_dir / f"limits_{array_name}_azimuth{azimuth_value}.png" + plt.savefig(output_file) + plt.close(fig) + output_files.append(output_file) + + return output_files diff --git a/src/simtools/production_configuration/plot_production_grid.py b/src/simtools/visualization/plot_production_grid.py similarity index 100% rename from src/simtools/production_configuration/plot_production_grid.py rename to src/simtools/visualization/plot_production_grid.py diff --git a/src/simtools/visualization/plot_simtel_event_histograms.py b/src/simtools/visualization/plot_simtel_event_histograms.py index 1a5cc18f5e..64376bc5b8 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): @@ -53,13 +48,25 @@ def _safe_value(limits, key): "energy": {"x": _safe_value(limits, "lower_energy_limit")}, "core_distance": {"x": _safe_value(limits, "upper_radius_limit")}, "angular_distance": {"x": _safe_value(limits, "viewcone_radius")}, - "core_vs_energy": { + "core_distance_vs_energy": { "x": _safe_value(limits, "upper_radius_limit"), "y": _safe_value(limits, "lower_energy_limit"), + "curve": limits.get("core_distance_vs_energy_curve"), + }, + "core_distance_vs_energy_cumulative": { + "x": _safe_value(limits, "upper_radius_limit"), + "y": _safe_value(limits, "lower_energy_limit"), + "curve": limits.get("core_distance_vs_energy_curve"), }, "angular_distance_vs_energy": { "x": _safe_value(limits, "viewcone_radius"), "y": _safe_value(limits, "lower_energy_limit"), + "curve": limits.get("angular_distance_vs_energy_curve"), + }, + "angular_distance_vs_energy_cumulative": { + "x": _safe_value(limits, "viewcone_radius"), + "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")}, } @@ -69,8 +76,8 @@ def _safe_value(limits, key): def _generate_plot_configurations(histograms, limits): """Generate plot configurations for all histogram types.""" hist_1d_params = {"color": "tab:green", "edgecolor": "tab:green", "lw": 1} - hist_2d_params = {"norm": "log", "cmap": "viridis", "show_contour": False} - hist_2d_normalized_params = {"norm": "linear", "cmap": "viridis", "show_contour": True} + hist_2d_params = {"norm": "log", "cmap": "viridis"} + hist_2d_normalized_params = {"norm": "linear", "cmap": "viridis"} plots = {} for name, hist in histograms.items(): if hist["histogram"] is None: @@ -80,7 +87,12 @@ def _generate_plot_configurations(histograms, limits): hist, name=name, plot_params=hist_1d_params, limits=limits ) else: - if "cumulative" in name or "efficiency" in name: + histogram_name = name.lower() + if ( + "cumulative" in histogram_name + or "efficiency" in histogram_name + or histogram_name.endswith("_eff") + ): plot_params = hist_2d_normalized_params else: plot_params = hist_2d_params @@ -143,7 +155,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") @@ -157,14 +169,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): @@ -186,66 +191,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, @@ -306,15 +251,21 @@ def _plot_data(ax, data, bins, plot_type, plot_params, colorbar_label): def _add_lines(ax, lines): """Add reference lines to the plot.""" - if lines.get("x") is not None: - ax.axvline(lines["x"], color="r", linestyle="--", linewidth=0.5) - if lines.get("y") is not None: - ax.axhline(lines["y"], color="r", linestyle="--", linewidth=0.5) if lines.get("r") is not None: ax.add_artist( plt.Circle((0, 0), lines["r"], color="r", fill=False, linestyle="--", linewidth=0.5) ) + for x_value in np.atleast_1d(lines.get("x", [])): + ax.axvline(x_value, color="r", linestyle="--", linewidth=0.5) + + for y_value in np.atleast_1d(lines.get("y", [])): + ax.axhline(y_value, color="r", linestyle="--", linewidth=0.5) + + curve = lines.get("curve") + if curve and curve.get("x") and curve.get("y"): + ax.plot(curve["x"], curve["y"], color="tab:orange", linestyle="-", linewidth=1.0) + def _create_2d_histogram_plot(data, bins, plot_params): """ @@ -327,50 +278,41 @@ def _create_2d_histogram_plot(data, bins, plot_params): bins : tuple of np.ndarray Bin edges for x and y axes plot_params : dict - Plot parameters including norm, cmap, and show_contour + Plot parameters including norm, cmap Returns ------- matplotlib.collections.QuadMesh The created pcolormesh object for colorbar attachment """ + cmap = plt.get_cmap(plot_params.get("cmap", "viridis")).copy() + cmap.set_bad(color="white", alpha=0.0) + if plot_params.get("norm") == "linear": + masked_data = np.ma.masked_equal(data.T, 0) pcm = plt.pcolormesh( bins[0], bins[1], - data.T, + masked_data, vmin=0, vmax=1, - cmap=plot_params.get("cmap", "viridis"), + cmap=cmap, ) - # Add contour line at value=1.0 for normalized histograms - if plot_params.get("show_contour", True): - x_centers = (bins[0][1:] + bins[0][:-1]) / 2 - y_centers = (bins[1][1:] + bins[1][:-1]) / 2 - x_mesh, y_mesh = np.meshgrid(x_centers, y_centers) - plt.contour( - x_mesh, - y_mesh, - data.T, - levels=[0.999999], # very close to 1 for floating point precision - colors=["tab:red"], - linestyles=["--"], - linewidths=[0.5], - ) else: + masked_data = np.ma.masked_less_equal(data.T, 0) # Handle empty or invalid data for logarithmic scaling data_max = data.max() if data_max <= 0: _logger.warning("No positive data found for logarithmic scaling, using linear scale") pcm = plt.pcolormesh( - bins[0], bins[1], data.T, vmin=0, vmax=max(1, data_max), cmap="viridis" + bins[0], bins[1], masked_data, vmin=0, vmax=max(1, data_max), cmap=cmap ) else: # Ensure vmin is less than vmax for LogNorm vmin = max(1, data[data > 0].min()) if np.any(data > 0) else 1 vmax = max(vmin + 1, data_max) pcm = plt.pcolormesh( - bins[0], bins[1], data.T, norm=LogNorm(vmin=vmin, vmax=vmax), cmap="viridis" + bins[0], bins[1], masked_data, norm=LogNorm(vmin=vmin, vmax=vmax), cmap=cmap ) return pcm 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/integration_tests/config/production_derive_corsika_limits_fits.yml b/tests/integration_tests/config/production_derive_corsika_limits_fits.yml index 4e27ffb13e..254145f0f1 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: + - core_distance,0.001,10 + - angular_distance,0.001,10 + energy_threshold_fraction: 0.0001 output_file: corsika_simulation_limits.ecsv output_path: simtools-output plot_histograms: true 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 d8d67ccbde..9e131c04a9 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,12 +11,16 @@ 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: + - core_distance,0.001,10 + - angular_distance,0.001,10 + energy_threshold_fraction: 0.01 model_version: 6.0.2 output_file: corsika_simulation_limits.ecsv output_path: simtools-output plot_histograms: true site: North + differential_loss_bins_per_decade: 5 integration_tests: - output_file: corsika_simulation_limits.ecsv - output_file: corsika_simulation_limits.meta.yml diff --git a/tests/integration_tests/config/production_merge_corsika_limits_check_only.yml b/tests/integration_tests/config/production_merge_corsika_limits_check_only.yml index 7109167bfb..f855872805 100644 --- a/tests/integration_tests/config/production_merge_corsika_limits_check_only.yml +++ b/tests/integration_tests/config/production_merge_corsika_limits_check_only.yml @@ -4,16 +4,8 @@ applications: configuration: merged_table: tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv grid_definition: tests/resources/corsika_simulation_limits/grid_definition.yaml - plot_grid_coverage: true - plot_limits: true output_path: simtools-output - integration_tests: - - output_file: grid_coverage_1.0_alpha.png - - output_file: grid_coverage_5.0_alpha.png - - output_file: limits_1lst_azimuth0.0.png - - output_file: limits_1lst_azimuth180.0.png - - output_file: limits_1mst_azimuth0.0.png - - output_file: limits_1mst_azimuth180.0.png + integration_tests: [] test_name: production_merge_corsika_limits_check_only schema_name: application_workflow.metaschema schema_version: 0.4.0 diff --git a/tests/integration_tests/config/production_merge_corsika_limits_merge_and_check.yml b/tests/integration_tests/config/production_merge_corsika_limits_merge_and_check.yml index 403f2ff0f5..32e7e90d24 100644 --- a/tests/integration_tests/config/production_merge_corsika_limits_merge_and_check.yml +++ b/tests/integration_tests/config/production_merge_corsika_limits_merge_and_check.yml @@ -7,15 +7,9 @@ applications: - tests/resources/corsika_simulation_limits/corsika_simulation_limits_lookup_02.ecsv grid_definition: tests/resources/corsika_simulation_limits/grid_definition.yaml output_file: merged_limits_merge_and_check.ecsv - plot_grid_coverage: true - plot_limits: true output_path: simtools-output integration_tests: - output_file: merged_limits_merge_and_check.ecsv - - output_file: grid_coverage_1.0_alpha.png - - output_file: grid_coverage_5.0_alpha.png - - output_file: limits_alpha_azimuth0.0.png - - output_file: limits_alpha_azimuth180.0.png test_name: production_merge_corsika_limits_merge_and_check schema_name: application_workflow.metaschema schema_version: 0.4.0 diff --git a/tests/integration_tests/config/production_merge_corsika_limits_merge_only.yml b/tests/integration_tests/config/production_merge_corsika_limits_merge_only.yml index a819e804c5..f3467fc5b3 100644 --- a/tests/integration_tests/config/production_merge_corsika_limits_merge_only.yml +++ b/tests/integration_tests/config/production_merge_corsika_limits_merge_only.yml @@ -6,7 +6,6 @@ applications: - tests/resources/corsika_simulation_limits/corsika_simulation_limits_lookup_01.ecsv - tests/resources/corsika_simulation_limits/corsika_simulation_limits_lookup_02.ecsv output_file: merged_limits_merge_only.ecsv - plot_limits: true output_path: simtools-output integration_tests: - output_file: merged_limits_merge_only.ecsv diff --git a/tests/integration_tests/config/production_plot_corsika_limits.yml b/tests/integration_tests/config/production_plot_corsika_limits.yml new file mode 100644 index 0000000000..db81e47b1d --- /dev/null +++ b/tests/integration_tests/config/production_plot_corsika_limits.yml @@ -0,0 +1,14 @@ +--- +applications: +- application: simtools-production-plot-corsika-limits + configuration: + input: tests/resources/corsika_simulation_limits/merged_corsika_limits_for_test.ecsv + output_path: simtools-output + integration_tests: + - output_file: limits_1lst_azimuth0.0.png + - output_file: limits_1lst_azimuth180.0.png + - output_file: limits_1mst_azimuth0.0.png + - output_file: limits_1mst_azimuth180.0.png + test_name: production_plot_corsika_limits +schema_name: application_workflow.metaschema +schema_version: 0.4.0 diff --git a/tests/unit_tests/applications/test_production_plot_corsika_limits.py b/tests/unit_tests/applications/test_production_plot_corsika_limits.py new file mode 100644 index 0000000000..1054d09a32 --- /dev/null +++ b/tests/unit_tests/applications/test_production_plot_corsika_limits.py @@ -0,0 +1,49 @@ +from pathlib import Path +from types import SimpleNamespace +from unittest.mock import MagicMock, patch + +from astropy.table import Column, Table + +import simtools.applications.production_plot_corsika_limits as app + + +def _create_merged_table(): + return Table( + [ + Column(data=[20.0, 40.0], name="zenith"), + Column(data=[0.0, 180.0], name="azimuth"), + Column(data=["dark", "moon"], name="nsb_level"), + Column(data=["alpha", "alpha"], name="array_name"), + Column(data=[0.1, 0.2], name="lower_energy_limit"), + Column(data=[1200.0, 1500.0], name="upper_radius_limit"), + Column(data=[8.0, 10.0], name="viewcone_radius"), + ] + ) + + +def test_main_reads_table_and_plots(tmp_test_directory): + """Test application orchestration from CLI input to plotting output.""" + output_dir = Path(tmp_test_directory) / "plots" + app_context = SimpleNamespace( + args={"input": "merged_limits.ecsv"}, + io_handler=MagicMock(), + ) + app_context.io_handler.get_output_directory.return_value = output_dir + + merged_table = _create_merged_table() + + with ( + patch( + "simtools.applications.production_plot_corsika_limits.build_application", + return_value=app_context, + ), + patch( + "simtools.applications.production_plot_corsika_limits.data_reader.read_table_from_file", + return_value=merged_table, + ) as mock_read_table, + patch("simtools.applications.production_plot_corsika_limits.plot_limits") as mock_limits, + ): + app.main() + + mock_read_table.assert_called_once_with("merged_limits.ecsv") + mock_limits.assert_called_once_with(merged_table, output_dir) 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 cc6c1911f4..e38519ddf5 100644 --- a/tests/unit_tests/production_configuration/test_derive_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_derive_corsika_limits.py @@ -4,7 +4,6 @@ from astropy.table import Table import simtools.production_configuration.derive_corsika_limits as derive_corsika_limits -from simtools.utils.names import get_array_element_name_from_common_identifier # Constants SIM_EVENTS_HISTOGRAMS_PATH = ( @@ -13,11 +12,12 @@ COMPUTE_LOWER_ENERGY_LIMIT_PATH = ( "simtools.production_configuration.derive_corsika_limits.compute_lower_energy_limit" ) -COMPUTE_UPPER_RADIUS_LIMIT_PATH = ( - "simtools.production_configuration.derive_corsika_limits.compute_upper_radius_limit" -) -COMPUTE_VIEWCONE_PATH = "simtools.production_configuration.derive_corsika_limits.compute_viewcone" +COMPUTE_LIMITS_PATH = "simtools.production_configuration.derive_corsika_limits._compute_limits" MOCK_FILE_PATH = "mock_file.fits" +DEFAULT_ALLOWED_LOSSES = { + "core_distance": {"loss_fraction": 0.2, "loss_min_events": 10}, + "angular_distance": {"loss_fraction": 0.2, "loss_min_events": 10}, +} def _pool_result( @@ -45,104 +45,34 @@ def _pool_result( } -@pytest.fixture -def mock_args_dict(): - """Create mock arguments dictionary.""" - return { - "event_data_file": "data_files.hdf5", - "telescope_ids": "telescope_ids.yml", - "loss_fraction": 0.2, - "plot_histograms": False, - "output_file": "test_output.ecsv", - } - - -@pytest.fixture -def mock_results(): - """Create mock results list.""" - return [ - { - "primary_particle": "gamma", - "telescope_ids": [1, 2], - "zenith": 20.0 * u.deg, - "azimuth": 180.0 * u.deg, - "nsb_level": 1.0, - "lower_energy_limit": 0.5 * u.TeV, - "upper_radius_limit": 400.0 * u.m, - "viewcone_radius": 5.0 * u.deg, - "array_name": "LST", - "layout": "LST", - } - ] - - -def test_generate_corsika_limits_grid(mocker, mock_args_dict, tmp_test_directory): - """Test generate_corsika_limits_grid function with single production.""" - # Mock dependencies - mock_collect = mocker.patch("simtools.io.ascii_handler.collect_data_from_file") - mock_collect.side_effect = [ - {"telescope_configs": {"LST": [1, 2], "MST": [3, 4]}}, - ] - - mock_pool = mocker.patch( - "simtools.production_configuration.derive_corsika_limits.process_pool_map_ordered" - ) - mock_pool.return_value = [ - {"primary_particle": "gamma", "array_name": "LST", "telescope_ids": [1, 2]}, - {"primary_particle": "gamma", "array_name": "MST", "telescope_ids": [3, 4]}, - ] - - mock_io = mocker.patch( - "simtools.production_configuration.derive_corsika_limits.io_handler.IOHandler" - ) - mock_io.return_value.get_output_directory.return_value = tmp_test_directory - - mock_write = mocker.patch( - "simtools.production_configuration.derive_corsika_limits.write_results" - ) - - # Run function (single production mode) - derive_corsika_limits.generate_corsika_limits_grid(mock_args_dict) - - # Verify calls - assert mock_collect.call_count == 1 - mock_pool.assert_called_once() - assert mock_write.call_count == 1 - - -def test_generate_corsika_limits_grid_normalizes_telescope_ids(mocker, mock_args_dict): - """Ensure numeric IDs are normalized to array-element names before processing.""" - mock_collect = mocker.patch("simtools.io.ascii_handler.collect_data_from_file") - mock_collect.return_value = {"telescope_configs": {"LST": [1, 2]}} - - mock_pool = mocker.patch( - "simtools.production_configuration.derive_corsika_limits.process_pool_map_ordered" - ) - mock_pool.return_value = [{}] - mocker.patch("simtools.production_configuration.derive_corsika_limits.write_results") - - derive_corsika_limits.generate_corsika_limits_grid(mock_args_dict) - - expected_telescopes = [ - get_array_element_name_from_common_identifier(1), - get_array_element_name_from_common_identifier(2), - ] - job_specs = mock_pool.call_args[0][1] - assert job_specs[0]["telescope_ids"] == expected_telescopes - - def test_process_file_passes_event_data_patterns_through(mocker): """Test _process_file passes glob patterns through to histogram resolution.""" 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_VIEWCONE_PATH, return_value=2.0 * u.deg) + mocker.patch( + COMPUTE_LIMITS_PATH, + return_value={ + "upper_radius_limit": 100.0 * u.m, + "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, + }, + ) - derive_corsika_limits._process_file("input/*.h5", "array_name", [1, 2], 0.2, False) + derive_corsika_limits._process_file( + "input/*.h5", + "array_name", + [1, 2], + DEFAULT_ALLOWED_LOSSES, + plot_histograms=False, + ) mock_histogram_class.assert_called_once_with( - "input/*.h5", array_name="array_name", telescope_list=[1, 2] + "input/*.h5", + "array_name", + [1, 2], + 10, ) @@ -153,7 +83,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, 0.1) # Verify metadata was written mock_dump.assert_called_once() @@ -163,14 +93,28 @@ 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, 0.1) table.info() 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 table.meta["loss_fraction"] == pytest.approx(0.2) + 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_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 table.meta["energy_threshold_fraction"] == pytest.approx(0.1) assert isinstance(table.meta["created"], str) assert "description" in table.meta @@ -293,12 +237,16 @@ def test_compute_limits_lower(): loss_fraction = 0.2 with pytest.raises(ValueError, match="limit_type must be 'lower' or 'upper'"): - derive_corsika_limits._compute_limits(hist, bin_edges, loss_fraction, limit_type="blabla") - - result = derive_corsika_limits._compute_limits( - hist, bin_edges, loss_fraction, limit_type="lower" + derive_corsika_limits._integral_limits(hist, bin_edges, loss_fraction, limit_type="blabla") + + result = derive_corsika_limits._integral_limits( + hist, + bin_edges, + loss_fraction, + loss_min_events=0, + limit_type="lower", ) - assert result == 2 + assert result == 3 def test_compute_limits_upper(): @@ -306,8 +254,12 @@ def test_compute_limits_upper(): bin_edges = np.array([0, 1, 2, 3, 4, 5]) loss_fraction = 0.2 - result = derive_corsika_limits._compute_limits( - hist, bin_edges, loss_fraction, limit_type="upper" + result = derive_corsika_limits._integral_limits( + hist, + bin_edges, + loss_fraction, + loss_min_events=0, + limit_type="upper", ) assert result == 3 @@ -317,35 +269,48 @@ def test_compute_limits_default_type(): bin_edges = np.array([0, 1, 2, 3, 4, 5]) loss_fraction = 0.2 - result = derive_corsika_limits._compute_limits(hist, bin_edges, loss_fraction) - assert result == 2 + result = derive_corsika_limits._integral_limits( + hist, + bin_edges, + loss_fraction, + loss_min_events=0, + ) + assert result == 3 -def test_compute_viewcone(mocker): - """Test compute_viewcone function with mocked histograms.""" - mock_hist = np.array([10, 8, 6, 4, 2]) - mock_bins = np.linspace(0, 20.0, 6) +def test_compute_limits_enforces_minimum_lost_events_upper(): + """Test _compute_limits enforces an absolute minimum number of lost events.""" + hist = np.array([5, 4, 3, 2, 1]) + bin_edges = np.array([0, 1, 2, 3, 4, 5]) - # Mock the histograms object - mock_histograms = mocker.MagicMock() - mock_histograms.histograms = {"angular_distance": {"histogram": mock_hist}} - mock_histograms.view_cone_bins = mock_bins + result = derive_corsika_limits._integral_limits( + hist, + bin_edges, + loss_fraction=0.2, + loss_min_events=10, + limit_type="upper", + ) + assert result == 1 - result = derive_corsika_limits.compute_viewcone(mock_histograms, 0.2) - assert isinstance(result, u.Quantity) - assert result.unit == u.deg - assert result.value > 0 +def test_compute_limits_enforces_minimum_lost_events_lower(): + """Test lower limits also honor the absolute minimum loss requirement.""" + hist = np.array([1, 2, 3, 4, 5]) + bin_edges = np.array([0, 1, 2, 3, 4, 5]) - expected = ( - derive_corsika_limits._compute_limits(mock_hist, mock_bins, 0.2, limit_type="upper") * u.deg + result = derive_corsika_limits._integral_limits( + hist, + bin_edges, + loss_fraction=0.2, + loss_min_events=10, + limit_type="lower", ) - assert result.value == pytest.approx(expected.value) + assert result == 5 def test_compute_lower_energy_limit(mocker): """Test compute_lower_energy_limit function with mocked histograms.""" - mock_hist = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + mock_hist = np.array([1.0, 12.0, 20.0, 12.0, 1.0]) mock_bins = np.logspace(-3, 3, 6) # Mock the histograms object @@ -361,32 +326,55 @@ def test_compute_lower_energy_limit(mocker): assert result.value > 0 expected = ( - derive_corsika_limits._compute_limits(mock_hist, mock_bins, 0.2, limit_type="lower") * u.TeV + derive_corsika_limits._find_low_energy_threshold_from_histogram( + mock_hist, + mock_bins, + threshold_fraction=0.2, + ) + * u.TeV ) assert result == expected -def test_compute_upper_radius_limit(mocker): - """Test compute_upper_radius_limit function with mocked histograms.""" - mock_hist = np.array([10.0, 20.0, 30.0, 40.0, 50.0]) - mock_bins = np.linspace(0, 500, 6) +def test_find_low_energy_threshold_from_histogram_basic(): + """Test nominal threshold finding from peak toward lower energies.""" + counts = np.array([1.0, 12.0, 20.0, 12.0, 1.0]) + bin_edges = np.array([0.1, 0.2, 0.4, 0.8, 1.6, 3.2]) - # Mock the histograms object - mock_histograms = mocker.MagicMock() - mock_histograms.histograms = {"core_distance": {"histogram": mock_hist}} - mock_histograms.core_distance_bins = mock_bins - mock_histograms.file_info = {} + # Peak index=2, N_peak=(12+20+12)/3=14.666..., threshold=1.466... + # Walking left from idx=2: 20,12,1 -> first below threshold at idx=0 + result = derive_corsika_limits._find_low_energy_threshold_from_histogram(counts, bin_edges) + assert result == pytest.approx(0.1) - 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 +def test_find_low_energy_threshold_from_histogram_peak_at_first_bin(): + """Test edge case where absolute maximum is at the first bin.""" + counts = np.array([10.0, 4.0, 1.0, 0.0]) + bin_edges = np.array([0.05, 0.1, 0.2, 0.4, 0.8]) - expected = ( - derive_corsika_limits._compute_limits(mock_hist, mock_bins, 0.2, limit_type="upper") * u.m - ) - assert result == expected + # No bins left of peak; fallback to first edge is expected. + result = derive_corsika_limits._find_low_energy_threshold_from_histogram(counts, bin_edges) + assert result == pytest.approx(0.05) + + +def test_find_low_energy_threshold_from_histogram_peak_at_last_bin(): + """Test edge case where absolute maximum is at the last bin.""" + counts = np.array([0.0, 0.2, 0.5, 10.0]) + bin_edges = np.array([0.1, 0.2, 0.4, 0.8, 1.6]) + + # Peak index=3, N_peak=(0.5+10)/2=5.25, threshold=0.525 + # Walking left from idx=3: 10,0.5 -> first below threshold at idx=2 + result = derive_corsika_limits._find_low_energy_threshold_from_histogram(counts, bin_edges) + assert result == pytest.approx(0.4) + + +def test_find_low_energy_threshold_from_histogram_raises_for_all_zero_counts(): + """Reject histograms without positive entries.""" + counts = np.array([0.0, 0.0, 0.0, 0.0]) + bin_edges = np.array([0.1, 0.2, 0.4, 0.8, 1.6]) + + with pytest.raises(ValueError, match="at least one positive entry"): + derive_corsika_limits._find_low_energy_threshold_from_histogram(counts, bin_edges) def test_is_close(caplog): @@ -420,20 +408,21 @@ def test_process_file_with_mocked_histograms(mocker): COMPUTE_LOWER_ENERGY_LIMIT_PATH, return_value=1.0 * u.TeV, ) - mock_compute_upper_radius_limit = mocker.patch( - COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, - ) - mock_compute_viewcone = mocker.patch( - COMPUTE_VIEWCONE_PATH, - return_value=2.0 * u.deg, + mock_compute_limits = mocker.patch( + COMPUTE_LIMITS_PATH, + return_value={ + "upper_radius_limit": 100.0 * u.m, + "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, + }, ) result = derive_corsika_limits._process_file( file_path=MOCK_FILE_PATH, array_name="MockArray", telescope_ids=[1, 2], - loss_fraction=0.2, + allowed_losses=DEFAULT_ALLOWED_LOSSES, plot_histograms=False, ) @@ -442,18 +431,249 @@ def test_process_file_with_mocked_histograms(mocker): "zenith": None, "azimuth": None, "nsb_level": None, + "br_energy_min": None, + "br_energy_max": None, + "br_core_scatter_max": None, + "br_viewcone_max": None, "lower_energy_limit": 1.0 * u.TeV, "upper_radius_limit": 100.0 * u.m, "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, } mock_histogram_class.assert_called_once_with( - MOCK_FILE_PATH, array_name="MockArray", telescope_list=[1, 2] + MOCK_FILE_PATH, + "MockArray", + [1, 2], + 10, ) mock_histograms.fill.assert_called_once() - mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.2) - mock_compute_upper_radius_limit.assert_called_once_with(mock_histograms, 0.2) - mock_compute_viewcone.assert_called_once_with(mock_histograms, 0.2) + mock_compute_lower_energy_limit.assert_called_once_with(mock_histograms, 0.01) + mock_compute_limits.assert_called_once_with(mock_histograms, DEFAULT_ALLOWED_LOSSES, 0) + + +def test_process_file_with_differential_loss_per_energy_bin(mocker): + """Test _process_file in differential-loss mode.""" + mock_histograms = mocker.MagicMock() + mock_histograms.fill.return_value = None + mock_histograms.file_info = {} + + mocker.patch( + SIM_EVENTS_HISTOGRAMS_PATH, + return_value=mock_histograms, + ) + + mock_compute_lower_energy_limit = mocker.patch( + COMPUTE_LOWER_ENERGY_LIMIT_PATH, + return_value=1.0 * u.TeV, + ) + mock_compute_limits = mocker.patch(COMPUTE_LIMITS_PATH) + mock_differential = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._compute_limits", + return_value={ + "upper_radius_limit": 120.0 * u.m, + "viewcone_radius": 3.0 * u.deg, + "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]}, + }, + ) + + result = derive_corsika_limits._process_file( + file_path=MOCK_FILE_PATH, + array_name="MockArray", + telescope_ids=[1, 2], + allowed_losses=DEFAULT_ALLOWED_LOSSES, + plot_histograms=False, + differential_loss_bins_per_decade=6, + ) + + assert result["lower_energy_limit"].value == pytest.approx(1.0) + assert result["upper_radius_limit"].value == pytest.approx(120.0) + assert result["viewcone_radius"].value == pytest.approx(3.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.01) + mock_compute_limits.assert_not_called() + mock_differential.assert_called_once_with(mock_histograms, DEFAULT_ALLOWED_LOSSES, 6) + + +@pytest.mark.parametrize( + ("file_info", "expected_core_scatter_max", "expected_viewcone_max"), + [ + ( + {"core_scatter_max": 120.0 * u.m, "viewcone_max": 3.0 * u.deg}, + 120.0 * u.m, + 3.0 * u.deg, + ), + ({}, None, None), + ], +) +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_distance_vs_energy": {"histogram": "core-hist"}, + "angular_distance_vs_energy": {"histogram": "viewcone-hist"}, + } + histograms.file_info = file_info + + mock_diff_limits = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._differential_upper_limits", + side_effect=[ + (120.0, [110.0, 120.0], [1.0, 10.0]), + (3.0, [2.5, 3.0], [1.0, 10.0]), + ], + ) + mock_is_close = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._is_close", + side_effect=[125.0 * u.m, 3.25 * u.deg], + ) + + 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) + 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[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[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 + assert mock_is_close.call_args_list[1].args[0].value == pytest.approx(3.0) + 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() + mock_histograms.file_info = {} + mock_event_histograms = mocker.patch( + "simtools.production_configuration.derive_corsika_limits.EventDataHistograms", + return_value=mock_histograms, + ) + mocker.patch(COMPUTE_LOWER_ENERGY_LIMIT_PATH, return_value=1.0 * u.TeV) + mocker.patch( + "simtools.production_configuration.derive_corsika_limits._compute_limits", + return_value={ + "upper_radius_limit": 120.0 * u.m, + "viewcone_radius": 3.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0], "y": [1.0]}, + "angular_distance_vs_energy_curve": {"x": [3.0], "y": [1.0]}, + }, + ) + + derive_corsika_limits._process_file( + file_path=MOCK_FILE_PATH, + array_name="MockArray", + telescope_ids=[1, 2], + allowed_losses=DEFAULT_ALLOWED_LOSSES, + plot_histograms=False, + differential_loss_bins_per_decade=6, + ) + + mock_event_histograms.assert_called_once_with( + MOCK_FILE_PATH, + "MockArray", + [1, 2], + 6, + ) + + +def test_differential_upper_limits(mocker): + """Test _differential_upper_limits slices energies and skips empty bins.""" + mock_integral_limits = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._integral_limits", + side_effect=[1.5, 2.5], + ) + mock_log = mocker.patch("simtools.production_configuration.derive_corsika_limits._logger.info") + + max_limit, limits, energy_centers = derive_corsika_limits._differential_upper_limits( + histogram2d=np.array([[1.0, 10.0], [2.0, 20.0], [3.0, 30.0]]), + 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]), + allowed_loss=DEFAULT_ALLOWED_LOSSES["core_distance"], + name="core_scatter", + unit="m", + ) + + np.testing.assert_array_equal( + mock_integral_limits.call_args_list[0].args[0], np.array([1.0, 2.0, 3.0]) + ) + np.testing.assert_array_equal( + mock_integral_limits.call_args_list[1].args[0], np.array([10.0, 20.0, 30.0]) + ) + assert max_limit == pytest.approx(2.5) + assert limits == [1.5, 2.5] + assert energy_centers == pytest.approx([np.sqrt(2.0), np.sqrt(7.5)]) + assert mock_log.call_count == 2 + + +def test_differential_upper_limits_falls_back_to_last_bin_edge(mocker): + """Test _differential_upper_limits falls back when all slices are empty.""" + mock_integral_limits = mocker.patch( + "simtools.production_configuration.derive_corsika_limits._integral_limits" + ) + mock_log = mocker.patch("simtools.production_configuration.derive_corsika_limits._logger.info") + + result = derive_corsika_limits._differential_upper_limits( + histogram2d=np.zeros((3, 2)), + 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]), + allowed_loss=DEFAULT_ALLOWED_LOSSES["angular_distance"], + name="viewcone", + unit="deg", + ) + + assert result == (3.0, [], []) + mock_integral_limits.assert_not_called() + mock_log.assert_not_called() def test_process_file_with_plot_histograms(mocker, tmp_test_directory): @@ -477,12 +697,13 @@ def test_process_file_with_plot_histograms(mocker, tmp_test_directory): return_value=1.0 * u.TeV, ) mocker.patch( - COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, - ) - mocker.patch( - COMPUTE_VIEWCONE_PATH, - return_value=2.0 * u.deg, + COMPUTE_LIMITS_PATH, + return_value={ + "upper_radius_limit": 100.0 * u.m, + "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, + }, ) mock_plot = mocker.patch( @@ -493,7 +714,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, ) @@ -507,9 +728,15 @@ def test_process_file_with_plot_histograms(mocker, tmp_test_directory): "zenith": None, "azimuth": None, "nsb_level": None, + "br_energy_min": None, + "br_energy_max": None, + "br_core_scatter_max": None, + "br_viewcone_max": None, "lower_energy_limit": 1.0 * u.TeV, "upper_radius_limit": 100.0 * u.m, "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, } assert kwargs["array_name"] == "MockArray" @@ -572,6 +799,82 @@ 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( + [ + "core_distance,2e-6,20", + "angular_distance,3e-6,30", + ] + ) + + 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", + "core_distance,5e-7,5", + ] + ) + + assert result["core_distance"]["loss_fraction"] == pytest.approx(5e-7) + assert result["core_distance"]["loss_min_events"] == 5 + 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( + [ + "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( + [ + "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_execute_production_job_single_job(mocker): """Test _execute_production_job executes one job correctly.""" mock_histograms = mocker.MagicMock() @@ -588,12 +891,13 @@ def test_execute_production_job_single_job(mocker): return_value=1.0 * u.TeV, ) mocker.patch( - COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, - ) - mocker.patch( - COMPUTE_VIEWCONE_PATH, - return_value=2.0 * u.deg, + COMPUTE_LIMITS_PATH, + return_value={ + "upper_radius_limit": 100.0 * u.m, + "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, + }, ) job_spec = { @@ -601,7 +905,8 @@ 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, + "energy_threshold_fraction": 0.1, "plot_histograms": False, "output_subdir": None, } @@ -655,7 +960,10 @@ 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": [ + "core_distance,0.2,10", + "angular_distance,0.2,10", + ], "plot_histograms": False, "output_file": "test_output.ecsv", "n_workers": 2, @@ -704,7 +1012,10 @@ 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": [ + "core_distance,0.2,10", + "angular_distance,0.2,10", + ], "plot_histograms": False, "output_file": "test_output.ecsv", "n_workers": 0, @@ -725,7 +1036,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, 0.1) # Should include production-origin columns assert "production_index" in table.colnames @@ -739,7 +1050,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, 0.1) # Production-origin columns are included and filled with None if missing assert "production_index" in table.colnames @@ -768,12 +1079,13 @@ def test_process_file_with_output_subdir(mocker, tmp_test_directory): return_value=1.0 * u.TeV, ) mocker.patch( - COMPUTE_UPPER_RADIUS_LIMIT_PATH, - return_value=100.0 * u.m, - ) - mocker.patch( - COMPUTE_VIEWCONE_PATH, - return_value=2.0 * u.deg, + COMPUTE_LIMITS_PATH, + return_value={ + "upper_radius_limit": 100.0 * u.m, + "viewcone_radius": 2.0 * u.deg, + "core_distance_vs_energy_curve": {"x": [100.0, 100.0], "y": [0.1, 1.0]}, + "angular_distance_vs_energy_curve": {"x": [2.0, 2.0], "y": [0.1, 1.0]}, + }, ) mock_plot = mocker.patch( @@ -786,7 +1098,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, ) @@ -795,3 +1107,46 @@ def test_process_file_with_output_subdir(mocker, tmp_test_directory): mock_plot.assert_called_once() call_kwargs = mock_plot.call_args[1] assert call_kwargs["output_path"] == output_subdir + + +@pytest.fixture +def mock_args_dict(): + """Fixture to provide mock arguments dictionary with required keys.""" + return { + "config_file": "dummy_config.yml", + "steps": None, + "ignore_runtime_environment": False, + "event_data_file": "dummy_event_data.h5", + "output_file": "corsika_limits.ecsv", + "allowed_losses": [ + "core_distance,0.2,10", + "angular_distance,0.2,10", + ], + "energy_threshold_fraction": 0.1, + "plot_histograms": False, + "n_workers": 1, + "array_layout_name": None, + "array_element_list": ["LSTN-01"], + "telescope_ids": None, + } + + +@pytest.fixture +def mock_results(): + """Fixture to provide one standard result row for table/writer tests.""" + return [ + { + "primary_particle": "gamma", + "array_name": "LST", + "zenith": 20.0 * u.deg, + "azimuth": 180.0 * u.deg, + "nsb_level": 1.0, + "lower_energy_limit": 0.5 * u.TeV, + "upper_radius_limit": 400.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/production_configuration/test_merge_corsika_limits.py b/tests/unit_tests/production_configuration/test_merge_corsika_limits.py index 082ae06b36..b35f034a8a 100644 --- a/tests/unit_tests/production_configuration/test_merge_corsika_limits.py +++ b/tests/unit_tests/production_configuration/test_merge_corsika_limits.py @@ -169,52 +169,6 @@ def test_check_grid_completeness(tmp_test_directory): assert not result -@patch("matplotlib.pyplot.savefig") -def test_plot_grid_coverage(mock_savefig, tmp_test_directory): - """Test generating grid coverage plots.""" - table = vstack( - [ - create_test_table(20, 0, "dark", "layout1"), - create_test_table(40, 0, "dark", "layout1"), - ] - ) - grid_definition = { - "zenith": [20, 40], - "azimuth": [0], - "nsb_level": ["dark"], - "array_name": ["layout1"], - } - merger = CorsikaMergeLimits(output_dir=tmp_test_directory) - - # Test plotting with grid definition - output_files = merger.plot_grid_coverage(table, grid_definition) - assert len(output_files) == 1 - mock_savefig.assert_called_once() - - # Test plotting without grid definition - mock_savefig.reset_mock() - output_files = merger.plot_grid_coverage(table, None) - assert not output_files - mock_savefig.assert_not_called() - - -@patch("matplotlib.pyplot.savefig") -def test_plot_limits(mock_savefig, tmp_test_directory): - """Test generating limit plots.""" - table = vstack( - [ - create_test_table(20, 0, "dark", "layout1"), - create_test_table(40, 0, "dark", "layout1"), - create_test_table(20, 0, "moon", "layout1"), - ] - ) - merger = CorsikaMergeLimits(output_dir=tmp_test_directory) - output_files = merger.plot_limits(table) - - assert len(output_files) == 1 - mock_savefig.assert_called_once() - - def test_write_merged_table(tmp_test_directory): """Test writing the merged table to file.""" table = create_test_table(20, 0, "dark", "layout1") @@ -341,12 +295,6 @@ def check_grid_completeness(self, merged_table, grid_definition): self.calls.append(("check_grid_completeness", merged_table, grid_definition)) return True, {"expected": 1, "found": 1, "missing": []} - def plot_grid_coverage(self, merged_table, grid_definition): - self.calls.append(("plot_grid_coverage", merged_table, grid_definition)) - - def plot_limits(self, merged_table): - self.calls.append(("plot_limits", merged_table)) - def write_merged_table(self, merged_table, output_file, input_files, grid_completeness): self.calls.append( ("write_merged_table", merged_table, output_file, input_files, grid_completeness) @@ -358,8 +306,6 @@ def write_merged_table(self, merged_table, output_file, input_files, grid_comple "input_files_list": None, "merged_table": None, "grid_definition": None, - "plot_grid_coverage": True, - "plot_limits": True, "output_file": "merged_output.ecsv", } @@ -368,8 +314,6 @@ def write_merged_table(self, merged_table, output_file, input_files, grid_comple call_names = [call[0] for call in merger.calls] assert "merge_tables" in call_names assert "check_grid_completeness" in call_names - assert "plot_grid_coverage" in call_names - assert "plot_limits" in call_names assert "write_merged_table" in call_names @@ -382,12 +326,6 @@ def __init__(self, output_dir): def check_grid_completeness(self, _merged_table, _grid_definition): return True, {} - def plot_grid_coverage(self, _merged_table, _grid_definition): - return None - - def plot_limits(self, _merged_table): - return None - def write_merged_table(self, *_args, **_kwargs): self.write_called = True @@ -403,8 +341,6 @@ def write_merged_table(self, *_args, **_kwargs): { "merged_table": "merged.ecsv", "grid_definition": None, - "plot_grid_coverage": False, - "plot_limits": False, "output_file": "unused.ecsv", }, merger=merger, diff --git a/tests/unit_tests/sim_events/test_histograms.py b/tests/unit_tests/sim_events/test_histograms.py index 0c84b164f0..bf26646979 100644 --- a/tests/unit_tests/sim_events/test_histograms.py +++ b/tests/unit_tests/sim_events/test_histograms.py @@ -80,7 +80,42 @@ def test_energy_bins(mock_reader, hdf5_file_name): mock_reader.return_value.triggered_shower_data.simulated_energy = np.array([1, 10, 100]) bins = histograms.energy_bins assert isinstance(bins, np.ndarray) - assert len(bins) == 100 + assert len(bins) == 61 + assert bins[0] == pytest.approx(1.0e-3) + assert bins[-1] == pytest.approx(1.0e3) + assert np.allclose(np.diff(np.log10(bins)), 0.1) + + +def test_energy_bins_use_file_info_energy_max(mock_reader, hdf5_file_name): + """Test energy_bins uses energy_max from file_info when available.""" + histograms = EventDataHistograms(hdf5_file_name) + histograms.file_info = { + "energy_min": 0.1 * u.TeV, + "energy_max": 30.0 * u.TeV, + } + + bins = histograms.energy_bins + + assert bins[0] == pytest.approx(0.1) + assert bins[-1] == pytest.approx(100.0) + assert len(bins) == 31 + assert np.allclose(np.diff(np.log10(bins)), 0.1) + + +def test_energy_bins_use_configured_bins_per_decade(mock_reader, hdf5_file_name): + """Test energy_bins respects the configured logarithmic resolution.""" + histograms = EventDataHistograms(hdf5_file_name, energy_bins_per_decade=5) + histograms.file_info = { + "energy_min": 0.1 * u.TeV, + "energy_max": 30.0 * u.TeV, + } + + bins = histograms.energy_bins + + assert bins[0] == pytest.approx(0.1) + assert bins[-1] == pytest.approx(100.0) + assert len(bins) == 16 + assert np.allclose(np.diff(np.log10(bins)), 0.2) def test_core_distance_bins(mock_reader, hdf5_file_name): @@ -588,9 +623,10 @@ def test_energy_bins_default(mock_reader, hdf5_file_name): bins = histograms.energy_bins assert isinstance(bins, np.ndarray) - assert len(bins) == 100 + assert len(bins) == 61 assert bins[0] == pytest.approx(1.0e-3) assert bins[-1] == pytest.approx(1.0e3) + assert np.allclose(np.diff(np.log10(bins)), 0.1) def test_core_distance_bins_with_file_info(mock_reader, hdf5_file_name): @@ -667,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"], @@ -679,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( @@ -699,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"], @@ -780,6 +820,27 @@ def test_energy_bins_with_histogram_edges(mock_reader, hdf5_file_name): assert np.array_equal(bins, mock_edges) +def test_update_file_info_stores_energy_max(mock_reader, hdf5_file_name): + """Test _update_file_info keeps energy_max from reduced file info.""" + histograms = EventDataHistograms(hdf5_file_name) + file_info_table = { + "primary_particle": "gamma", + "zenith": 20.0 * u.deg, + "azimuth": 0.0 * u.deg, + "nsb_level": 1.0, + "energy_min": 0.03 * u.TeV, + "energy_max": 30.0 * u.TeV, + "core_scatter_max": 500.0 * u.m, + "viewcone_max": 10.0 * u.deg, + "solid_angle": 1.0 * u.sr, + "scatter_area": 1.0 * u.cm**2, + } + + histograms._update_file_info(file_info_table) + + assert_quantity_allclose(histograms.file_info["energy_max"], 30.0 * u.TeV) + + def test_print_summary(mock_histograms, mocker, caplog): """Test the print_summary method.""" histograms = mock_histograms diff --git a/tests/unit_tests/visualization/test_plot_corsika_limits.py b/tests/unit_tests/visualization/test_plot_corsika_limits.py new file mode 100644 index 0000000000..3b1f466518 --- /dev/null +++ b/tests/unit_tests/visualization/test_plot_corsika_limits.py @@ -0,0 +1,143 @@ +from unittest.mock import patch + +import matplotlib.pyplot as plt +from astropy.table import Column, Table, vstack + +from simtools.visualization import plot_corsika_limits + + +def create_test_table( + zenith, + azimuth, + nsb_level, + array_name="test_array", + lower_energy_limit=0.01, + br_energy_min=None, + br_core_scatter_max=None, + br_viewcone_max=None, +): + """Create a minimal CORSIKA limits table row for plotting tests.""" + columns = [ + Column(data=["gamma"], name="primary_particle"), + Column(data=[array_name], name="array_name"), + Column(data=[[1, 2, 3]], name="telescope_ids"), + Column(data=[zenith], name="zenith"), + Column(data=[azimuth], name="azimuth"), + Column(data=[nsb_level], name="nsb_level"), + Column(data=[lower_energy_limit], name="lower_energy_limit"), + Column(data=[2000], name="upper_radius_limit"), + Column(data=[10], name="viewcone_radius"), + ] + + if br_energy_min is not None: + columns.append(Column(data=[br_energy_min], name="br_energy_min")) + if br_core_scatter_max is not None: + columns.append(Column(data=[br_core_scatter_max], name="br_core_scatter_max")) + if br_viewcone_max is not None: + columns.append(Column(data=[br_viewcone_max], name="br_viewcone_max")) + + return Table(columns) + + +@patch("simtools.visualization.plot_corsika_limits.plt.savefig") +def test_plot_grid_coverage(mock_savefig, tmp_test_directory): + """Test generating grid coverage plots.""" + table = vstack( + [ + create_test_table(20, 0, "dark", "layout1"), + create_test_table(40, 0, "dark", "layout1"), + ] + ) + grid_definition = { + "zenith": [20, 40], + "azimuth": [0], + "nsb_level": ["dark"], + "array_name": ["layout1"], + } + + output_files = plot_corsika_limits.plot_grid_coverage( + table, grid_definition, tmp_test_directory + ) + assert len(output_files) == 1 + mock_savefig.assert_called_once() + + _, ax = plt.subplots() + found_combinations_str = {("20", "0", "dark", "layout1"), ("40", "0", "dark", "layout1")} + plot_corsika_limits._plot_single_grid_coverage( + ax, + [20, 40], + [0], + "dark", + "layout1", + found_combinations_str, + "gamma", + ) + assert "Primary Particle=gamma" in ax.get_title() + plt.close(ax.figure) + + mock_savefig.reset_mock() + output_files = plot_corsika_limits.plot_grid_coverage(table, None, tmp_test_directory) + assert not output_files + mock_savefig.assert_not_called() + + +@patch("simtools.visualization.plot_corsika_limits.plt.savefig") +@patch("simtools.visualization.plot_corsika_limits.plt.suptitle") +def test_plot_limits(mock_suptitle, mock_savefig, tmp_test_directory): + """Test generating CORSIKA limits plots.""" + table = vstack( + [ + create_test_table(20, 0, "dark", "layout1"), + create_test_table(40, 0, "dark", "layout1"), + create_test_table(20, 0, "moon", "layout1"), + ] + ) + + output_files = plot_corsika_limits.plot_limits(table, tmp_test_directory) + assert len(output_files) == 1 + mock_savefig.assert_called_once() + suptitle_text = mock_suptitle.call_args.args[0] + assert "Primary Particle=gamma" in suptitle_text + + +@patch("simtools.visualization.plot_corsika_limits.plt.savefig") +def test_plot_limits_with_broad_range_overlay(mock_savefig, tmp_test_directory): + """Test plotting broad-range limits as gray dashed overlays.""" + table = vstack( + [ + create_test_table( + 20, + 0, + "dark", + "layout1", + lower_energy_limit=0.01, + br_energy_min=0.008, + br_core_scatter_max=2200, + br_viewcone_max=11, + ), + create_test_table( + 40, + 0, + "dark", + "layout1", + lower_energy_limit=0.02, + br_energy_min=0.015, + br_core_scatter_max=2300, + br_viewcone_max=12, + ), + ] + ) + + with patch("simtools.visualization.plot_corsika_limits.plt.close"): + output_files = plot_corsika_limits.plot_limits(table, tmp_test_directory) + + assert len(output_files) == 1 + assert mock_savefig.called + fig = plt.gcf() + axes = fig.axes[:3] + assert any( + line.get_linestyle() == "--" and line.get_color() == "gray" + for axis in axes + for line in axis.get_lines() + ) + plt.close(fig) diff --git a/tests/unit_tests/production_configuration/test_plot_production_grid.py b/tests/unit_tests/visualization/test_plot_production_grid.py similarity index 99% rename from tests/unit_tests/production_configuration/test_plot_production_grid.py rename to tests/unit_tests/visualization/test_plot_production_grid.py index a8a72c4b4a..f89fb8f663 100644 --- a/tests/unit_tests/production_configuration/test_plot_production_grid.py +++ b/tests/unit_tests/visualization/test_plot_production_grid.py @@ -10,7 +10,7 @@ from astropy.time import Time from astropy.utils import iers -from simtools.production_configuration.plot_production_grid import ( +from simtools.visualization.plot_production_grid import ( DEFAULT_OUTPUT_FILE_STEM, ProductionGridPlotter, ) 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 10972b4363..ef44011662 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" @@ -75,6 +72,21 @@ def test_create_2d_histogram_plot_log(sample_data): plt.close(fig) +def test_create_2d_histogram_plot_masks_zero_bins(): + data = np.array([[0, 2], [3, 0]]) + bins = (np.array([0, 1, 2]), np.array([0, 1, 2])) + plot_params = {"norm": "linear", "cmap": "viridis", "show_contour": False} + + fig, _ = plt.subplots() + pcm = _create_2d_histogram_plot(data, bins, plot_params) + + plotted_array = pcm.get_array() + assert np.ma.isMaskedArray(plotted_array) + assert np.any(np.ma.getmaskarray(plotted_array)) + assert pcm.cmap.get_bad()[-1] == pytest.approx(0.0) + plt.close(fig) + + def test_create_2d_histogram_plot_no_positive_data(): data = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) bins_x = np.array([0, 1, 2, 3]) @@ -93,8 +105,9 @@ def test_create_2d_histogram_plot_no_positive_data(): @pytest.mark.parametrize( ("lines", "expect_lines", "expect_circles"), [ - ({"x": 1, "y": 2}, True, 0), + ({"x": 1, "y": 2}, False, 0), ({"r": 3}, False, 1), + ({"curve": {"x": [1, 2], "y": [3, 4]}}, True, 0), ({}, False, 0), ], ) @@ -102,10 +115,10 @@ def test_add_lines(lines, expect_lines, expect_circles): fig, ax = plt.subplots() plot_simtel_event_histograms._add_lines(ax, lines) if expect_lines: - if "x" in lines: - assert any(line.get_xdata() == [lines["x"], lines["x"]] for line in ax.get_lines()) - if "y" in lines: - assert any(line.get_ydata() == [lines["y"], lines["y"]] for line in ax.get_lines()) + if "curve" in lines: + plotted = ax.get_lines()[-1] + np.testing.assert_array_equal(plotted.get_xdata(), np.array([1, 2])) + np.testing.assert_array_equal(plotted.get_ydata(), np.array([3, 4])) else: if not lines or ("x" not in lines and "y" not in lines): remaining = [ln for ln in ax.get_lines() if ln.get_label() == "_nolegend_"] @@ -254,8 +267,6 @@ def test_create_plot(): np.testing.assert_array_equal(bar_args[1], data) np.testing.assert_array_equal(bar_kwargs["width"], np.diff(bins)) assert bar_kwargs["color"] == plot_params["color"] - mock_ax.axvline.assert_called_once_with(1, color="r", linestyle="--", linewidth=0.5) - mock_ax.axhline.assert_called_once_with(2, color="r", linestyle="--", linewidth=0.5) mock_ax.set.assert_called_once_with( xlabel="X-axis", ylabel="Y-axis", @@ -337,7 +348,7 @@ def test_create_plot_histogram2d_colorbar(): assert result is mock_fig # Ensure internal 2d creation helpers used mock_pcolor.assert_called_once() - mock_contour.assert_called_once() + mock_contour.assert_not_called() mock_colorbar.assert_called_once() mock_show.assert_called_once() mock_fig.savefig.assert_not_called() @@ -357,79 +368,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" @@ -464,17 +402,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 @@ -496,65 +432,17 @@ 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 = { "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"}, @@ -570,7 +458,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, ) @@ -583,16 +471,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(): @@ -734,6 +622,29 @@ def test_get_limits(): result = _get_limits("angular_distance", limits) assert result == {"x": 5} + 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_distance_vs_energy", limits) + assert result["x"] == 100 + assert result["y"] == pytest.approx(0.1) + assert result["curve"] == limits["core_distance_vs_energy_curve"] + + 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_distance_vs_energy_curve"] + + result = _get_limits("angular_distance_vs_energy", limits) + assert result["x"] == 5 + assert result["y"] == pytest.approx(0.1) + assert result["curve"] == limits["angular_distance_vs_energy_curve"] + + result = _get_limits("angular_distance_vs_energy_cumulative", limits) + assert result["x"] == 5 + assert result["y"] == pytest.approx(0.1) + assert result["curve"] == limits["angular_distance_vs_energy_curve"] + @pytest.fixture def mock_histograms(): @@ -745,7 +656,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 ( @@ -758,19 +668,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 ( @@ -783,14 +691,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(): @@ -833,3 +738,10 @@ def test_generate_plot_configurations(): _, kwargs = mock_create_2d.call_args assert "plot_params" in kwargs assert kwargs["plot_params"]["norm"] == "linear" + + histos = {"core_distance_vs_energy_eff": {"histogram": "abc", "1d": False}} + with patch(f"{MOD}._create_2d_plot_config") as mock_create_2d: + plot_simtel_event_histograms._generate_plot_configurations(histos, None) + _, kwargs = mock_create_2d.call_args + assert "plot_params" in kwargs + assert kwargs["plot_params"]["norm"] == "linear"