diff --git a/docs/changes/2457.maintenance.md b/docs/changes/2457.maintenance.md new file mode 100644 index 0000000000..93ccd3f7df --- /dev/null +++ b/docs/changes/2457.maintenance.md @@ -0,0 +1 @@ +Exporting trigger threshold as parameter and updating plotting. diff --git a/docs/source/user-guide/applications/simtools-derive-bias-curves.md b/docs/source/user-guide/applications/simtools-derive-bias-curves.md index 4188e141ec..f2e85bff92 100644 --- a/docs/source/user-guide/applications/simtools-derive-bias-curves.md +++ b/docs/source/user-guide/applications/simtools-derive-bias-curves.md @@ -16,37 +16,44 @@ The tool: 1. Extracts NSB trigger rates from reduced event-data HDF5 files 2. Calculates proton trigger rates from proton reduced event-data HDF5 files 3. Plots both curves on the same figure for comparison -4. Outputs ecsv files for runwise nsb simulation, -runwise proton simulation, nsb rate and proton rate vs threshold +4. Outputs ECSV tables for runwise NSB simulation, runwise proton simulation, and combined bias curves +5. Calculates the trigger threshold as the intersection point between the NSB curve and the scaled proton curve +6. Exports the trigger threshold as a model parameter (e.g., ``asum_threshold`` or ``dsum_threshold`` depending on the telescope's default trigger type) The input directory should contain both: -- NSB reduced event-data HDF5 files -- Proton simulation reduced event-data HDF5 files +- NSB reduced event-data HDF5 files (e.g., ``gamma*.reduced_event_data.hdf5``) +- Proton simulation reduced event-data HDF5 files (e.g., ``proton*.reduced_event_data.hdf5``) -The input files can be generated using simtools-generate-bias-curve-submissions. +The input files can be generated using ``simtools-generate-bias-curve-submissions``. **Command line arguments** data_dir (str, required) Directory containing NSB/proton reduced event-data HDF5 files (e.g. gamma* and proton*). +scaling_factor (float, optional) + Scaling factor to account for ions not simulated in the proton dataset. Default: 1.35 figure_file (str, optional) Output plot file path or output directory. Default: bias_curve.png nsb_table_file (str, optional) Output ECSV table file for NSB trigger rates. If not specified, no table is written. proton_table_file (str, optional) Output ECSV table file for proton rates. If not specified, no table is written. +parameter_version (str, required) + Model parameter version for the exported trigger threshold (e.g., ``1.0.0``). +title (str, optional) + Title for the bias curve plot. Default: "Trigger Rate Bias Curves". site (str, required) Site name (North/South) for telescope configuration. model_version (str, required) Model version for telescope configuration. telescope (str, required) Telescope name for configuration. -title (str, optional) - Plot title. Default: "Trigger Rate Bias Curves" -ymin (float, optional) - Minimum y-axis value for plot. Default: 1e2 -ymax (float, optional) - Maximum y-axis value for plot. Default: 5e5 + +**Notes** + +- The trigger threshold is calculated as the intersection between the NSB trigger rate curve and the scaled proton trigger rate curve (scaled by ``scaling_factor``). +- The exported model parameter (``asum_threshold`` or ``dsum_threshold``) is written to the standard model data output directory under ``//``. +- If no intersection point is found, the application raises an error. **Example** @@ -57,6 +64,8 @@ ymax (float, optional) --site North \\ --model_version 7.0.0 \\ --telescope LSTN-01 \\ + --parameter_version 1.0.0 \\ + --scaling_factor 1.35 \\ --figure_file bias_curves.png ``` diff --git a/src/simtools/applications/derive_bias_curves.py b/src/simtools/applications/derive_bias_curves.py index fab1594ea6..66b8d0e410 100644 --- a/src/simtools/applications/derive_bias_curves.py +++ b/src/simtools/applications/derive_bias_curves.py @@ -18,6 +18,13 @@ "(e.g. gamma* and proton*)." ), ), + cli.ArgumentDefinition( + "scaling_factor", + type=float, + required=True, + help=("Scaling factor to account for ions we didn't simulate"), + default=1.35, + ), cli.ArgumentDefinition( "figure_file", type=Path, @@ -37,19 +44,10 @@ help="Output ECSV table file for proton rates. If not specified, no table is written.", ), cli.ArgumentDefinition( - "title", type=str, default="Trigger Rate Bias Curves", help="Plot title." - ), - cli.ArgumentDefinition( - "ymin", - type=float, - default=100.0, - help="Minimum trigger rate value for plotting. Default: 1e2", - ), - cli.ArgumentDefinition( - "ymax", - type=float, - default=500000.0, - help="Maximum trigger rate value for plotting. Default: 5e5", + "title", + type=str, + default="Trigger Rate Bias Curves", + help="Title for the bias curve plot. Default: 'Trigger Rate Bias Curves'", ), ) @@ -59,6 +57,7 @@ arguments=( *_ARGUMENTS, cli.MODEL_VERSION, + cli.PARAMETER_VERSION(required=True), cli.OVERWRITE_MODEL_PARAMETERS, cli.SITE, cli.TELESCOPE, diff --git a/src/simtools/simtel/bias_curve_generator.py b/src/simtools/simtel/bias_curve_generator.py index d03d0abc72..7ae0846011 100644 --- a/src/simtools/simtel/bias_curve_generator.py +++ b/src/simtools/simtel/bias_curve_generator.py @@ -7,6 +7,8 @@ from astropy import units as u from astropy.table import Table +from simtools.data_model import model_data_writer +from simtools.io import io_handler from simtools.model.telescope_model import TelescopeModel from simtools.simtel.nsb_trigger_calculator import ( derive_nsb_triggers, @@ -42,11 +44,25 @@ def generate_bias_curves(args): plot_output_path = plot_tables.resolve_plot_output_path(args["figure_file"]) bias_curve_table_output = plot_output_path.with_suffix(".ecsv") + trigger_threshold = _calculate_trigger_threshold(args, nsb_stats, proton_stats) - _logger.info("Plotting bias curves...") - plot_tables.plot_bias_curves(nsb_stats, proton_stats, args, plot_output_path) + # Log the data points + _logger.info("Trigger threshold calculation data:") + thresholds = sorted(set(nsb_stats.keys()) | set(proton_stats.keys())) + scaling_factor = args.get("scaling_factor", 1.35) + for thresh in thresholds: + nsb_rate = nsb_stats[thresh]["rate_hz"] if thresh in nsb_stats else None + proton_rate = proton_stats[thresh]["rate_hz"] if thresh in proton_stats else None + scaled_proton = scaling_factor * proton_rate if proton_rate is not None else None + _logger.info( + f" Threshold {thresh}: NSB={nsb_rate} Hz, " + f"Proton={proton_rate} Hz, Scaled={scaled_proton} Hz" + ) + _logger.info("Plotting bias curves...") + plot_tables.plot_bias_curves(nsb_stats, proton_stats, args, plot_output_path, trigger_threshold) _write_bias_curve_ecsv(nsb_stats, proton_stats, bias_curve_table_output) + _export_trigger_threshold_as_model_parameter(args, trigger_threshold) _logger.info(f"Bias curve plot written to {plot_output_path}") _logger.info(f"Bias curve table written to {bias_curve_table_output}") @@ -335,3 +351,182 @@ def _write_bias_curve_ecsv(nsb_stats, proton_stats, output_file): output_file.parent.mkdir(parents=True, exist_ok=True) table.write(output_file, format="ascii.ecsv", overwrite=True) + + +def _calculate_trigger_threshold(args, nsb_stats, proton_stats): + """ + Calculate trigger threshold from bias curve intersection. + + Trigger threshold is calculated as the intersection between NSB curve and + scaled proton curve (using the scaling factor from args). + + Parameters + ---------- + args : dict + Dictionary with configuration parameters. + nsb_stats : dict + NSB statistics by threshold. + proton_stats : dict + Proton statistics by threshold. + + Returns + ------- + float + The calculated trigger threshold. + + Raises + ------ + ValueError + If no valid threshold points exist or intersection cannot be found. + """ + # Get all unique thresholds from both NSB and proton stats + thresholds = sorted(set(nsb_stats.keys()) | set(proton_stats.keys())) + # Extract rates for each threshold + nsb_rates = [] + proton_rates = [] + for threshold in thresholds: + nsb_rate = nsb_stats[threshold]["rate_hz"] if threshold in nsb_stats else np.nan + proton_rate = proton_stats[threshold]["rate_hz"] if threshold in proton_stats else np.nan + nsb_rates.append(nsb_rate) + proton_rates.append(proton_rate) + nsb_rates = np.array(nsb_rates) + proton_rates = np.array(proton_rates) + thresholds = np.array(thresholds) + + # Remove NaN values (keep only thresholds where both NSB and proton data exist) + valid_mask = ~(np.isnan(nsb_rates) | np.isnan(proton_rates)) + nsb_rates = nsb_rates[valid_mask] + proton_rates = proton_rates[valid_mask] + thresholds = thresholds[valid_mask] + + if len(thresholds) == 0: + raise ValueError( + "No valid threshold points with both NSB and proton data. " + "Cannot calculate trigger threshold." + ) + # Scale proton rates to account for ions we didn't simulate + scaling_factor = args.get("scaling_factor", 1.35) + scaled_proton_rates = scaling_factor * proton_rates + trigger_threshold = _find_intersection_point(thresholds, nsb_rates, scaled_proton_rates) + if trigger_threshold is not None: + _logger.info(f"Calculated trigger threshold: {trigger_threshold}") + return trigger_threshold + raise ValueError("Could not find intersection point between NSB and scaled proton curves.") + + +def _find_intersection_point(thresholds, nsb_rates, scaled_proton_rates): + """ + Find the threshold value where NSB trigger rate intersects with scaled proton trigger rate. + + Uses linear interpolation between the two data points that bracket the intersection. + + Parameters + ---------- + thresholds : numpy.ndarray + Threshold values from bias curve. + nsb_rates : numpy.ndarray + NSB trigger rates at each threshold. + scaled_proton_rates : numpy.ndarray + Scaled proton trigger rates at each threshold. + + Returns + ------- + float or None + Threshold value at intersection point, or None if no intersection found. + """ + # Sort by threshold to ensure ordering + sort_idx = np.argsort(thresholds) + x = thresholds[sort_idx] + y_nsb = nsb_rates[sort_idx] + y_proton = scaled_proton_rates[sort_idx] + + # Find where NSB crosses below scaled proton + # Look for: y_nsb[i] > y_proton[i] and y_nsb[i+1] <= y_proton[i+1] + for i in range(len(x) - 1): + if y_nsb[i] > y_proton[i] and y_nsb[i + 1] <= y_proton[i + 1]: + # Found the bracket + x1, x2 = float(x[i]), float(x[i + 1]) + y1_nsb, y2_nsb = float(y_nsb[i]), float(y_nsb[i + 1]) + y1_proton, y2_proton = float(y_proton[i]), float(y_proton[i + 1]) + + # Linear interpolation + # At intersection: y1_nsb + t*(y2_nsb-y1_nsb) = y1_proton + t*(y2_proton-y1_proton) + # Solve for t: t = (y1_proton - y1_nsb) / ((y2_nsb - y1_nsb) - (y2_proton - y1_proton)) + numerator = y1_proton - y1_nsb + denominator = (y2_nsb - y1_nsb) - (y2_proton - y1_proton) + t = numerator / denominator + # Clamp t to [0, 1] to stay within bracket + t = max(0.0, min(1.0, t)) + return x1 + t * (x2 - x1) + + _logger.debug("No intersection found") + return None + + +def _export_trigger_threshold_as_model_parameter(args, trigger_threshold): + """ + Export trigger threshold as a model parameter. + + Determines whether to use asum_threshold or dsum_threshold based on the + telescope's default_trigger parameter. + + Parameters + ---------- + args : dict + Dictionary with configuration parameters. + trigger_threshold : float + The calculated trigger threshold value. + """ + try: + # Get telescope name from args + telescope_name = args.get("telescope") + if not telescope_name: + _logger.warning("No telescope name provided. Using 'unknown' as telescope name.") + telescope_name = "unknown" + parameter_version = args["parameter_version"] + + # Determine which threshold parameter to use based on default_trigger + telescope_model = TelescopeModel( + site=args["site"], + telescope_name=telescope_name, + model_version=args["model_version"], + ) + default_trigger = telescope_model.get_parameter_value("default_trigger") + + if default_trigger == "AnalogSum": + parameter_name = "asum_threshold" + # asum_threshold expects float64 in mV + value = round(trigger_threshold, 2) + unit = "mV" + elif default_trigger == "DigitalSum": + parameter_name = "dsum_threshold" + # dsum_threshold expects int64 in counts + value = round(trigger_threshold) + unit = "count" + else: + _logger.warning( + f"Unknown default_trigger '{default_trigger}' for telescope {telescope_name}. " + "Cannot export trigger threshold." + ) + return + + output_path = io_handler.IOHandler().get_output_directory() + output_file = f"{parameter_name}-{parameter_version}.json" + model_data_writer.ModelDataWriter.write_model_parameter( + parameter_name=parameter_name, + value=value, + instrument=telescope_name, + parameter_version=parameter_version, + output_file=output_file, + output_path=output_path / telescope_name / parameter_name, + metadata_input_dict={**args, "source": "bias_curve_analysis"}, + unit=unit, + check_db_for_existing_parameter=False, + ) + + _logger.info( + f"Exported trigger threshold as {parameter_name} for {telescope_name}: {value}" + ) + + except (OSError, ValueError, KeyError) as exc: + _logger.warning(f"Failed to export trigger threshold as model parameter: {exc}") diff --git a/src/simtools/simtel/nsb_trigger_calculator.py b/src/simtools/simtel/nsb_trigger_calculator.py index 7c36f16b6d..a5338aedbd 100644 --- a/src/simtools/simtel/nsb_trigger_calculator.py +++ b/src/simtools/simtel/nsb_trigger_calculator.py @@ -18,8 +18,8 @@ def extract_run_number(file_path): """Extract run number from the FILE_INFO table in the HDF5 file.""" try: - file_info = Table.read(file_path, path="FILE_INFO") - except OSError, ValueError, KeyError: + file_info = Table.read(file_path, path="FILE_INFO", format="hdf5") + except OSError, ValueError, KeyError, ImportError: _logger.exception(f"Failed to read FILE_INFO from {file_path}") return None @@ -40,8 +40,8 @@ def extract_run_number(file_path): def extract_threshold(file_path): """Extract threshold from the FILE_INFO table in the HDF5 file.""" try: - file_info = Table.read(file_path, path="FILE_INFO") - except OSError, ValueError, KeyError: + file_info = Table.read(file_path, path="FILE_INFO", format="hdf5") + except OSError, ValueError, KeyError, ImportError: _logger.exception(f"Failed to read FILE_INFO from {file_path}") return None diff --git a/src/simtools/visualization/plot_tables.py b/src/simtools/visualization/plot_tables.py index 01f4fd1175..9e5506f732 100644 --- a/src/simtools/visualization/plot_tables.py +++ b/src/simtools/visualization/plot_tables.py @@ -329,7 +329,7 @@ def resolve_plot_output_path(output, file_name="bias_curve.png"): return output_path / file_name -def plot_bias_curves(nsb_stats, proton_stats, config, output_path): +def plot_bias_curves(nsb_stats, proton_stats, config, output_path, trigger_threshold=None): """ Plot NSB and proton bias curves. @@ -343,12 +343,27 @@ def plot_bias_curves(nsb_stats, proton_stats, config, output_path): Plot configuration with title, ymin, and ymax. output_path : Path or str Output path for plot image. + trigger_threshold : float, optional + Trigger threshold value to annotate on the plot. """ fig, axis = plt.subplots(figsize=(10, 7)) _plot_nsb_curve(axis, nsb_stats) _plot_proton_curve(axis, proton_stats) - _configure_bias_curve_axis(axis, config) + _plot_scaled_proton_curve(axis, proton_stats) + + # Add vertical line at trigger threshold if provided (before legend creation) + if trigger_threshold is not None: + axis.axvline( + x=trigger_threshold, + color="grey", + linestyle="--", + linewidth=1.5, + alpha=0.5, + label=f"Trigger Threshold: {trigger_threshold:.1f}", + ) + + _configure_bias_curve_axis(axis, config, nsb_stats, proton_stats) fig.tight_layout() output_path = Path(output_path) @@ -401,6 +416,32 @@ def _plot_proton_curve(axis, proton_stats): _plot_log_linear_trend(axis, proton_thresholds, proton_rates, color="tab:orange") +def _plot_scaled_proton_curve(axis, proton_stats): + """Plot 1.35 * proton trigger rates for bias curve intersection.""" + if not proton_stats: + return + + proton_thresholds = sorted(proton_stats.keys()) + proton_rates = [proton_stats[t]["rate_hz"] for t in proton_thresholds] + proton_errors = [proton_stats[t]["error_hz"] for t in proton_thresholds] + + # Scale by 1.35 + scaled_proton_rates = [1.35 * r for r in proton_rates] + scaled_proton_errors = [1.35 * e for e in proton_errors] + + axis.errorbar( + proton_thresholds, + scaled_proton_rates, + yerr=scaled_proton_errors, + fmt="^", + label="1.35 x Proton", + color="tab:red", + capsize=3, + ) + + _plot_log_linear_trend(axis, proton_thresholds, scaled_proton_rates, color="tab:red") + + def _plot_log_linear_trend(axis, thresholds, rates, color): """Plot a log-linear trend line when at least two positive rates are available.""" if len(thresholds) < 2: @@ -420,13 +461,35 @@ def _plot_log_linear_trend(axis, thresholds, rates, color): axis.plot(x_fit, y_fit, "--", color=color, alpha=0.5, linewidth=1) -def _configure_bias_curve_axis(axis, config): +def _configure_bias_curve_axis(axis, config, nsb_stats, proton_stats): """Configure bias-curve axis labels, scaling, and legend.""" axis.set_title(config["title"], fontsize=14, fontweight="bold") axis.set_xlabel("Threshold", fontsize=12) axis.set_ylabel("Trigger Rate [Hz]", fontsize=12) axis.set_yscale("log") - axis.set_ylim(config["ymin"], config["ymax"]) + + # Dynamically set y-axis limits based on data + all_rates = [] + scaling_factor = config.get("scaling_factor", 1.35) + if nsb_stats: + all_rates.extend(stats["rate_hz"] for stats in nsb_stats.values() if stats["rate_hz"] > 0) + if proton_stats: + proton_rates = [stats["rate_hz"] for stats in proton_stats.values() if stats["rate_hz"] > 0] + all_rates.extend(proton_rates) + all_rates.extend(scaling_factor * rate for rate in proton_rates) + + if all_rates: + # Add 10% padding above and below + y_min = min(all_rates) * 0.9 + y_max = max(all_rates) * 1.1 + # Ensure reasonable minimum for log scale + if y_min <= 0: + y_min = 0.1 + axis.set_ylim(y_min, y_max) + else: + # Fallback to config values if no data + axis.set_ylim(config.get("ymin", 1), config.get("ymax", 1e6)) + axis.grid(which="both", alpha=0.3, linestyle=":") handles, _ = axis.get_legend_handles_labels() diff --git a/tests/unit_tests/simtel/test_bias_curve_generator.py b/tests/unit_tests/simtel/test_bias_curve_generator.py index 400ffbaf8f..866b5dcc00 100644 --- a/tests/unit_tests/simtel/test_bias_curve_generator.py +++ b/tests/unit_tests/simtel/test_bias_curve_generator.py @@ -17,9 +17,10 @@ def _base_args(tmp_path): "site": "North", "model_version": "7.0.0", "telescope": "LSTN-01", - "title": "Bias curve", - "ymin": 1, - "ymax": 1e6, + "scaling_factor": 1.35, + "parameter_version": "test", + "overwrite_model_parameters": False, + "title": "Trigger Rate Bias Curves", } @@ -271,7 +272,243 @@ def test_write_bias_curve_ecsv_writes_combined_table(tmp_path): assert "220 100.12 1.23 5.99 0.46" in output_text +def test_generate_bias_curves_raises_when_no_nsb_inputs_exist(tmp_path): + args = _base_args(tmp_path) + args["nsb_table_file"] = tmp_path / "nsb.ecsv" + + with ( + patch("simtools.simtel.bias_curve_generator._calculate_time_window", return_value=0.001), + patch("simtools.simtel.bias_curve_generator._extract_nsb_rates", return_value={}), + patch("simtools.simtel.bias_curve_generator._extract_proton_rates") as mock_extract_proton, + patch( + "simtools.simtel.bias_curve_generator.plot_tables.resolve_plot_output_path", + return_value=tmp_path / "bias.png", + ), + patch("simtools.simtel.bias_curve_generator.plot_tables.plot_bias_curves") as mock_plot, + patch("simtools.simtel.bias_curve_generator._write_bias_curve_ecsv") as mock_write_bias, + ): + with pytest.raises(FileNotFoundError, match="No NSB input files found"): + bias_curve_generator.generate_bias_curves(args) + + mock_extract_proton.assert_not_called() + mock_plot.assert_not_called() + mock_write_bias.assert_not_called() + + +def test_generate_bias_curves_raises_when_no_proton_inputs_exist(tmp_path): + args = _base_args(tmp_path) + + with ( + patch("simtools.simtel.bias_curve_generator._calculate_time_window", return_value=0.001), + patch( + "simtools.simtel.bias_curve_generator._extract_nsb_rates", + return_value={220: {"rate_hz": 100.0}}, + ), + patch("simtools.simtel.bias_curve_generator._extract_proton_rates", return_value={}), + patch("simtools.simtel.bias_curve_generator.plot_tables.plot_bias_curves") as mock_plot, + patch("simtools.simtel.bias_curve_generator._write_bias_curve_ecsv") as mock_write_bias, + ): + with pytest.raises(FileNotFoundError, match="No proton input files found"): + bias_curve_generator.generate_bias_curves(args) + + mock_plot.assert_not_called() + mock_write_bias.assert_not_called() + + +def test_calculate_trigger_threshold_success(): + """Test successful trigger threshold calculation with intersecting curves.""" + + # Create sample data where NSB and 1.35*proton curves intersect + # Need NSB > scaled_proton at low threshold and NSB < scaled_proton at high threshold + nsb_stats = { + 220: {"rate_hz": 1500.0, "error_hz": 50.0}, + 250: {"rate_hz": 1000.0, "error_hz": 30.0}, + 280: {"rate_hz": 400.0, "error_hz": 20.0}, # Lower NSB to ensure sign change + } + proton_stats = { + 220: {"rate_hz": 1100.0, "error_hz": 40.0}, + 250: {"rate_hz": 740.0, "error_hz": 25.0}, + 280: {"rate_hz": 370.0, "error_hz": 15.0}, + } + + args = {"scaling_factor": 1.35} + threshold = bias_curve_generator._calculate_trigger_threshold(args, nsb_stats, proton_stats) + # Should find intersection between 250 and 280 + assert threshold is not None + assert 250 <= threshold <= 280 + + +def test_calculate_trigger_threshold_no_valid_points(): + """Test that ValueError is raised when no valid threshold points exist.""" + args = {"scaling_factor": 1.35} + nsb_stats = {} + proton_stats = {} + + with pytest.raises(ValueError, match="No valid threshold points with both NSB and proton data"): + bias_curve_generator._calculate_trigger_threshold(args, nsb_stats, proton_stats) + + +def test_calculate_trigger_threshold_no_intersection(): + """Test that ValueError is raised when curves don't intersect.""" + args = {"scaling_factor": 1.35} + # NSB rates are always higher than 1.35*proton rates + nsb_stats = { + 220: {"rate_hz": 1000.0, "error_hz": 50.0}, + 250: {"rate_hz": 800.0, "error_hz": 30.0}, + 280: {"rate_hz": 600.0, "error_hz": 20.0}, + } + proton_stats = { + 220: {"rate_hz": 500.0, "error_hz": 40.0}, + 250: {"rate_hz": 400.0, "error_hz": 25.0}, + 280: {"rate_hz": 300.0, "error_hz": 15.0}, + } + + with pytest.raises(ValueError, match="Could not find intersection point"): + bias_curve_generator._calculate_trigger_threshold(args, nsb_stats, proton_stats) + + +def test_calculate_trigger_threshold_with_nan_values(): + """Test that NaN values are filtered out correctly.""" + + # Mix of valid and NaN values - need at least 2 valid points with sign change + nsb_stats = { + 220: {"rate_hz": np.nan, "error_hz": 50.0}, # Invalid + 250: {"rate_hz": 1000.0, "error_hz": 30.0}, # Valid + 280: {"rate_hz": 400.0, "error_hz": 20.0}, # Valid - adjusted to create sign change + } + proton_stats = { + 220: {"rate_hz": 700.0, "error_hz": 40.0}, + 250: {"rate_hz": 740.0, "error_hz": 25.0}, + 280: {"rate_hz": 370.0, "error_hz": 15.0}, + } + + # Should work with the valid points (250, 280) and find intersection + args = {"scaling_factor": 1.35} + threshold = bias_curve_generator._calculate_trigger_threshold(args, nsb_stats, proton_stats) + + assert threshold is not None + # Intersection should be between 250 and 280 + assert 250 <= threshold <= 280 + + +def test_export_trigger_threshold_as_model_parameter(tmp_path): + """Test successful export of trigger threshold as model parameter.""" + args = { + "telescope": "LSTN-01", + "site": "North", + "model_version": "7.0.0", + "parameter_version": "1.0.0", + } + + with patch("simtools.simtel.bias_curve_generator.io_handler.IOHandler") as mock_io: + mock_io.return_value.get_output_directory.return_value = tmp_path + + with patch( + "simtools.simtel.bias_curve_generator.model_data_writer.ModelDataWriter.write_model_parameter" + ) as mock_write: + with patch("simtools.simtel.bias_curve_generator.TelescopeModel") as mock_tel: + mock_model = MagicMock() + mock_model.get_parameter_value.return_value = "DigitalSum" + mock_tel.return_value = mock_model + + bias_curve_generator._export_trigger_threshold_as_model_parameter(args, 250.0) + + # Verify the model parameter was written with correct values + mock_write.assert_called_once() + call_kwargs = mock_write.call_args.kwargs + assert call_kwargs["parameter_name"] == "dsum_threshold" + assert call_kwargs["value"] == 250 # int for DigitalSum + assert call_kwargs["instrument"] == "LSTN-01" + assert call_kwargs["parameter_version"] == "1.0.0" + assert call_kwargs["unit"] == "count" + assert call_kwargs["check_db_for_existing_parameter"] is False + + +def test_export_trigger_threshold_uses_asum_for_analog_sum(tmp_path): + """Test that asum_threshold is used when default_trigger is AnalogSum.""" + args = { + "telescope": "LSTN-01", + "site": "North", + "model_version": "7.0.0", + "parameter_version": "1.0.0", + } + + with patch("simtools.simtel.bias_curve_generator.io_handler.IOHandler") as mock_io: + mock_io.return_value.get_output_directory.return_value = tmp_path + + with patch( + "simtools.simtel.bias_curve_generator.model_data_writer.ModelDataWriter.write_model_parameter" + ) as mock_write: + with patch("simtools.simtel.bias_curve_generator.TelescopeModel") as mock_tel: + mock_model = MagicMock() + mock_model.get_parameter_value.return_value = "AnalogSum" + mock_tel.return_value = mock_model + + bias_curve_generator._export_trigger_threshold_as_model_parameter(args, 250.5) + + # Verify the model parameter was written with correct values for AnalogSum + mock_write.assert_called_once() + call_kwargs = mock_write.call_args.kwargs + assert call_kwargs["parameter_name"] == "asum_threshold" + assert call_kwargs["value"] == pytest.approx(250.5, rel=0.01) # Float for AnalogSum + assert call_kwargs["instrument"] == "LSTN-01" + assert call_kwargs["parameter_version"] == "1.0.0" + assert call_kwargs["unit"] == "mV" + assert call_kwargs["check_db_for_existing_parameter"] is False + + +def test_export_trigger_threshold_handles_missing_telescope(tmp_path): + """Test that missing telescope name defaults to 'unknown' and function returns early.""" + args = {"parameter_version": "1.0.0"} # No telescope, no site, no model_version + + with patch("simtools.simtel.bias_curve_generator.io_handler.IOHandler") as mock_io: + mock_io.return_value.get_output_directory.return_value = tmp_path + + with patch( + "simtools.simtel.bias_curve_generator.model_data_writer.ModelDataWriter.write_model_parameter" + ) as mock_write: + with patch("simtools.simtel.bias_curve_generator._logger") as mock_logger: + bias_curve_generator._export_trigger_threshold_as_model_parameter(args, 250.0) + + # Should warn about missing telescope and then fail when trying to create TelescopeModel + assert mock_logger.warning.call_count >= 1 + # Write should not be called since we can't create TelescopeModel without site/model_version + mock_write.assert_not_called() + + +def test_export_trigger_threshold_handles_exception(tmp_path): + """Test that exceptions during export are caught and logged.""" + args = { + "telescope": "LSTN-01", + "parameter_version": "1.0.0", + "site": "North", + "model_version": "7.0.0", + } + + with patch("simtools.simtel.bias_curve_generator.io_handler.IOHandler") as mock_io: + mock_io.return_value.get_output_directory.return_value = tmp_path + + with patch( + "simtools.simtel.bias_curve_generator.model_data_writer.ModelDataWriter.write_model_parameter" + ) as mock_write: + mock_write.side_effect = OSError("Disk full") + + with patch("simtools.simtel.bias_curve_generator.TelescopeModel") as mock_tel: + mock_model = MagicMock() + mock_model.get_parameter_value.return_value = "DigitalSum" + mock_tel.return_value = mock_model + + with patch("simtools.simtel.bias_curve_generator._logger") as mock_logger: + bias_curve_generator._export_trigger_threshold_as_model_parameter(args, 250.0) + + # Should catch the exception and log warning + mock_logger.warning.assert_called_once_with( + "Failed to export trigger threshold as model parameter: Disk full" + ) + + def test_generate_bias_curves_runs_full_pipeline(tmp_path): + """Test that generate_bias_curves runs the full pipeline successfully.""" args = _base_args(tmp_path) args["proton_table_file"] = tmp_path / "proton.ecsv" args["nsb_table_file"] = tmp_path / "nsb.ecsv" @@ -295,52 +532,143 @@ def test_generate_bias_curves_runs_full_pipeline(tmp_path): ), patch("simtools.simtel.bias_curve_generator.plot_tables.plot_bias_curves") as mock_plot, patch("simtools.simtel.bias_curve_generator._write_bias_curve_ecsv") as mock_write_bias, + patch( + "simtools.simtel.bias_curve_generator._calculate_trigger_threshold", return_value=250.0 + ) as mock_calc_threshold, + patch( + "simtools.simtel.bias_curve_generator._export_trigger_threshold_as_model_parameter" + ) as mock_export, ): bias_curve_generator.generate_bias_curves(args) - mock_write_proton.assert_called_once() + mock_write_proton.assert_called_once_with( + {220: {"runs": {1: 10.0}, "rate_hz": 10.0, "error_hz": 0.0, "num_runs": 1}}, + tmp_path / "proton.ecsv", + ) mock_plot.assert_called_once() mock_write_bias.assert_called_once() + mock_calc_threshold.assert_called_once() + mock_export.assert_called_once() -def test_generate_bias_curves_raises_when_no_nsb_inputs_exist(tmp_path): +def test_run_nsb_trigger_derivation(tmp_path): + """Test _run_nsb_trigger_derivation function.""" args = _base_args(tmp_path) args["nsb_table_file"] = tmp_path / "nsb.ecsv" with ( - patch("simtools.simtel.bias_curve_generator._calculate_time_window", return_value=0.001), - patch("simtools.simtel.bias_curve_generator._extract_nsb_rates", return_value={}), - patch("simtools.simtel.bias_curve_generator._extract_proton_rates") as mock_extract_proton, patch( - "simtools.simtel.bias_curve_generator.plot_tables.resolve_plot_output_path", - return_value=tmp_path / "bias.png", - ), - patch("simtools.simtel.bias_curve_generator.plot_tables.plot_bias_curves") as mock_plot, - patch("simtools.simtel.bias_curve_generator._write_bias_curve_ecsv") as mock_write_bias, + "simtools.simtel.bias_curve_generator.derive_nsb_triggers", + return_value={220: {"rate_hz": 100.0}}, + ) as mock_derive, ): - with pytest.raises(FileNotFoundError, match="No NSB input files found"): - bias_curve_generator.generate_bias_curves(args) + result = bias_curve_generator._run_nsb_trigger_derivation(tmp_path, args, 0.001) + + assert result == {220: {"rate_hz": 100.0}} + mock_derive.assert_called_once_with( + { + "root_dir": tmp_path, + "pattern": "gamma*.reduced_event_data.hdf5", + "output": tmp_path / "nsb.ecsv", + "time_window": 0.001, + "verbose": False, + } + ) - mock_extract_proton.assert_not_called() - mock_plot.assert_not_called() - mock_write_bias.assert_not_called() +def test_calculate_proton_statistics_single_run(tmp_path): + """Test _calculate_proton_statistics_for_threshold with a single run.""" + files = {1: tmp_path / "run1.hdf5"} -def test_generate_bias_curves_raises_when_no_proton_inputs_exist(tmp_path): + with patch( + "simtools.simtel.bias_curve_generator._calculate_proton_rate_for_file", + return_value=10.0, + ): + stats = bias_curve_generator._calculate_proton_statistics_for_threshold(files, {}) + + assert stats["runs"] == {1: 10.0} + assert stats["rate_hz"] == pytest.approx(10.0) + assert stats["error_hz"] == pytest.approx(0.0) # Zero error for single run + assert stats["num_runs"] == 1 + + +def test_calculate_proton_rate_for_file_handles_oserror(tmp_path): + """Test _calculate_proton_rate_for_file handles OSError.""" args = _base_args(tmp_path) - with ( - patch("simtools.simtel.bias_curve_generator._calculate_time_window", return_value=0.001), - patch( - "simtools.simtel.bias_curve_generator._extract_nsb_rates", - return_value={220: {"rate_hz": 100.0}}, - ), - patch("simtools.simtel.bias_curve_generator._extract_proton_rates", return_value={}), - patch("simtools.simtel.bias_curve_generator.plot_tables.plot_bias_curves") as mock_plot, - patch("simtools.simtel.bias_curve_generator._write_bias_curve_ecsv") as mock_write_bias, + with patch( + "simtools.simtel.bias_curve_generator.telescope_trigger_rates", + side_effect=OSError("File not found"), ): - with pytest.raises(FileNotFoundError, match="No proton input files found"): - bias_curve_generator.generate_bias_curves(args) + result = bias_curve_generator._calculate_proton_rate_for_file( + tmp_path / "events.hdf5", args + ) + assert result is None - mock_plot.assert_not_called() - mock_write_bias.assert_not_called() + +def test_calculate_proton_rate_for_file_handles_keyerror(tmp_path): + """Test _calculate_proton_rate_for_file handles KeyError.""" + args = _base_args(tmp_path) + + with patch( + "simtools.simtel.bias_curve_generator.telescope_trigger_rates", + side_effect=KeyError("missing_key"), + ): + result = bias_curve_generator._calculate_proton_rate_for_file( + tmp_path / "events.hdf5", args + ) + assert result is None + + +def test_calculate_proton_rate_for_file_handles_valueerror(tmp_path): + """Test _calculate_proton_rate_for_file handles ValueError.""" + args = _base_args(tmp_path) + + with patch( + "simtools.simtel.bias_curve_generator.telescope_trigger_rates", + side_effect=ValueError("Invalid value"), + ): + result = bias_curve_generator._calculate_proton_rate_for_file( + tmp_path / "events.hdf5", args + ) + assert result is None + + +def test_calculate_proton_rate_for_file_handles_attributeerror(tmp_path): + """Test _calculate_proton_rate_for_file handles AttributeError.""" + args = _base_args(tmp_path) + + with patch( + "simtools.simtel.bias_curve_generator.telescope_trigger_rates", + side_effect=AttributeError("Missing attribute"), + ): + result = bias_curve_generator._calculate_proton_rate_for_file( + tmp_path / "events.hdf5", args + ) + assert result is None + + +def test_group_hdf5_files_skips_non_proton_files(tmp_test_directory): + """Test that _group_hdf5_files_by_threshold_and_run skips non-proton files.""" + from pathlib import Path + + tmp_path = Path(str(tmp_test_directory)) + _write_file_info_hdf5( + tmp_path / "gamma_run000001_asum220.reduced_event_data.hdf5", + "gamma_run000001_asum220.simtel.zst", + ) + + result = bias_curve_generator._group_hdf5_files_by_threshold_and_run(tmp_path) + assert result == {} # No proton files, should return empty dict + + +def test_group_hdf5_files_skips_missing_metadata(tmp_test_directory): + """Test that _group_hdf5_files_by_threshold_and_run skips files with missing metadata.""" + from pathlib import Path + + tmp_path = Path(str(tmp_test_directory)) + # Create a file without proper metadata + (tmp_path / "proton_run000001.reduced_event_data.hdf5").touch() + + result = bias_curve_generator._group_hdf5_files_by_threshold_and_run(tmp_path) + assert result == {} # File skipped due to missing threshold/run diff --git a/tests/unit_tests/visualization/test_plot_tables.py b/tests/unit_tests/visualization/test_plot_tables.py index 233ea7fb6c..cd5fa5f330 100644 --- a/tests/unit_tests/visualization/test_plot_tables.py +++ b/tests/unit_tests/visualization/test_plot_tables.py @@ -607,6 +607,7 @@ def test_resolve_plot_output_path_file_and_directory(tmp_path): @mock.patch("simtools.visualization.plot_tables._configure_bias_curve_axis") +@mock.patch("simtools.visualization.plot_tables._plot_scaled_proton_curve") @mock.patch("simtools.visualization.plot_tables._plot_proton_curve") @mock.patch("simtools.visualization.plot_tables._plot_nsb_curve") @mock.patch("simtools.visualization.plot_tables.plt") @@ -614,6 +615,7 @@ def test_plot_bias_curves_saves_and_closes_figure( mock_plt, mock_plot_nsb_curve, mock_plot_proton_curve, + mock_plot_scaled_proton_curve, mock_configure_axis, tmp_path, ): @@ -623,15 +625,16 @@ def test_plot_bias_curves_saves_and_closes_figure( output_path = tmp_path / "nested" / "bias.png" config = {"title": "Bias", "ymin": 1, "ymax": 1e6} + nsb_stats = {220: {"rate_hz": 10, "error_hz": 1}} + proton_stats = {220: {"rate_hz": 5}} - plot_tables.plot_bias_curves( - {220: {"rate_hz": 10, "error_hz": 1}}, {220: {"rate_hz": 5}}, config, output_path - ) + plot_tables.plot_bias_curves(nsb_stats, proton_stats, config, output_path) mock_plt.subplots.assert_called_once_with(figsize=(10, 7)) - mock_plot_nsb_curve.assert_called_once_with(mock_axis, {220: {"rate_hz": 10, "error_hz": 1}}) - mock_plot_proton_curve.assert_called_once_with(mock_axis, {220: {"rate_hz": 5}}) - mock_configure_axis.assert_called_once_with(mock_axis, config) + mock_plot_nsb_curve.assert_called_once_with(mock_axis, nsb_stats) + mock_plot_proton_curve.assert_called_once_with(mock_axis, proton_stats) + mock_plot_scaled_proton_curve.assert_called_once_with(mock_axis, proton_stats) + mock_configure_axis.assert_called_once_with(mock_axis, config, nsb_stats, proton_stats) mock_fig.tight_layout.assert_called_once_with() mock_fig.savefig.assert_called_once_with(output_path, dpi=200, bbox_inches="tight") mock_plt.close.assert_called_once_with(mock_fig) @@ -702,6 +705,39 @@ def test_plot_proton_curve_returns_early_when_empty(mock_plot_trend): mock_plot_trend.assert_not_called() +@mock.patch("simtools.visualization.plot_tables._plot_log_linear_trend") +def test_plot_scaled_proton_curve_draws_points_and_trend(mock_plot_trend): + axis = mock.MagicMock() + proton_stats = { + 240: {"rate_hz": 20.0, "error_hz": 2.0}, + 220: {"rate_hz": 10.0, "error_hz": 1.0}, + } + + plot_tables._plot_scaled_proton_curve(axis, proton_stats) + + # Rates should be scaled by 1.35 + axis.errorbar.assert_called_once_with( + [220, 240], + [13.5, 27.0], # 10.0 * 1.35, 20.0 * 1.35 + yerr=[1.35, 2.7], # 1.0 * 1.35, 2.0 * 1.35 + fmt="^", + label="1.35 x Proton", + color="tab:red", + capsize=3, + ) + mock_plot_trend.assert_called_once_with(axis, [220, 240], [13.5, 27.0], color="tab:red") + + +@mock.patch("simtools.visualization.plot_tables._plot_log_linear_trend") +def test_plot_scaled_proton_curve_returns_early_when_empty(mock_plot_trend): + axis = mock.MagicMock() + + plot_tables._plot_scaled_proton_curve(axis, {}) + + axis.errorbar.assert_not_called() + mock_plot_trend.assert_not_called() + + def test_plot_log_linear_trend_returns_without_plot_for_insufficient_data(): axis = mock.MagicMock() @@ -729,14 +765,15 @@ def test_configure_bias_curve_axis_with_legend(): axis = mock.MagicMock() axis.get_legend_handles_labels.return_value = ([object()], ["NSB"]) config = {"title": "Bias", "ymin": 1, "ymax": 1e6} + nsb_stats = {220: {"rate_hz": 10, "error_hz": 1}} + proton_stats = {220: {"rate_hz": 5}} - plot_tables._configure_bias_curve_axis(axis, config) + plot_tables._configure_bias_curve_axis(axis, config, nsb_stats, proton_stats) axis.set_title.assert_called_once_with("Bias", fontsize=14, fontweight="bold") axis.set_xlabel.assert_called_once_with("Threshold", fontsize=12) axis.set_ylabel.assert_called_once_with("Trigger Rate [Hz]", fontsize=12) axis.set_yscale.assert_called_once_with("log") - axis.set_ylim.assert_called_once_with(1, 1e6) axis.grid.assert_called_once_with(which="both", alpha=0.3, linestyle=":") axis.legend.assert_called_once_with(fontsize=11, loc="best") @@ -746,10 +783,88 @@ def test_configure_bias_curve_axis_without_legend_logs_warning(mock_logger): axis = mock.MagicMock() axis.get_legend_handles_labels.return_value = ([], []) config = {"title": "Bias", "ymin": 1, "ymax": 1e6} + nsb_stats = {} + proton_stats = {} - plot_tables._configure_bias_curve_axis(axis, config) + plot_tables._configure_bias_curve_axis(axis, config, nsb_stats, proton_stats) axis.legend.assert_not_called() mock_logger.warning.assert_called_once_with( "No NSB or proton rates found; writing empty bias-curve plot" ) + + +@mock.patch("simtools.visualization.plot_tables._configure_bias_curve_axis") +@mock.patch("simtools.visualization.plot_tables._plot_scaled_proton_curve") +@mock.patch("simtools.visualization.plot_tables._plot_proton_curve") +@mock.patch("simtools.visualization.plot_tables._plot_nsb_curve") +@mock.patch("simtools.visualization.plot_tables.plt") +def test_plot_bias_curves_with_trigger_threshold( + mock_plt, + mock_plot_nsb_curve, + mock_plot_proton_curve, + mock_plot_scaled_proton_curve, + mock_configure_axis, + tmp_path, +): + """Test that trigger threshold vertical line is added.""" + mock_fig = mock.MagicMock() + mock_axis = mock.MagicMock() + mock_plt.subplots.return_value = (mock_fig, mock_axis) + mock_axis.axvline.return_value = mock.MagicMock() + + output_path = tmp_path / "bias.png" + config = {"title": "Bias", "ymin": 1, "ymax": 1e6} + trigger_threshold = 250.0 + nsb_stats = {220: {"rate_hz": 10, "error_hz": 1}} + proton_stats = {220: {"rate_hz": 5}} + + plot_tables.plot_bias_curves( + nsb_stats, + proton_stats, + config, + output_path, + trigger_threshold, + ) + + mock_axis.axvline.assert_called_once() + call_kwargs = mock_axis.axvline.call_args[1] + assert call_kwargs["x"] == pytest.approx(trigger_threshold) + assert call_kwargs["color"] == "grey" + assert call_kwargs["linestyle"] == "--" + assert call_kwargs["linewidth"] == pytest.approx(1.5) + assert call_kwargs["alpha"] == pytest.approx(0.5) + + +@mock.patch("simtools.visualization.plot_tables._configure_bias_curve_axis") +@mock.patch("simtools.visualization.plot_tables._plot_scaled_proton_curve") +@mock.patch("simtools.visualization.plot_tables._plot_proton_curve") +@mock.patch("simtools.visualization.plot_tables._plot_nsb_curve") +@mock.patch("simtools.visualization.plot_tables.plt") +def test_plot_bias_curves_without_trigger_threshold( + mock_plt, + mock_plot_nsb_curve, + mock_plot_proton_curve, + mock_plot_scaled_proton_curve, + mock_configure_axis, + tmp_path, +): + """Test that vertical line is skipped when no trigger threshold is provided.""" + mock_fig = mock.MagicMock() + mock_axis = mock.MagicMock() + mock_plt.subplots.return_value = (mock_fig, mock_axis) + + output_path = tmp_path / "bias.png" + config = {"title": "Bias", "ymin": 1, "ymax": 1e6} + nsb_stats = {220: {"rate_hz": 10, "error_hz": 1}} + proton_stats = {220: {"rate_hz": 5}} + + plot_tables.plot_bias_curves( + nsb_stats, + proton_stats, + config, + output_path, + trigger_threshold=None, # No threshold provided + ) + + mock_axis.axvline.assert_not_called()