diff --git a/ctapipe/atmosphere.py b/ctapipe/atmosphere.py new file mode 100644 index 00000000000..f544c634dd7 --- /dev/null +++ b/ctapipe/atmosphere.py @@ -0,0 +1,396 @@ +#!/usr/bin/env python3 + +"""Atmosphere density models and functions to transform between column density +(X in grammage units) and height (meters) units. + +Zenith angle is taken into account in the line-of-sight integral to compute the +column density X assuming Earth as a flat plane (the curvature is not taken into +account) + +""" + +import abc +from dataclasses import dataclass +from functools import partial +from typing import Dict + +import numpy as np +from astropy import units as u +from astropy.table import Table +from scipy.interpolate import interp1d + +__all__ = [ + "AtmosphereDensityProfile", + "ExponentialAtmosphereDensityProfile", + "TableAtmosphereDensityProfile", + "FiveLayerAtmosphereDensityProfile", +] + +SUPPORTED_TABLE_VERSIONS = { + 1, +} + + +class AtmosphereDensityProfile(abc.ABC): + """ + Base class for models of atmosphere density. + """ + + @abc.abstractmethod + def __call__(self, height: u.Quantity) -> u.Quantity: + """ + Returns + ------- + u.Quantity["g cm-3"] + the density at height h + """ + + @abc.abstractmethod + def integral(self, height: u.Quantity) -> u.Quantity: + r"""Integral of the profile along the height axis, i.e. the *atmospheric + depth* :math:`X`. + + .. math:: X(h) = \int_{h}^{\infty} \rho(h') dh' + + Returns + ------- + u.Quantity["g/cm2"]: + Integral of the density from height h to infinity + + """ + + def line_of_sight_integral( + self, distance: u.Quantity, zenith_angle=0 * u.deg, output_units=u.g / u.cm**2 + ): + r"""Line-of-sight integral from the shower distance to infinity, along + the direction specified by the zenith angle. This is sometimes called + the *slant depth*. The atmosphere here is assumed to be Cartesian, the + curvature of the Earth is not taken into account. + + .. math:: X(h, \Psi) = \int_{h}^{\infty} \rho(h' \cos{\Psi}) dh' + + Parameters + ---------- + distance: u.Quantity["length"] + line-of-site distance from observer to point + zenith_angle: u.Quantity["angle"] + zenith angle of observation + output_units: u.Unit + unit to output (must be convertible to g/cm2) + + """ + + return ( + self.integral(distance * np.cos(zenith_angle)) / np.cos(zenith_angle) + ).to(output_units) + + def peek(self): + """ + Draw quick plot of profile + """ + # pylint: disable=import-outside-toplevel + import matplotlib.pyplot as plt + + fig, axis = plt.subplots(1, 3, constrained_layout=True, figsize=(10, 3)) + + fig.suptitle(self.__class__.__name__) + height = np.linspace(1, 100, 500) * u.km + density = self(height) + axis[0].set_xscale("linear") + axis[0].set_yscale("log") + axis[0].plot(height, density) + axis[0].set_xlabel(f"Height / {height.unit.to_string('latex')}") + axis[0].set_ylabel(f"Density / {density.unit.to_string('latex')}") + axis[0].grid(True) + + distance = np.linspace(1, 100, 500) * u.km + for zenith_angle in [0, 40, 50, 70] * u.deg: + column_density = self.line_of_sight_integral(distance, zenith_angle) + axis[1].plot(distance, column_density, label=f"$\\Psi$={zenith_angle}") + + axis[1].legend(loc="best") + axis[1].set_xlabel(f"Distance / {distance.unit.to_string('latex')}") + axis[1].set_ylabel(f"Column Density / {column_density.unit.to_string('latex')}") + axis[1].set_yscale("log") + axis[1].grid(True) + + zenith_angle = np.linspace(0, 80, 20) * u.deg + for distance in [0, 5, 10, 20] * u.km: + column_density = self.line_of_sight_integral(distance, zenith_angle) + axis[2].plot(zenith_angle, column_density, label=f"Height={distance}") + + axis[2].legend(loc="best") + axis[2].set_xlabel( + f"Zenith Angle $\\Psi$ / {zenith_angle.unit.to_string('latex')}" + ) + axis[2].set_ylabel(f"Column Density / {column_density.unit.to_string('latex')}") + axis[2].set_yscale("log") + axis[2].grid(True) + + plt.show() + return fig, axis + + @classmethod + def from_table(cls, table: Table): + """return a subclass of AtmosphereDensityProfile from a serialized + table""" + + if "TAB_TYPE" not in table.meta: + raise ValueError("expected a TAB_TYPE metadata field") + + version = table.meta.get("TAB_VER", "") + if version not in SUPPORTED_TABLE_VERSIONS: + raise ValueError(f"Table version not supported: '{version}'") + + tabtype = table.meta.get("TAB_TYPE") + + if tabtype == "ctapipe.atmosphere.TableAtmosphereDensityProfile": + return TableAtmosphereDensityProfile(table) + if tabtype == "ctapipe.atmosphere.FiveLayerAtmosphereDensityProfile": + return FiveLayerAtmosphereDensityProfile(table) + + raise TypeError(f"Unknown AtmosphereDensityProfile type: '{tabtype}'") + + +@dataclass +class ExponentialAtmosphereDensityProfile(AtmosphereDensityProfile): + """ + A simple functional density profile modeled as an exponential. + + The is defined following the form: + + .. math:: \\rho(h) = \\rho_0 e^{-h/h_0} + + + .. code-block:: python + + from ctapipe.atmosphere import ExponentialAtmosphereDensityProfile + density_profile = ExponentialAtmosphereDensityProfile() + density_profile.peek() + + + Attributes + ---------- + scale_height: u.Quantity["m"] + scale height (h0) + scale_density: u.Quantity["g cm-3"] + scale density (rho0) + """ + + scale_height: u.Quantity = 8 * u.km + scale_density: u.Quantity = 0.00125 * u.g / (u.cm**3) + + @u.quantity_input(height=u.m) + def __call__(self, height) -> u.Quantity: + return self.scale_density * np.exp(-height / self.scale_height) + + @u.quantity_input(height=u.m) + def integral( + self, + height, + ) -> u.Quantity: + return ( + self.scale_density * self.scale_height * np.exp(-height / self.scale_height) + ) + + +class TableAtmosphereDensityProfile(AtmosphereDensityProfile): + """Tabular profile from a table that has both the density and it's integral + pre-computed. The table is interpolated to return the density and its integral. + + .. code-block:: python + + from astropy.table import Table + from astropy import units as u + + from ctapipe.atmosphere import TableAtmosphereDensityProfile + + table = Table( + dict( + height=[1,10,20] * u.km, + density=[0.00099,0.00042, 0.00009] * u.g / u.cm**3 + column_density=[1044.0, 284.0, 57.0] * u.g / u.cm**2 + ) + ) + + profile = TableAtmosphereDensityProfile(table=table) + print(profile(10 * u.km)) + + + Attributes + ---------- + table: Table + Points that define the model + + See Also + -------- + ctapipe.io.eventsource.EventSource.atmosphere_density_profile: + load a TableAtmosphereDensityProfile from a supported EventSource + """ + + def __init__(self, table: Table): + """ + Parameters + ---------- + table: Table + Table with columns `height`, `density`, and `column_density` + """ + + for col in ["height", "density", "column_density"]: + if col not in table.colnames: + raise ValueError(f"Missing expected column: {col} in table") + + self.table = table[ + (table["height"] >= 0) + & (table["density"] > 0) + & (table["column_density"] > 0) + ] + + # interpolation is done in log-y to minimize spline wobble + + self._density_interp = interp1d( + self.table["height"].to("km").value, + np.log10(self.table["density"].to("g cm-3").value), + kind="cubic", + ) + self._col_density_interp = interp1d( + self.table["height"].to("km").value, + np.log10(self.table["column_density"].to("g cm-2").value), + kind="cubic", + ) + + # ensure it can be read back + self.table.meta["TAB_TYPE"] = "ctapipe.atmosphere.TableAtmosphereDensityProfile" + self.table.meta["TAB_VER"] = 1 + + @u.quantity_input(height=u.m) + def __call__(self, height) -> u.Quantity: + return u.Quantity( + 10 ** self._density_interp(height.to_value(u.km)), u.g / u.cm**3 + ) + + @u.quantity_input(height=u.m) + def integral(self, height) -> u.Quantity: + return u.Quantity( + 10 ** self._col_density_interp(height.to_value(u.km)), u.g / u.cm**2 + ) + + def __repr__(self): + return ( + f"{self.__class__.__name__}(meta={self.table.meta}, rows={len(self.table)})" + ) + + +# Here we define some utility functions needed to build the piece-wise 5-layer +# model. + +# pylint: disable=invalid-name,unused-argument +def _exponential(h, a, b, c): + """exponential atmosphere""" + return a + b * np.exp(-h / c) + + +def _d_exponential(h, a, b, c): + """derivative of exponential atmosphere""" + return -b / c * np.exp(-h / c) + + +def _linear(h, a, b, c): + """linear atmosphere""" + return a - b * h / c + + +def _d_linear(h, a, b, c): + """derivative of linear atmosphere""" + return -b / c + + +class FiveLayerAtmosphereDensityProfile(AtmosphereDensityProfile): + r""" + CORSIKA 5-layer atmosphere model + + Layers 1-4 are modeled with: + + .. math:: T(h) = a_i + b_i \exp{-h/c_i} + + Layer 5 is modeled with: + + ..math:: T(h) = a_5 - b_5 \frac{h}{c_5} + + References + ---------- + [corsika-user] D. Heck and T. Pierog, "Extensive Air Shower Simulation with CORSIKA: + A User’s Guide", 2021, Appendix F + """ + + def __init__(self, table: Table): + self.table = table + + param_a = self.table["a"].to("g/cm2") + param_b = self.table["b"].to("g/cm2") + param_c = self.table["c"].to("km") + + # build list of column density functions and their derivatives: + self._funcs = [ + partial(f, a=param_a[i], b=param_b[i], c=param_c[i]) + for i, f in enumerate([_exponential] * 4 + [_linear]) + ] + self._d_funcs = [ + partial(f, a=param_a[i], b=param_b[i], c=param_c[i]) + for i, f in enumerate([_d_exponential] * 4 + [_d_linear]) + ] + + @classmethod + def from_array(cls, array: np.ndarray, metadata: Dict = None): + """construct from a 5x5 array as provided by eventio""" + + if metadata is None: + metadata = {} + + if array.shape != (5, 5): + raise ValueError("expected ndarray with shape (5,5)") + + table = Table( + array, + names=["height", "a", "b", "c", "1/c"], + units=["cm", "g/cm2", "g/cm2", "cm", "cm-1"], + meta=metadata, + ) + + table.meta.update( + dict( + TAB_VER=1, + TAB_TYPE="ctapipe.atmosphere.FiveLayerAtmosphereDensityProfile", + ) + ) + return cls(table) + + @u.quantity_input(height=u.m) + def __call__(self, height) -> u.Quantity: + which_func = np.digitize(height, self.table["height"]) - 1 + condlist = [which_func == i for i in range(5)] + return u.Quantity( + -1 + * np.piecewise( + height, + condlist=condlist, + funclist=self._d_funcs, + ) + ).to(u.g / u.cm**3) + + @u.quantity_input(height=u.m) + def integral(self, height) -> u.Quantity: + which_func = np.digitize(height, self.table["height"]) - 1 + condlist = [which_func == i for i in range(5)] + return u.Quantity( + np.piecewise( + x=height, + condlist=condlist, + funclist=self._funcs, + ) + ).to(u.g / u.cm**2) + + def __repr__(self): + return ( + f"{self.__class__.__name__}(meta={self.table.meta}, rows={len(self.table)})" + ) diff --git a/ctapipe/instrument/atmosphere.py b/ctapipe/instrument/atmosphere.py index e0f85c6002e..76a9e467554 100644 --- a/ctapipe/instrument/atmosphere.py +++ b/ctapipe/instrument/atmosphere.py @@ -3,6 +3,7 @@ """ import numpy as np from astropy.units import Quantity +from astropy.utils.decorators import deprecated from scipy.interpolate import interp1d from ctapipe.utils import get_table_dataset @@ -10,6 +11,9 @@ __all__ = ["get_atmosphere_profile_table", "get_atmosphere_profile_functions"] +@deprecated( + since="v0.16.0", message="use ctapipe.atmosphere profiles from EventSource instead" +) def get_atmosphere_profile_table(atmosphere_name="paranal"): """ Get an atmosphere profile table @@ -30,6 +34,9 @@ def get_atmosphere_profile_table(atmosphere_name="paranal"): return table +@deprecated( + since="v0.16.0", message="use ctapipe.atmosphere profiles from EventSource instead" +) def get_atmosphere_profile_functions(atmosphere_name="paranal", with_units=True): """ Gives atmospheric profile as a continuous function thickness( diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py index ecc295b8912..14f33a22537 100644 --- a/ctapipe/io/datawriter.py +++ b/ctapipe/io/datawriter.py @@ -23,6 +23,7 @@ from ..instrument import SubarrayDescription from . import EventSource, HDF5TableWriter, TableWriter from . import metadata as meta +from .astropy_helpers import write_table from .datalevels import DataLevel from .simteleventsource import SimTelEventSource from .tableio import FixedPointColumnTransform, TelListToMaskTransform @@ -273,6 +274,9 @@ def _setup_outputfile(self): self._write_scheduling_and_observation_blocks() if self._is_simulation: self._write_simulation_configuration() + self._write_atmosphere_profile( + "/simulation/service/atmosphere_density_profile" + ) def __enter__(self): return self @@ -797,3 +801,31 @@ def _write_context_metadata_headers(self): context_dict[key] = value meta.write_to_hdf5(context_dict, self._writer.h5file) + + def _write_atmosphere_profile(self, path): + """ + write atmosphere profiles if they are in a tabular format + + Parameters + ---------- + path: str + path in the HDF5 file where to place the profile + + """ + + profile = self.event_source.atmosphere_density_profile + + if profile: + if hasattr(profile, "table"): + write_table( + table=profile.table, + h5file=self._writer.h5file, + path=path, + append=False, + ) + else: + self.logger.warning( + f"The AtmosphereDensityProfile type '{profile.__class__.__name__}' " + "is not serializable. No atmosphere profile will be stored in the " + "output file" + ) diff --git a/ctapipe/io/eventsource.py b/ctapipe/io/eventsource.py index 2ce3bf9f25a..1debd1ca144 100644 --- a/ctapipe/io/eventsource.py +++ b/ctapipe/io/eventsource.py @@ -7,6 +7,8 @@ from traitlets.config.loader import LazyConfigValue +from ctapipe.atmosphere import AtmosphereDensityProfile + from ..containers import ( ArrayEventContainer, ObservationBlockContainer, @@ -280,6 +282,20 @@ def obs_ids(self) -> List[int]: """ return list(self.observation_blocks.keys()) + @property + def atmosphere_density_profile(self) -> AtmosphereDensityProfile: + """atmosphere density profile that can be integrated to + convert between h_max and X_max. This should correspond + either to what was used in a simualtion, or a measurment + for use with observed data. + + Returns + ------- + AtmosphereDensityProfile: + profile to be used + """ + return None + @abstractmethod def _generator(self) -> Generator[ArrayEventContainer, None, None]: """ diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py index d92f7c9c67a..ccae1d762e7 100644 --- a/ctapipe/io/hdf5eventsource.py +++ b/ctapipe/io/hdf5eventsource.py @@ -8,6 +8,7 @@ import tables from astropy.utils.decorators import lazyproperty +from ctapipe.atmosphere import AtmosphereDensityProfile from ctapipe.instrument.optics import FocalLengthKind from ..containers import ( @@ -42,6 +43,7 @@ from ..core.traits import UseEnum from ..instrument import SubarrayDescription from ..utils import IndexFinder +from .astropy_helpers import read_table from .datalevels import DataLevel from .eventsource import EventSource from .hdf5tableio import HDF5TableReader, get_column_attrs @@ -86,6 +88,33 @@ def get_hdf5_datalevels(h5file): return tuple(datalevels) +def read_atmosphere_density_profile( + h5file: tables.File, path="/simulation/service/atmosphere_density_profile" +): + """return a subclass of AtmosphereDensityProfile by + reading a table in a h5 file + + Parameters + ---------- + h5file: tables.File + file handle of HDF5 file + path: str + path in the file where the serialized model is stored as an + astropy.table.Table + + Returns + ------- + AtmosphereDensityProfile: + subclass depending on type of table + """ + + if path not in h5file: + return None + + table = read_table(h5file=h5file, path=path) + return AtmosphereDensityProfile.from_table(table) + + class HDF5EventSource(EventSource): """ Event source for files in the ctapipe DL1 format. @@ -248,6 +277,10 @@ def subarray(self): def datalevels(self): return get_hdf5_datalevels(self.file_) + @lazyproperty + def atmosphere_density_profile(self) -> AtmosphereDensityProfile: + return read_atmosphere_density_profile(self.file_) + @lazyproperty def obs_ids(self): return list(np.unique(self.file_.root.dl1.event.subarray.trigger.col("obs_id"))) diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py index 4dd57cc20ad..ad7b35c565e 100644 --- a/ctapipe/io/simteleventsource.py +++ b/ctapipe/io/simteleventsource.py @@ -1,17 +1,25 @@ import enum import warnings +from contextlib import nullcontext +from enum import Enum, auto, unique from gzip import GzipFile from io import BufferedReader from pathlib import Path -from typing import Dict +from typing import Dict, Optional, Union import numpy as np from astropy import units as u from astropy.coordinates import Angle, EarthLocation +from astropy.table import Table from astropy.time import Time from eventio.file_types import is_eventio from eventio.simtel.simtelfile import SimTelFile +from ..atmosphere import ( + AtmosphereDensityProfile, + FiveLayerAtmosphereDensityProfile, + TableAtmosphereDensityProfile, +) from ..calib.camera.gainselection import GainSelector from ..containers import ( ArrayEventContainer, @@ -39,6 +47,7 @@ ) from ..coordinates import CameraFrame from ..core import Map +from ..core.provenance import Provenance from ..core.traits import Bool, Float, Undefined, UseEnum, create_class_enum_trait from ..instrument import ( CameraDescription, @@ -61,7 +70,12 @@ from .datalevels import DataLevel from .eventsource import EventSource -__all__ = ["SimTelEventSource"] +__all__ = [ + "SimTelEventSource", + "read_atmosphere_profile_from_simtel", + "AtmosphereProfileKind", + "MirrorClass", +] # Mapping of SimTelArray Calibration trigger types to EventType: # from simtelarray: type Dark (0), pedestal (1), in-lid LED (2) or laser/LED (3+) data. @@ -274,6 +288,113 @@ def apply_simtel_r1_calibration( return r1_waveforms, selected_gain_channel +@unique +class AtmosphereProfileKind(Enum): + """ + choice for which kind of atmosphere density profile to load if more than + one is present""" + + #: Don't load a profile + NONE = auto() + + #: Use a TableAtmosphereDensityProfile + TABLE = auto() + + #: Use a FiveLayerAtmosphereDensityProfile + FIVELAYER = auto() + + #: Try TABLE first, and if it doesn't exist, use FIVELAYER + AUTO = auto() + + +def read_atmosphere_profile_from_simtel( + simtelfile: Union[str, Path, SimTelFile], kind=AtmosphereProfileKind.AUTO +) -> Optional[TableAtmosphereDensityProfile]: + """Read an atmosphere profile from a SimTelArray file as an astropy Table + + Parameters + ---------- + simtelfile: str | SimTelFile + filename of a SimTelArray file containing an atmosphere profile + kind: AtmosphereProfileKind | str + which type of model to load. In AUTO mode, table is tried first, + and if it doesn't exist, fivelayer is used. + + Returns + ------- + Optional[TableAtmosphereDensityProfile]: + Profile read from a table, with interpolation + + """ + if not isinstance(kind, AtmosphereProfileKind): + raise TypeError( + "read_atmosphere_profile_from_simtel: kind should be a AtmosphereProfileKind" + ) + + profiles = [] + + if kind == AtmosphereProfileKind.NONE: + return None + + if isinstance(simtelfile, (str, Path)): + context_manager = SimTelFile(simtelfile) + Provenance().add_input_file( + filename=simtelfile, role="ctapipe.atmosphere.AtmosphereDensityProfile" + ) + + else: + context_manager = nullcontext(simtelfile) + + with context_manager as simtel: + + if ( + not hasattr(simtel, "atmospheric_profiles") + or len(simtel.atmospheric_profiles) == 0 + ): + return [] + + for atmo in simtel.atmospheric_profiles: + + metadata = dict( + observation_level=atmo["obslevel"] * u.cm, + atmosphere_id=atmo["id"], + atmosphere_name=atmo["name"].decode("utf-8"), + atmosphere_height=atmo["htoa"] * u.cm, + ) + + if "altitude_km" in atmo and kind in { + AtmosphereProfileKind.TABLE, + AtmosphereProfileKind.AUTO, + }: + table = Table( + dict( + height=atmo["altitude_km"] * u.km, + density=atmo["rho"] * u.g / u.cm**3, + column_density=atmo["thickness"] * u.g / u.cm**2, + ), + meta=metadata, + ) + profiles.append(TableAtmosphereDensityProfile(table=table)) + + elif ( + kind == AtmosphereProfileKind.FIVELAYER + and "five_layer_atmosphere" in atmo + ): + profiles.append( + FiveLayerAtmosphereDensityProfile.from_array( + atmo["five_layer_atmosphere"], metadata=metadata + ) + ) + + else: + raise ValueError(f"Couldn't load requested profile '{kind}'") + + if profiles: + return profiles[0] + + return None + + class SimTelEventSource(EventSource): """ Read events from a SimTelArray data file (in EventIO format). @@ -348,6 +469,17 @@ class SimTelEventSource(EventSource): ), ).tag(config=True) + atmosphere_profile_choice = UseEnum( + AtmosphereProfileKind, + default_value=AtmosphereProfileKind.AUTO, + help=( + "Which type of atmosphere density profile to load " + "from the file, in case more than one exists. If set " + "to AUTO, TABLE will be attempted first and if missing, " + "FIVELAYER will be loaded." + ), + ).tag(config=True) + def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs): """ EventSource for simtelarray files using the pyeventio library. @@ -390,6 +522,11 @@ def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs): self.gain_selector = GainSelector.from_name( self.gain_selector_type, parent=self ) + + self._atmosphere_density_profile = read_atmosphere_profile_from_simtel( + self.file_, kind=self.atmosphere_profile_choice + ) + self.log.debug(f"Using gain selector {self.gain_selector}") def __exit__(self, exc_type, exc_val, exc_tb): @@ -410,6 +547,10 @@ def datalevels(self): def simulation_config(self) -> Dict[int, SimulationConfigContainer]: return self._simulation_config + @property + def atmosphere_density_profile(self) -> AtmosphereDensityProfile: + return self._atmosphere_density_profile + @property def observation_blocks(self) -> Dict[int, ObservationBlockContainer]: """ diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py index 80a23be0fa8..04cf133d414 100644 --- a/ctapipe/io/tests/test_simteleventsource.py +++ b/ctapipe/io/tests/test_simteleventsource.py @@ -1,3 +1,5 @@ +""" tests of SimTelEventSource """ +# pylint: disable=import-outside-toplevel import copy from itertools import zip_longest from pathlib import Path @@ -13,7 +15,12 @@ from ctapipe.instrument.camera.geometry import UnknownPixelShapeWarning from ctapipe.instrument.optics import ReflectorShape from ctapipe.io import DataLevel -from ctapipe.io.simteleventsource import SimTelEventSource, apply_simtel_r1_calibration +from ctapipe.io.simteleventsource import ( + AtmosphereProfileKind, + SimTelEventSource, + apply_simtel_r1_calibration, + read_atmosphere_profile_from_simtel, +) from ctapipe.utils import get_dataset_path gamma_test_large_path = get_dataset_path("gamma_test_large.simtel.gz") @@ -517,6 +524,45 @@ def test_simtel_no_metadata(monkeypatch): assert all([t.optics.name.startswith("UNKNOWN") for t in subarray.tel.values()]) +def test_load_atmosphere_from_simtel(prod5_gamma_simtel_path): + """ + Load atmosphere from a SimTelEventSource + """ + from ctapipe.atmosphere import ( + FiveLayerAtmosphereDensityProfile, + TableAtmosphereDensityProfile, + ) + + profile = read_atmosphere_profile_from_simtel( + prod5_gamma_simtel_path, kind=AtmosphereProfileKind.AUTO + ) + assert isinstance(profile, TableAtmosphereDensityProfile) + + profile = read_atmosphere_profile_from_simtel( + prod5_gamma_simtel_path, kind=AtmosphereProfileKind.TABLE + ) + assert isinstance(profile, TableAtmosphereDensityProfile) + + profile = read_atmosphere_profile_from_simtel( + prod5_gamma_simtel_path, kind=AtmosphereProfileKind.FIVELAYER + ) + assert isinstance(profile, FiveLayerAtmosphereDensityProfile) + + # old simtel files don't have the profile in them, so a null list should be + # returned + simtel_path_old = get_dataset_path("gamma_test_large.simtel.gz") + profile = read_atmosphere_profile_from_simtel(simtel_path_old) + assert not profile + + +def test_atmosphere_profile(prod5_gamma_simtel_path): + """check that for a file with a profile in it that we get it back""" + from ctapipe.atmosphere import AtmosphereDensityProfile + + with SimTelEventSource(prod5_gamma_simtel_path) as source: + assert isinstance(source.atmosphere_density_profile, AtmosphereDensityProfile) + + @pytest.mark.parametrize("sign", (-1, 1)) def test_float32_pihalf(sign): float32_pihalf = np.float32(sign * np.pi / 2) diff --git a/ctapipe/tests/test_atmosphere.py b/ctapipe/tests/test_atmosphere.py new file mode 100644 index 00000000000..2fc0648488c --- /dev/null +++ b/ctapipe/tests/test_atmosphere.py @@ -0,0 +1,135 @@ +""" +Checks for the atmosphere module +""" + +# pylint: disable=import-outside-toplevel +from unittest.mock import patch + +import astropy.units as u +import numpy as np +import pytest + +import ctapipe.atmosphere as atmo +from ctapipe.utils import get_dataset_path + +SIMTEL_PATH = get_dataset_path( + "gamma_20deg_0deg_run2___cta-prod5-paranal_desert" + "-2147m-Paranal-dark_cone10-100evts.simtel.zst" +) + + +def get_simtel_profile_from_eventsource(): + """get a TableAtmosphereDensityProfile from a simtel file""" + from ctapipe.io import EventSource + + with EventSource(SIMTEL_PATH) as source: + return source.atmosphere_density_profile + + +@pytest.fixture(scope="session", name="table_profile") +def fixture_table_profile(): + """a table profile for testing""" + return get_simtel_profile_from_eventsource() + + +def get_simtel_fivelayer_profile(): + """ + get a sample 3-layer profile + """ + from ctapipe.io.simteleventsource import ( + AtmosphereProfileKind, + read_atmosphere_profile_from_simtel, + ) + + return read_atmosphere_profile_from_simtel( + SIMTEL_PATH, kind=AtmosphereProfileKind.FIVELAYER + ) + + +@pytest.mark.parametrize( + "density_model", + [ + atmo.ExponentialAtmosphereDensityProfile(), + get_simtel_profile_from_eventsource(), + get_simtel_fivelayer_profile(), + ], +) +def test_models(density_model): + """check that a set of model classes work""" + + # test we can convert to correct units + density_model(10 * u.km).to(u.kg / u.m**3) + + # ensure units are properly taken into account + assert np.isclose(density_model(1 * u.km), density_model(1000 * u.m)) + + # check we can also compute the integral + column_density = density_model.integral(10 * u.km) + assert column_density.unit.is_equivalent(u.g / u.cm**2) + + assert np.isclose( + density_model.integral(1 * u.km), density_model.integral(1000 * u.m) + ) + + with patch("matplotlib.pyplot.show"): + density_model.peek() + + +def test_exponential_model(): + """check exponential models""" + + density_model = atmo.ExponentialAtmosphereDensityProfile( + scale_height=10 * u.m, scale_density=0.00125 * u.g / u.cm**3 + ) + assert np.isclose(density_model(1_000_000 * u.km), 0 * u.g / u.cm**3) + assert np.isclose(density_model(0 * u.km), density_model.scale_density) + + +def test_table_model_interpolation(table_profile): + """check that interpolation is reasonable""" + + np.testing.assert_allclose( + table_profile(table_profile.table["height"].to("km")), + table_profile.table["density"].to("g cm-3"), + ) + + # check that fine interpolation up to 100 km : + height_fine = np.linspace(0, 100, 1000) * u.km + assert np.isfinite(table_profile.integral(height_fine)).all() + + +def test_against_reference(): + """ + Test five-layer and table methods against a reference analysis from + SimTelArray. Data from communication with K. Bernloehr. + + See https://github.com/cta-observatory/ctapipe/pull/2000 + """ + from ctapipe.utils import get_table_dataset + + reference_table = get_table_dataset( + "atmosphere_profile_comparison_from_simtelarray" + ) + + fit_reference = np.array( + [ + [0.00 * 100000, -140.508, 1178.05, 994186, 0], + [9.75 * 100000, -18.4377, 1265.08, 708915, 0], + [19.00 * 100000, 0.217565, 1349.22, 636143, 0], + [46.00 * 100000, -0.000201796, 703.745, 721128, 0], + [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0], + ] + ) + + profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference) + + height = reference_table["Altitude_km"].to("km") + + np.testing.assert_allclose( + 1.0 - profile_5(height) / reference_table["rho_5"], 0, atol=1e-5 + ) + np.testing.assert_allclose( + 1.0 - profile_5.line_of_sight_integral(height) / reference_table["thick_5"], + 0, + atol=1e-5, + ) diff --git a/docs/conf.py b/docs/conf.py index 2d507da0c60..bf22faa3c70 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -116,6 +116,7 @@ def setup(app): ("py:obj", "name"), ("py:class", "astropy.coordinates.baseframe.BaseCoordinateFrame"), ("py:class", "astropy.table.table.Table"), + ("py:class", "eventio.simtel.simtelfile.SimTelFile"), ] # The suffix(es) of source filenames. diff --git a/docs/ctapipe_api/atmosphere/index.rst b/docs/ctapipe_api/atmosphere/index.rst new file mode 100644 index 00000000000..78ca82592d0 --- /dev/null +++ b/docs/ctapipe_api/atmosphere/index.rst @@ -0,0 +1,16 @@ +.. _atmosphere: + +================================================== +Atmosphere Models (`~ctapipe.atmosphere`) +================================================== + +.. currentmodule:: ctapipe.atmosphere + +Models of the atmosphere useful for transforming between *column densities* +(atmosphere depths :math:`X`, in units of mass per area) and *heights* (in distance +above sea-level units). + +Reference/API +============= + +.. automodapi:: ctapipe.atmosphere