Improve derivation of CORSIKA limits and allow to define a loss fraction per energy bin ('differential loss').

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').

docs/changes/2196.feature.md

@@ -0,0 +1 @@
+Add application to plot CORSIKA limits. Rearrange limits plotting for `production_merge_corsika_limits.py`.

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}

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:
+```

docs/source/user-guide/applications.md

@@ -94,6 +94,7 @@ simtools-production-derive-corsika-limits <applications/simtools-production-deri
 simtools-production-derive-statistics <applications/simtools-production-derive-statistics>
 simtools-production-generate-grid <applications/simtools-production-generate-grid>
 simtools-production-merge-corsika-limits <applications/simtools-production-merge-corsika-limits>
+simtools-production-plot-corsika-limits <applications/simtools-production-plot-corsika-limits>
 simtools-run-application <applications/simtools-run-application>
 simtools-simulate-flasher <applications/simtools-simulate-flasher>
 simtools-simulate-illuminator <applications/simtools-simulate-illuminator>

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:

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"

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

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

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,15 +53,15 @@
       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:
 
    .. code-block:: console
 
       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():