"""Light emission simulation (e.g. illuminators or flashers)."""
import logging
import shlex
from pathlib import Path
import astropy.units as u
import numpy as np
from simtools.atmosphere.transmission import AtmosphericTransmission
from simtools.context import RunContext
from simtools.job_execution import execution
from simtools.job_execution.process_accounting import build_accounting_command
from simtools.model.calibration_model import CalibrationModel
from simtools.model.model_utils import (
initialize_simulation_models,
read_overwrite_model_parameter_dict,
)
from simtools.model_repository.asset_names import get_simtel_table_file_name
from simtools.ray_tracing.backends.interface import LightSourceRequest, get_backend
from simtools.runners import runner_services
from simtools.runners.simtel_runner import SimtelRunner, sim_telarray_env_as_string
from simtools.sim_events import validate_run_outputs
from simtools.simtel.photon_bunch_file import (
bypass_optics_level,
camera_bunches,
write_photon_bunches,
)
from simtools.simtel.simtel_seeds import instance_seed_pairs
from simtools.simulation.configuration import get_light_source_writer, get_model_writer
from simtools.simulation.light_source_profiles import angular_profile, pulse_profile
from simtools.utils import value_conversion
from simtools.utils.geometry import (
fiducial_radius_from_shape,
ground_direction_to_telescope_frame,
ground_to_telescope_frame,
)
AXES_OFFSETS = "axes_offsets"
CORSIKA_OBSERVATION_LEVEL = "corsika_observation_level"
FLASHER_PHOTONS = "flasher_photons"
FLASHER_POSITION = "flasher_position"
FLASHER_WAVELENGTH = "flasher_wavelength"
FLAT_FIELDING = "flat_fielding"
LIGHT_SOURCE_TYPE = "light_source_type"
LIGHT_SOURCE_POSITION = "light_source_position"
LIGHT_SOURCE_POINTING = "light_source_pointing"
NUMBER_OF_EVENTS = "number_of_events"
RUN_NUMBER = "run_number"
TELESCOPE_SPHERE_RADIUS = "telescope_sphere_radius"
[docs]
class SimulatorLightEmission(SimtelRunner):
"""
Light emission simulation (e.g. illuminators or flashers).
Uses the sim_telarray LightEmission package to simulate the light emission.
Parameters
----------
light_emission_config : dict, optional
Configuration for the light emission (e.g. number of events, model names)
telescope : str, optional
Telescope name.
label : str, optional
Label for the simulation
"""
[docs]
@staticmethod
def get_available_wavelengths_from_config(config, model_reader=None, context=None):
"""
Get available wavelengths from model configuration without full initialization.
This is a lightweight method that only initializes the calibration model
to retrieve wavelengths, without telescope/site models or file I/O.
Parameters
----------
config : dict
Configuration dictionary with site, light_source, and model_version
Returns
-------
list of astropy.units.Quantity
List of available wavelengths with units
Raises
------
ValueError
If required configuration keys are missing
"""
site = config.get("site")
model_version = config.get("model_version")
light_source = config.get("light_source")
if site is None or model_version is None or light_source is None:
raise ValueError(
"Missing required configuration keys for wavelength query: "
f"site={site}, model_version={model_version}, light_source={light_source}"
)
calibration_model = CalibrationModel(
site=site,
calibration_device_model_name=light_source,
model_version=model_version,
label=f"temp_wavelength_query_{light_source}",
overwrite_model_parameter_dict=read_overwrite_model_parameter_dict(context=context),
model_reader=model_reader,
context=context,
)
wavelengths = calibration_model.get_parameter_value_with_unit(FLASHER_WAVELENGTH)
return value_conversion.ensure_list(wavelengths)
def __init__(
self, light_emission_config, telescope=None, label=None, model_reader=None, context=None
):
"""Initialize SimulatorLightEmission."""
light_emission_config = dict(light_emission_config)
self._logger = logging.getLogger(__name__)
self.configuration_writers = {}
self._seed_pair = None
self.context = context or RunContext.from_environment()
self.model_reader = model_reader or self.context.model_reader
self.io_handler = self.context.paths
telescope = telescope or light_emission_config.get("telescope")
label = f"{label}_{telescope}" if label else telescope
wavelength = light_emission_config.get("wavelength")
if wavelength is not None: # always include in the label for consistency
label = f"{label}_{wavelength.to(u.nm).value:.0f}nm"
super().__init__(label=label, config=light_emission_config, context=self.context)
self.submission_files = runner_services.RunnerServices(
light_emission_config, run_type="sub", label=label, context=self.context
)
model_kwargs = {
"label": label,
"site": light_emission_config.get("site"),
"telescope_name": telescope,
"calibration_device_name": light_emission_config.get("light_source"),
"calibration_device_type": light_emission_config.get(LIGHT_SOURCE_TYPE),
"model_version": light_emission_config.get("model_version"),
"model_reader": self.model_reader,
"context": self.context,
}
if light_emission_config.get("model_directory") is not None:
model_kwargs["model_directory"] = light_emission_config["model_directory"]
self.telescope_model, self.site_model, self.calibration_model = (
initialize_simulation_models(**model_kwargs)
)
get_model_writer(self.telescope_model).write_sim_telarray_config_file(
additional_models=self.site_model
)
self.light_emission_config = self._initialize_light_emission_configuration(
light_emission_config
)
def _initialize_light_emission_configuration(self, config):
"""Initialize light emission configuration from the calibration model."""
flasher_type = self.calibration_model.get_parameter_value("flasher_type")
if flasher_type:
config[LIGHT_SOURCE_TYPE] = flasher_type.lower()
if config.get("run_mode") is None:
config["run_mode"] = config[LIGHT_SOURCE_TYPE]
if config.get(FLASHER_PHOTONS) is None:
config[FLASHER_PHOTONS] = self.calibration_model.get_parameter_value(FLASHER_PHOTONS)
else:
photons = value_conversion.parse_typed_sequence(config[FLASHER_PHOTONS], int)
if len(photons) == 1:
self.calibration_model.overwrite_model_parameter(FLASHER_PHOTONS, photons[0])
photons = photons[0]
config[FLASHER_PHOTONS] = photons
if config.get("wavelength") is not None:
requested = config["wavelength"].to(u.nm)
allowed = [
wl.to(u.nm)
for wl in value_conversion.ensure_list(
self.calibration_model.get_parameter_value_with_unit(FLASHER_WAVELENGTH)
)
]
# Small tolerance for floating point comparison
if not any(abs(requested - wl) < 0.5 * u.nm for wl in allowed):
raise ValueError(
f"Wavelength {requested.value} nm is not supported by the "
f"illuminator. Allowed wavelengths are: {[wl.value for wl in allowed]} nm"
)
closest_wavelength = min(allowed, key=lambda wl: abs(wl - requested))
self.calibration_model.overwrite_model_parameter(FLASHER_WAVELENGTH, closest_wavelength)
config["wavelength"] = closest_wavelength
self._logger.info(f"Using wavelength: {closest_wavelength.value} nm")
if config.get(LIGHT_SOURCE_POSITION) is not None:
config[LIGHT_SOURCE_POSITION] = (
np.array(config[LIGHT_SOURCE_POSITION], dtype=float) * u.m
)
return config
@staticmethod
def _repeat_or_map_events(events, n_photon_levels):
"""Apply ff-1m event mapping rules to event counts."""
if n_photon_levels == 1:
return [events[0]]
if len(events) == n_photon_levels:
return events
if len(events) == 1:
return [events[0]] * n_photon_levels
raise ValueError(
"Invalid number_of_events list length. Use one value or one value per photon intensity."
)
[docs]
def build_flasher_event_and_photon_sequences(self):
"""Build ff-1m-compatible events/photons sequences."""
photons_int = value_conversion.parse_typed_sequence(
self.light_emission_config.get(FLASHER_PHOTONS), int
)
events = value_conversion.parse_typed_sequence(
self.light_emission_config.get(NUMBER_OF_EVENTS, 1), int
)
events_int = self._repeat_or_map_events(events, len(photons_int))
return events_int, photons_int
[docs]
def simulate(self):
"""Simulate light emission."""
run_script = self.prepare_run()
execution.run_command(
run_script,
stdout=self.submission_files.get_file_name("sub_out"),
stderr=self.submission_files.get_file_name("sub_err"),
)
[docs]
def prepare_run(self):
"""Write the bash run script with the light-emission command and return its path."""
script_file = self.submission_files.get_file_name(file_type="sub_script")
output_file = self.runner_service.get_file_name(file_type="sim_telarray_output")
if output_file.exists():
raise FileExistsError(
f"sim_telarray output file exists, cancelling simulation: {output_file}"
)
if self._backend() is not None:
self._write_camera_photons(self.runner_service.get_file_name(file_type="iact_output"))
lines = self.make_run_command()
script_file.write_text("".join(lines), encoding="utf-8")
return script_file
def _backend(self):
"""Return the ray-tracing backend class if it traces light, else None."""
backend = get_backend(context=self.context)
return backend if backend.traces_light_sources else None
def _seeds(self):
"""Return the (instrument seed, simulation seed) of this run."""
if self._seed_pair is None:
self._seed_pair = instance_seed_pairs(context=self.context)[0]
return self._seed_pair
def _dual_mirror(self):
"""Return True for dual-mirror telescopes."""
return self.telescope_model.get_parameter_value("mirror_class") == 2
def _bypass_optics(self):
"""Return sim_telarray's Bypass_Optics value: photons arrive at the camera entrance."""
if self._backend() is not None:
return bypass_optics_level(self._dual_mirror())
return 1 if self.light_emission_config[LIGHT_SOURCE_TYPE] == FLAT_FIELDING else 0
def _write_camera_photons(self, iact_output):
"""Write the photons the ray-tracing backend brings to the camera entrance."""
level_m = (
self.site_model.get_parameter_value_with_unit(CORSIKA_OBSERVATION_LEVEL).to(u.m).value
)
cal = self.calibration_model
shape, width_ns, decay_ns = cal.get_parameter_value("flasher_pulse_shape")
source = {
"angular_profile": angular_profile(
cal.get_parameter_value("flasher_angular_distribution"),
cal.get_parameter_value_with_unit("flasher_angular_distribution_width")
.to(u.deg)
.value,
),
"pulse": pulse_profile(shape, width_ns, decay_ns),
"wavelength_nm": cal.get_parameter_value_with_unit(FLASHER_WAVELENGTH).to(u.nm).value,
"seeds": self._seeds(),
}
if self.light_emission_config[LIGHT_SOURCE_TYPE] == FLAT_FIELDING:
request, radius, transmission, emission_altitude_m = self._flasher_request(
source, level_m
)
else:
request, radius, transmission, emission_altitude_m = self._illuminator_request(
source, level_m
)
backend = self._backend()(
self.telescope_model, self.site_model, self.label, request.telescope_zenith, False
)
dual_mirror = self._dual_mirror()
camera_rotation = (
self.telescope_model.get_parameter_value_with_unit("camera_rotate").to(u.deg).value
if dual_mirror
else 0.0
)
events = (
np.concatenate(
[
camera_bunches(
block | {"weight": block["weight"] * transmission},
dual_mirror=dual_mirror,
emission_altitude_m=emission_altitude_m,
camera_rotation_deg=camera_rotation,
)
for block in pulse
]
)
for pulse in backend.camera_photons(request)
)
write_photon_bunches(
iact_output,
events,
observation_level_m=level_m,
telescope_radius_m=radius,
atmospheric_profile=self.site_model.get_parameter_table("atmospheric_profile"),
run_number=self.light_emission_config.get(RUN_NUMBER) or 1,
)
def _photons_per_pulse(self):
"""Return the photons emitted in each light pulse (one entry per event)."""
events, photons = self.build_flasher_event_and_photon_sequences()
return tuple(n for count, n in zip(events, photons) for _ in range(count))
def _flasher_request(self, source, level_m):
"""Flasher light to the camera (as ff-1m); return (request, radius, 1, altitude)."""
x, y, z = (
c.to(u.m).value
for c in self.calibration_model.get_parameter_value_with_unit(FLASHER_POSITION)[:3]
)
# Primary-mirror frame to telescope frame.
z -= self.telescope_model.get_parameter_value_with_unit("mirror_offset").to(u.m).value
radius = fiducial_radius_from_shape(
self.telescope_model.get_parameter_value_with_unit("camera_body_diameter")
.to(u.m)
.value,
self.telescope_model.get_parameter_value("camera_body_shape"),
)
request = LightSourceRequest(
position_m=(x, y, z),
axis=(0.0, 0.0, 1.0),
photons_per_pulse=self._photons_per_pulse(),
telescope_zenith=0.0,
inside_telescope=True,
camera_radius_m=radius,
**source,
)
distance = self.calculate_distance_focal_plane_calibration_device().to(u.m).value
return request, radius, 1.0, level_m + distance
def _illuminator_request(self, source, level_m):
"""Illuminator light via the atmosphere; return (request, radius, transmission, level)."""
theta, phi = self._get_telescope_pointing()
position = self.get_illuminator_position()
telescope = self._get_telescope_position_ground_with_axis_offset(*position)
radius = self.telescope_model.get_parameter_value_with_unit(TELESCOPE_SPHERE_RADIUS)
telescope, radius = [c.to(u.m).value for c in telescope], radius.to(u.m).value
position = np.array([c.to(u.m).value for c in position])
relative = position - np.array(telescope)
if relative[2] == 0:
raise ValueError("Horizontal illuminator paths require a nonzero altitude difference")
pointing = np.asarray(self.get_illuminator_pointing_vector(position * u.m), dtype=float)
transmission = AtmosphericTransmission(
self.site_model.get_parameter_table("atmospheric_transmission")
).transmission(
source["wavelength_nm"],
level_m + position[2],
level_m + telescope[2],
np.linalg.norm(relative) / abs(relative[2]),
)
request = LightSourceRequest(
position_m=tuple(ground_to_telescope_frame(relative, theta, phi)),
axis=tuple(
ground_direction_to_telescope_frame(pointing / np.linalg.norm(pointing), theta, phi)
),
photons_per_pulse=self._photons_per_pulse(),
telescope_zenith=theta,
**source,
)
return request, radius, float(transmission), level_m
[docs]
def make_run_command(self, run_number=None, input_file=None): # pylint: disable=unused-argument
"""Light emission and sim_telarray run command (a backend writes the photons beforehand)."""
iact_output = self.runner_service.get_file_name(file_type="iact_output")
light_emission = (
""
if self._backend() is not None
else f"{get_light_source_writer(self).make_command(iact_output)}\n\n"
)
return [
"#!/usr/bin/env bash\n",
"set -e -o pipefail\n",
light_emission,
(
f"[ -s '{iact_output}' ] || "
f"{{ echo 'LightEmission did not produce IACT file' >&2; exit 1; }}\n\n"
),
f"{self._make_simtel_script()}\n\n",
f"rm -f '{iact_output}'\n\n",
]
def _get_telescope_pointing(self):
"""Return (theta, phi) telescope pointing for the light source type and position."""
if self.light_emission_config[LIGHT_SOURCE_TYPE] == FLAT_FIELDING:
return 0.0, 0.0
if self.light_emission_config.get(LIGHT_SOURCE_POSITION) is not None:
self._logger.info("Using fixed (vertical up) telescope pointing.")
return 0.0, 0.0
illuminator_position = self.get_illuminator_position()
_, angles = self._calibration_pointing_direction(*illuminator_position)
return angles[0], angles[1]
def _calibration_pointing_direction(self, x_cal=None, y_cal=None, z_cal=None):
"""Return [pointing vector, telescope and beam angles] for an off-telescope device."""
if x_cal is None or y_cal is None or z_cal is None:
x_cal, y_cal, z_cal = self.calibration_model.get_parameter_value_with_unit(
"array_element_position_ground"
)
x_cal, y_cal, z_cal = [coord.to(u.m).value for coord in (x_cal, y_cal, z_cal)]
cal_vect = np.array([x_cal, y_cal, z_cal])
x_tel, y_tel, z_tel = self._get_telescope_position_ground_with_axis_offset(
x_cal=x_cal * u.m,
y_cal=y_cal * u.m,
z_cal=z_cal * u.m,
)
x_tel, y_tel, z_tel = [coord.to(u.m).value for coord in (x_tel, y_tel, z_tel)]
tel_vect = np.array([x_tel, y_tel, z_tel])
direction_vector = tel_vect - cal_vect
# pointing vector from calibration device to telescope
pointing_vector = np.round(direction_vector / np.linalg.norm(direction_vector), 6)
# Calculate telescope theta and phi angles
tel_theta = 180 - np.round(
np.rad2deg(np.arccos(direction_vector[2] / np.linalg.norm(direction_vector))), 6
)
tel_phi = 180 - np.round(
np.rad2deg(np.arctan2(direction_vector[1], direction_vector[0])), 6
)
# Calculate source beam theta and phi angles
direction_vector_inv = direction_vector * -1
source_theta = np.round(
np.rad2deg(np.arccos(direction_vector_inv[2] / np.linalg.norm(direction_vector_inv))),
6,
)
source_phi = np.round(
np.rad2deg(np.arctan2(direction_vector_inv[1], direction_vector_inv[0])), 6
)
return pointing_vector.tolist(), [tel_theta, tel_phi, source_theta, source_phi]
def _get_telescope_position_ground_with_axis_offset(self, x_cal=None, y_cal=None, z_cal=None):
"""Return the telescope ground position corrected for the illuminator-axis offsets."""
if self.light_emission_config.get(LIGHT_SOURCE_POSITION) is not None and not (
self.uses_telescope_position_file(self.light_emission_config.get(LIGHT_SOURCE_POINTING))
):
return (0.0 * u.m,) * 3
x_tel, y_tel, z_tel = self.telescope_model.get_parameter_value_with_unit(
"array_element_position_ground"
)
axes_offsets = self.telescope_model.get_parameter_value_with_unit(AXES_OFFSETS)
offset_1 = axes_offsets[0].to(u.m).value
offset_2 = axes_offsets[1].to(u.m).value if len(axes_offsets) > 1 else 0.0
if np.isclose(offset_1, 0.0) and np.isclose(offset_2, 0.0):
return x_tel, y_tel, z_tel
if x_cal is None or y_cal is None or z_cal is None:
x_cal, y_cal, z_cal = self.calibration_model.get_parameter_value_with_unit(
"array_element_position_ground"
)
x_tel_m = x_tel.to(u.m).value
y_tel_m = y_tel.to(u.m).value
z_tel_m = z_tel.to(u.m).value
delta_x = x_cal.to(u.m).value - x_tel_m
delta_y = y_cal.to(u.m).value - y_tel_m
delta_z = z_cal.to(u.m).value - z_tel_m
norm_horizontal = np.hypot(delta_x, delta_y)
# Angle from North (y-axis) towards East (x-axis)
phi = np.arctan2(delta_x, delta_y)
el = np.arctan2(delta_z, norm_horizontal)
# Unit vectors in telescope-pointing convention (phi, el).
u_pointing = np.array(
[
np.cos(el) * np.sin(phi),
np.cos(el) * np.cos(phi),
np.sin(el),
],
dtype=float,
)
u_xy = np.array([np.sin(phi), np.cos(phi), 0.0], dtype=float)
altitude_axis = np.array([np.cos(phi), -np.sin(phi), 0.0], dtype=float)
v_perpendicular = np.cross(altitude_axis, u_pointing)
norm_v = np.linalg.norm(v_perpendicular)
v_perpendicular /= norm_v
corrected_position = np.array([x_tel_m, y_tel_m, z_tel_m], dtype=float)
corrected_position += offset_1 * u_xy + offset_2 * v_perpendicular
return (
corrected_position[0] * u.m,
corrected_position[1] * u.m,
corrected_position[2] * u.m,
)
[docs]
def write_telescope_position_file(self, illuminator_position=None):
"""
Write the telescope positions to a telescope_position file.
The file will contain lines in the format: x y z r in cm
Returns
-------
Path
The path to the generated telescope_position file.
"""
if illuminator_position is None:
illuminator_position = self.get_illuminator_position()
x_cal, y_cal, z_cal = illuminator_position
x_tel, y_tel, z_tel = self._get_telescope_position_ground_with_axis_offset(
x_cal=x_cal,
y_cal=y_cal,
z_cal=z_cal,
)
x_tel, y_tel, z_tel = [coord.to(u.cm).value for coord in (x_tel, y_tel, z_tel)]
radius = self.telescope_model.get_parameter_value_with_unit(TELESCOPE_SPHERE_RADIUS)
radius = radius.to(u.cm).value # Convert radius to cm
tel = self.get_file_name_token(self.light_emission_config.get("telescope") or "telescope")
cal = self.get_file_name_token(
self.light_emission_config.get("light_source") or "calibration"
)
telescope_position_file = (
self.io_handler.get_output_directory("light_emission")
/ f"telescope_position_{tel}_{cal}.dat"
)
telescope_position_file.write_text(f"{x_tel} {y_tel} {z_tel} {radius}\n", encoding="utf-8")
return telescope_position_file
[docs]
def get_illuminator_position(self):
"""Return the illuminator emission position in ground coordinates."""
pos = self.light_emission_config.get(LIGHT_SOURCE_POSITION)
if pos is not None:
return pos
x_pos, y_pos, z_pos = self.calibration_model.get_parameter_value_with_unit(
"array_element_position_ground"
)
tower_height = self.calibration_model.get_parameter_value_with_unit(
"illuminator_tower_height"
)
return x_pos, y_pos, z_pos + tower_height
[docs]
def get_illuminator_pointing_vector(self, pos=None):
"""Return illuminator pointing vector; prefer explicit config if available."""
pointing_vector = self.light_emission_config.get(LIGHT_SOURCE_POINTING)
if pointing_vector is not None:
return pointing_vector
if pos is None:
pos = self.get_illuminator_position()
x_cal, y_cal, z_cal = pos
return self._calibration_pointing_direction(x_cal, y_cal, z_cal)[0]
[docs]
@staticmethod
def uses_telescope_position_file(pointing_vector):
"""Return False for the downward beam used with an artificial telescope at the origin."""
try:
vec = np.asarray(pointing_vector, dtype=float)
except TypeError, ValueError:
return True
if vec.size < 3 or not np.all(np.isfinite(vec[:3])):
return True
is_default_down = np.allclose(vec[:3], [0.0, 0.0, -1.0], atol=1e-6)
return not is_default_down
[docs]
@staticmethod
def get_file_name_token(value):
"""Return a filename token for a light source or telescope name."""
return "".join(ch if (ch.isalnum() or ch in ("-", "_")) else "_" for ch in str(value))
def _make_simtel_script(self):
"""Return the command to run sim_telarray on the output of the previous step."""
bypass = self._bypass_optics()
theta, phi = (0.0, 0.0) if bypass else self._get_telescope_pointing()
simtel_bin = str(self.context.software.sim_telarray_exe)
parts = [
simtel_bin,
f"-I{self.telescope_model.config_file_directory}",
"-c",
str(self.telescope_model.config_file_path),
"-DNUM_TELESCOPES=1",
]
altitude = (
self.site_model.get_parameter_value_with_unit(CORSIKA_OBSERVATION_LEVEL).to(u.m).value
)
atmosphere = self.site_model.get_parameter_value("atmospheric_transmission")
if str(atmosphere).lower().endswith(".ecsv"): # generated sim_telarray table
parameter_data = self.site_model.parameters.get("atmospheric_transmission", {})
atmosphere = get_simtel_table_file_name(parameter_data) or (
f"atmospheric_transmission-{Path(self.telescope_model.config_file_path).stem}.dat"
)
options = [
("altitude", altitude),
("atmospheric_transmission", atmosphere),
("TRIGGER_TELESCOPES", "1"),
("TELTRIG_MIN_SIGSUM", "2"),
("PULSE_ANALYSIS", "-30"),
("MAXIMUM_TELESCOPES", 1),
("telescope_theta", f"{theta}"),
("telescope_phi", f"{phi}"),
]
if bypass:
options.append(("Bypass_Optics", f"{bypass}"))
options.append(("random_seed", ",".join(str(seed) for seed in self._seeds())))
input_file = self.runner_service.get_file_name(file_type="iact_output")
output_file = self.runner_service.get_file_name(file_type="sim_telarray_output")
histo_file = self.runner_service.get_file_name(file_type="sim_telarray_histogram")
options += [
("power_law", "2.68"),
("input_file", f"{input_file}"),
("output_file", f"{output_file}"),
("histogram_file", f"{histo_file}"),
]
for key, value in options:
parts.extend(["-C", f"{key}={value}"])
log_file = self.runner_service.get_file_name(file_type="sim_telarray_log")
command = build_accounting_command(
parts,
self.runner_service.get_file_name(file_type="sim_telarray_resources"),
"sim_telarray",
self.light_emission_config.get(RUN_NUMBER, 1),
log_file=log_file,
)
return sim_telarray_env_as_string() + shlex.join(command) + "\n"
[docs]
def calculate_distance_focal_plane_calibration_device(self):
"""
Return the distance between the focal plane and a flasher (focal length + focus offset - z).
Returns
-------
astropy.units.Quantity
Distance between calibration device and focal plane.
"""
focal_length = self.telescope_model.get_parameter_value_with_unit("focal_length").to(u.m)
flasher_z = self.calibration_model.get_parameter_value_with_unit(FLASHER_POSITION)[2].to(
u.m
)
return focal_length + self._vertical_focus_offset() - flasher_z
def _vertical_focus_offset(self):
"""Return the vertical-telescope focus offset (zenith terms only if both are non-zero)."""
values = self.telescope_model.get_parameter_value_with_unit("focus_offset")
if values is None:
return 0.0 * u.m
values = value_conversion.ensure_list(values)
offset = values[0].to(u.m)
if len(values) < 4 or not (values[2].value and values[3].value):
return offset
reference = values[1].to(u.rad).value
return (
offset
+ (values[2].value * (1.0 - np.cos(reference)) - values[3].value * np.sin(reference))
* values[0].unit.to(u.m)
* u.m
)
[docs]
def validate_simulations(self):
"""Validate that the simulations were successful."""
validate_run_outputs.validate_sim_telarray(
data_files=Path(self.runner_service.get_file_name(file_type="sim_telarray_output")),
log_files=Path(self.runner_service.get_file_name(file_type="sim_telarray_log")),
array_models=None,
)
