diff --git a/ctapipe/calib/camera/calibrator.py b/ctapipe/calib/camera/calibrator.py index 5383e0d81ce..f5ea3fad65c 100644 --- a/ctapipe/calib/camera/calibrator.py +++ b/ctapipe/calib/camera/calibrator.py @@ -284,7 +284,7 @@ def __call__(self, event): Parameters ---------- event : container - A `~ctapipe.containers.ArrayEventContainer` event container + A `~ctapipe.containers.SubarrayEventContainer` event container """ for tel_event in event.tel.values(): self.calibrate_tel_event(tel_event) diff --git a/ctapipe/calib/camera/flatfield.py b/ctapipe/calib/camera/flatfield.py index 3a3139ebcf0..d461e69a04a 100644 --- a/ctapipe/calib/camera/flatfield.py +++ b/ctapipe/calib/camera/flatfield.py @@ -108,7 +108,7 @@ def calculate_relative_gain(self, event): Parameters ---------- - event: ctapipe.containers.ArrayEventContainer + event: ctapipe.containers.SubarrayEventContainer Returns: True if the mon.tel[tel_id].flatfield is updated, False otherwise diff --git a/ctapipe/calib/camera/pedestals.py b/ctapipe/calib/camera/pedestals.py index 56a859db62c..4c7ff3385d2 100644 --- a/ctapipe/calib/camera/pedestals.py +++ b/ctapipe/calib/camera/pedestals.py @@ -134,7 +134,7 @@ def calculate_pedestals(self, event): Parameters ---------- - event: ctapipe.containers.ArrayEventContainer + event: ctapipe.containers.SubarrayEventContainer Returns: True if the mon.tel[tel_id].pedestal is updated, False otherwise @@ -203,12 +203,12 @@ def _extract_charge(self, event) -> DL1TelescopeContainer: Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer general event container Returns ------- - DL1CameraContainer + DL1TelescopeContainer """ waveforms = event.tel[self.tel_id].r1.waveform selected_gain_channel = event.tel[self.tel_id].r1.selected_gain_channel diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 4c8f74b4831..ca33a10105f 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -12,49 +12,51 @@ from .core import Container, Field, Map __all__ = [ - "SubarrayEventContainer", + "BaseHillasParametersContainer", + "BaseTimingParametersContainer", + "CameraHillasParametersContainer", + "CameraTimingParametersContainer", "ConcentrationContainer", + "CoreParametersContainer", "DL0TelescopeContainer", "DL1CameraCalibrationContainer", "DL1TelescopeContainer", "DL2SubarrayContainer", - "TelescopeCalibrationContainer", - "ArrayEventIndexContainer", + "DispContainer", "EventType", "FlatFieldContainer", "HillasParametersContainer", - "CoreParametersContainer", "ImageParametersContainer", + "IntensityStatisticsContainer", "LeakageContainer", - "TelescopeMonitoringContainer", "MorphologyContainer", - "BaseHillasParametersContainer", - "CameraHillasParametersContainer", - "CameraTimingParametersContainer", + "ObservationBlockContainer", + "ObservationBlockState", + "ObservingMode", "ParticleClassificationContainer", + "PeakTimeStatisticsContainer", "PedestalContainer", "PixelStatusContainer", "R0TelescopeContainer", "R1TelescopeContainer", + "ReconstructedShowerContainer", "ReconstructedEnergyContainer", "ReconstructedGeometryContainer", - "DispContainer", - "SimulationTelescopeContainer", + "SchedulingBlockContainer", "SimulatedShowerContainer", "SimulatedShowerDistribution", "SimulationConfigContainer", + "SimulationTelescopeContainer", + "StatisticsContainer", + "SubarrayEventContainer", + "SubarrayEventIndexContainer", + "SubarrayTriggerContainer", + "TelescopeCalibrationContainer", + "TelescopeEventContainer", "TelescopeEventIndexContainer", - "BaseTimingParametersContainer", + "TelescopeMonitoringContainer", "TimingParametersContainer", - "SubarrayTriggerContainer", "WaveformCalibrationContainer", - "StatisticsContainer", - "IntensityStatisticsContainer", - "PeakTimeStatisticsContainer", - "SchedulingBlockContainer", - "ObservationBlockContainer", - "ObservingMode", - "ObservationBlockState", ] @@ -182,7 +184,7 @@ class TelescopeEventType(enum.Enum): UNKNOWN = 255 -class ArrayEventIndexContainer(Container): +class SubarrayEventIndexContainer(Container): """index columns to include in event lists, common to all data levels""" default_prefix = "" # don't want to prefix these @@ -609,7 +611,7 @@ class SimulatedShowerContainer(Container): class SimulationTelescopeContainer(Container): """ - True images and parameters derived from them, analgous to the `DL1CameraContainer` + True images and parameters derived from them, analgous to the `DL1TelescopeContainer` but for simulated data. """ @@ -1221,7 +1223,7 @@ class SubarrayEventContainer(Container): count = Field(0, description="number of events processed") index = Field( - default_factory=ArrayEventIndexContainer, + default_factory=SubarrayEventIndexContainer, description="event indexing information", ) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index fd9f7c2ab68..f84289860d1 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -14,6 +14,7 @@ IntensityStatisticsContainer, PeakTimeStatisticsContainer, SubarrayEventContainer, + TelescopeEventContainer, TimingParametersContainer, ) from ..core import QualityQuery, TelescopeComponent @@ -210,7 +211,7 @@ def _parameterize_image( # parameterization return default - def _process_telescope_event(self, tel_event): + def _process_telescope_event(self, tel_event: TelescopeEventContainer): """ Loop over telescopes and process the calibrated images into parameters """ diff --git a/ctapipe/image/muon/processor.py b/ctapipe/image/muon/processor.py index 1fac1270f6f..5eb11ac01de 100644 --- a/ctapipe/image/muon/processor.py +++ b/ctapipe/image/muon/processor.py @@ -124,7 +124,7 @@ def _process_telescope_event(self, tel_event): Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer Collection of all event information tel_id: int Telescope ID of the instrument that has measured the image diff --git a/ctapipe/io/eventsource.py b/ctapipe/io/eventsource.py index 10e3d4bf1e5..1129a7265a7 100644 --- a/ctapipe/io/eventsource.py +++ b/ctapipe/io/eventsource.py @@ -28,7 +28,7 @@ class EventSource(Component): """ Parent class for EventSources. - EventSources read input files and generate `~ctapipe.containers.ArrayEventContainer` + EventSources read input files and generate `~ctapipe.containers.SubarrayEventContainer` instances when iterated over. A new EventSource should be created for each type of event file read @@ -72,7 +72,7 @@ class EventSource(Component): 0 1 - **NOTE**: EventSource implementations should not reuse the same ArrayEventContainer, + **NOTE**: EventSource implementations should not reuse the same SubarrayEventContainer, as these are mutable and may lead to errors when analyzing multiple events. diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py index f2088cdae89..23b9035ed3a 100644 --- a/ctapipe/io/hdf5eventsource.py +++ b/ctapipe/io/hdf5eventsource.py @@ -12,7 +12,6 @@ from ctapipe.instrument.optics import FocalLengthKind from ..containers import ( - ArrayEventIndexContainer, CameraHillasParametersContainer, CameraTimingParametersContainer, ConcentrationContainer, @@ -41,6 +40,7 @@ SimulationSubarrayContainer, SimulationTelescopeContainer, SubarrayEventContainer, + SubarrayEventIndexContainer, SubarrayTriggerContainer, TelescopeEventContainer, TelescopeEventIndexContainer, @@ -138,8 +138,8 @@ class HDF5EventSource(EventSource): specifying the file to be read. Looping over the EventSource yields events from the _generate_events - method. An event equals an ArrayEventContainer instance. - See ctapipe.containers.ArrayEventContainer for details. + method. An event equals an SubarrayEventContainer instance. + See ctapipe.containers.SubarrayEventContainer for details. Attributes ---------- @@ -393,7 +393,7 @@ def _is_hillas_in_camera_frame(self): def _generator(self): """ - Yield ArrayEventContainer to iterate through events. + Yield SubarrayEventContainer to iterate through events. """ self.reader = HDF5TableReader(self.file_) @@ -562,7 +562,7 @@ def _generator(self): # Setup iterators for the array events events = HDF5TableReader(self.file_).read( "/dl0/event/subarray/trigger", - [SubarrayTriggerContainer, ArrayEventIndexContainer], + [SubarrayTriggerContainer, SubarrayEventIndexContainer], ignore_columns={"tel"}, ) telescope_trigger_reader = HDF5TableReader(self.file_).read( diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py index a7153278cb1..ddc1134785d 100644 --- a/ctapipe/io/simteleventsource.py +++ b/ctapipe/io/simteleventsource.py @@ -22,7 +22,6 @@ ) from ..calib.camera.gainselection import GainSelector from ..containers import ( - ArrayEventIndexContainer, CoordinateFrameType, DL0SubarrayContainer, DL0TelescopeContainer, @@ -42,6 +41,7 @@ SimulationSubarrayContainer, SimulationTelescopeContainer, SubarrayEventContainer, + SubarrayEventIndexContainer, SubarrayPointingContainer, SubarrayTriggerContainer, TelescopeCalibrationContainer, @@ -740,7 +740,7 @@ def _generate_events(self): simulation=SimulationSubarrayContainer(shower=shower), dl0=DL0SubarrayContainer(trigger=array_trigger), pointing=self._fill_array_pointing(), - index=ArrayEventIndexContainer(obs_id=obs_id, event_id=event_id), + index=SubarrayEventIndexContainer(obs_id=obs_id, event_id=event_id), ) event.meta["origin"] = "hessio" event.meta["input_url"] = self.input_url diff --git a/ctapipe/io/toymodel.py b/ctapipe/io/toymodel.py index e5abfceb1e4..bfc150552da 100644 --- a/ctapipe/io/toymodel.py +++ b/ctapipe/io/toymodel.py @@ -9,11 +9,11 @@ import numpy as np from ..containers import ( - ArrayEventIndexContainer, DL1TelescopeContainer, ObservationBlockContainer, SchedulingBlockContainer, SubarrayEventContainer, + SubarrayEventIndexContainer, ) from ..core import TelescopeComponent, traits from ..image import toymodel @@ -99,7 +99,7 @@ def _generator(self): def generate_event(self): event = SubarrayEventContainer( - index=ArrayEventIndexContainer(obs_id=1, event_id=self.event_id), + index=SubarrayEventIndexContainer(obs_id=1, event_id=self.event_id), dl0=None, dl2=None, simulation=None, diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py index 793207c5f1f..7c0a9822b8d 100644 --- a/ctapipe/reco/hillas_intersection.py +++ b/ctapipe/reco/hillas_intersection.py @@ -93,7 +93,7 @@ def __call__(self, event): Parameters ---------- - event : `~ctapipe.containers.ArrayEventContainer` + event : `~ctapipe.containers.SubarrayEventContainer` The event, needs to have dl1 parameters. Will be filled with the corresponding dl2 containers, reconstructed stereo geometry and telescope-wise impact position. diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py index 588029157d8..55a2e32ed5b 100644 --- a/ctapipe/reco/hillas_reconstructor.py +++ b/ctapipe/reco/hillas_reconstructor.py @@ -134,7 +134,7 @@ def __call__(self, event): Parameters ---------- - event : `~ctapipe.containers.ArrayEventContainer` + event : `~ctapipe.containers.SubarrayEventContainer` The event, needs to have dl1 parameters. Will be filled with the corresponding dl2 containers, reconstructed stereo geometry and telescope-wise impact position. diff --git a/ctapipe/reco/preprocessing.py b/ctapipe/reco/preprocessing.py index d7fd8f5f958..2822a164783 100644 --- a/ctapipe/reco/preprocessing.py +++ b/ctapipe/reco/preprocessing.py @@ -75,7 +75,7 @@ def collect_features( Parameters ---------- - event : ArrayEventContainer + event : SubarrayEventContainer The event container from which to collect the features tel_id : int The telscope id for which to collect the features diff --git a/ctapipe/reco/reconstructor.py b/ctapipe/reco/reconstructor.py index 1803a13b4c4..af32ae73020 100644 --- a/ctapipe/reco/reconstructor.py +++ b/ctapipe/reco/reconstructor.py @@ -87,7 +87,7 @@ def __call__(self, event: SubarrayEventContainer): Parameters ---------- - event : `ctapipe.containers.ArrayEventContainer` + event : `ctapipe.containers.SubarrayEventContainer` The event, needs to have dl1 parameters. Will be filled with the corresponding dl2 containers, reconstructed stereo geometry and telescope-wise impact position. diff --git a/ctapipe/reco/shower_processor.py b/ctapipe/reco/shower_processor.py index 86e3b267d37..b9e0b64e416 100644 --- a/ctapipe/reco/shower_processor.py +++ b/ctapipe/reco/shower_processor.py @@ -75,7 +75,7 @@ def __call__(self, event: SubarrayEventContainer): Parameters ---------- - event : ctapipe.containers.ArrayEventContainer + event : ctapipe.containers.SubarrayEventContainer Top-level container for all event information. """ for reconstructor in self.reconstructors: diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py index 77afdb3e632..0f1be197270 100644 --- a/ctapipe/reco/sklearn.py +++ b/ctapipe/reco/sklearn.py @@ -161,7 +161,7 @@ def __call__(self, event: SubarrayEventContainer) -> None: Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer """ @abstractmethod @@ -677,7 +677,7 @@ def __call__(self, event: SubarrayEventContainer) -> None: Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer """ for tel_id, tel_event in event.tel.items(): table = collect_features(event, tel_id, self.instrument_table) diff --git a/ctapipe/tools/tests/test_apply_models.py b/ctapipe/tools/tests/test_apply_models.py index f9f95a69739..13090595c40 100644 --- a/ctapipe/tools/tests/test_apply_models.py +++ b/ctapipe/tools/tests/test_apply_models.py @@ -2,10 +2,10 @@ import pytest from ctapipe.containers import ( - ArrayEventIndexContainer, ParticleClassificationContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, + SubarrayEventIndexContainer, ) from ctapipe.core import run_tool from ctapipe.core.tool import ToolConfigurationError @@ -42,7 +42,8 @@ def test_apply_energy_regressor( table = read_table(output_path, f"/dl2/event/subarray/energy/{prefix}") for col in "obs_id", "event_id": assert ( - table[col].description == ArrayEventIndexContainer.fields[col].description + table[col].description + == SubarrayEventIndexContainer.fields[col].description ) for name, field in ReconstructedEnergyContainer.fields.items(): diff --git a/docs/ctapipe_api/containers/index.rst b/docs/ctapipe_api/containers/index.rst index 5be56a863e5..15cbdfc15f1 100644 --- a/docs/ctapipe_api/containers/index.rst +++ b/docs/ctapipe_api/containers/index.rst @@ -13,7 +13,7 @@ The `ctapipe.containers` module contains the data model definition of all ctapipe `~ctapipe.core.Container` classes, which provide the container definitions for all ctapipe data. -The base Container for an event is in `ctapipe.containers.ArrayEventContainer`. +The base Container for an event is in `ctapipe.containers.SubarrayEventContainer`. Reference/API diff --git a/docs/ctapipe_api/io/index.rst b/docs/ctapipe_api/io/index.rst index 82f18d7c908..73bdbfe7b9d 100644 --- a/docs/ctapipe_api/io/index.rst +++ b/docs/ctapipe_api/io/index.rst @@ -131,7 +131,7 @@ Writing Output Files: ===================== The `DataWriter` Component allows one to write a series of events (stored in -`ctapipe.containers.ArrayEventContainer`) to a standardized HDF5 format file +`ctapipe.containers.SubarrayEventContainer`) to a standardized HDF5 format file following the data model (see :ref:`datamodels`). This includes all related datasets such as the instrument and simulation configuration information, simulated shower and image information, observed images and parameters and reconstruction diff --git a/docs/data_models/dl1.rst b/docs/data_models/dl1.rst index f792d495074..126f14acd34 100644 --- a/docs/data_models/dl1.rst +++ b/docs/data_models/dl1.rst @@ -35,16 +35,16 @@ The following datasets will be written to the group ``/dl1/event/`` in the outp - (group) * - ``subarray/trigger`` - subarray trigger information - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.TriggerContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.SubarrayTriggerContainer` * - ``telescope/`` - Per-telescope Per-event information - (group) * - ``telescope/parameters/tel_{TEL_ID:03d}`` - tables of image parameters (one per telescope) - - :py:class:`~ctapipe.containers.TelEventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` + - :py:class:`~ctapipe.containers.TelescopeEventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` * - ``telescope/images/tel_{TEL_ID:03d}`` - tables of telescope images (one per telescope) - - :py:class:`~ctapipe.containers.TelEventIndexContainer` +, :py:class:`~ctapipe.containers.DL1CameraContainer` + - :py:class:`~ctapipe.containers.TelescopeEventIndexContainer` +, :py:class:`~ctapipe.containers.DL1TelescopeContainer` -------------- DL2 Data Model @@ -65,26 +65,26 @@ output file, where ```` is the identifier of the algorithm (e.g. - Contents * - /geometry - shower geometry reconstruction - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedGeometryContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedGeometryContainer` * - /energy - shower energy reconstruction - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedEnergyContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedEnergyContainer` * - /classification - shower classification parameters - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ParticleClassificationContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer`, :py:class:`~ctapipe.containers.ParticleClassificationContainer` --------------------- Simulation Data Model --------------------- * - ``/simulation/event/subarray/shower`` - true shower parameters from Monte-Carlo simulation - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedShowerContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedShowerContainer` * - ``/simulation/event/telescope/images/tel_{TEL_ID:03d}`` - simulated camera images - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedCameraContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedCameraContainer` * - ``/simulation/event/telescope/parameters/tel_{TEL_ID:03d}`` - Parameters derived form the simulated camera images - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` * - ``/simulation/service/shower_distribution`` - simulated shower distribution histograms - :py:class:`~ctapipe.containers.SimulatedShowerDistribution` diff --git a/docs/examples/camera_display.ipynb b/docs/examples/camera_display.ipynb index 3d1274ee515..4b7adff21b1 100644 --- a/docs/examples/camera_display.ipynb +++ b/docs/examples/camera_display.ipynb @@ -13,7 +13,7 @@ "metadata": {}, "outputs": [], "source": [ - "import astropy.coordinates as c\n", + "from astropy.coordinates import SkyCoord\n", "import astropy.units as u\n", "import matplotlib.pylab as plt\n", "import numpy as np\n", @@ -50,7 +50,9 @@ ")\n", "\n", "image, sig, bg = model.generate_image(geom, intensity=1500, nsb_level_pe=10)\n", - "mask = tailcuts_clean(geom, image, picture_thresh=15, boundary_thresh=5)" + "mask = tailcuts_clean(\n", + " geom, image, picture_thresh=15, boundary_thresh=5, min_number_picture_neighbors=2\n", + ")" ] }, { @@ -133,7 +135,7 @@ "outputs": [], "source": [ "geom_camframe = geom\n", - "geom = geom_camframe.transform_to(EngineeringCameraFrame())" + "geom = geom_camframe.transform_to(TelescopeFrame())" ] }, { @@ -158,7 +160,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(1, 3, figsize=(15, 4))\n", - "for ii, cmap in enumerate([\"PuOr_r\", \"rainbow\", \"twilight\"]):\n", + "for ii, cmap in enumerate([\"inferno\", \"viridis\", \"cividis\"]):\n", " disp = CameraDisplay(geom, image=image, ax=ax[ii], title=cmap)\n", " disp.add_colorbar()\n", " disp.cmap = cmap" @@ -229,8 +231,11 @@ "fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n", "norms = [\"lin\", \"log\", \"symlog\", PowerNorm(0.5)]\n", "\n", + "# show invalid values in gray (e.g. values < 0 for log color map)\n", + "cmap = plt.get_cmap(\"inferno\").with_extremes(bad=\"gray\")\n", + "\n", "for norm, ax in zip(norms, axes.flatten()):\n", - " disp = CameraDisplay(geom, image=image, ax=ax)\n", + " disp = CameraDisplay(geom, image=image, ax=ax, cmap=cmap)\n", " disp.norm = norm\n", " disp.add_colorbar()\n", " ax.set_title(str(norm))\n", @@ -265,15 +270,15 @@ "disp = CameraDisplay(\n", " geom, image=image, cmap=\"gray\", ax=ax[0], title=\"Image mask in green\"\n", ")\n", - "disp.highlight_pixels(mask, alpha=0.8, linewidth=2, color=\"green\")\n", + "disp.highlight_pixels(mask, alpha=0.8, linewidth=1, color=\"green\")\n", "\n", "disp = CameraDisplay(\n", " geom, image=image, cmap=\"gray\", ax=ax[1], title=\"Image mask in green (zoom)\"\n", ")\n", "disp.highlight_pixels(mask, alpha=1, linewidth=3, color=\"green\")\n", "\n", - "ax[1].set_ylim(-0.5, 0.5)\n", - "ax[1].set_xlim(-0.5, 0.5)" + "ax[1].set_xlim(-1, 1.5)\n", + "ax[1].set_ylim(-0.2, 1.8)" ] }, { @@ -296,14 +301,14 @@ "metadata": {}, "outputs": [], "source": [ - "clean_image = image.copy()\n", - "clean_image[~mask] = 0\n", - "hillas = hillas_parameters(geom, clean_image)\n", + "hillas = hillas_parameters(geom[mask], image[mask])\n", "\n", "plt.figure(figsize=(6, 6))\n", "disp = CameraDisplay(geom, image=image, cmap=\"gray_r\")\n", "disp.highlight_pixels(mask, alpha=0.5, color=\"dodgerblue\")\n", - "disp.overlay_moments(hillas, color=\"red\", linewidth=3, with_label=False)" + "\n", + "disp.overlay_moments(hillas, color=\"red\", linewidth=1, with_label=False)\n", + "disp.overlay_moments(hillas, color=\"red\", linewidth=2, n_sigma=2, keep_old=True)" ] }, { @@ -317,6 +322,18 @@ "Note that the parameter `keep_old` is False by default, meaning adding a new point will clear the previous ones (useful for animations, but perhaps unexpected for a static plot). Set it to `True` to plot multiple markers." ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# for transformations to work, we need the telescope pointing attached to the frame\n", + "telescope_pointing = SkyCoord(alt=70 * u.deg, az=20 * u.deg, frame=\"altaz\")\n", + "\n", + "geom.frame = geom.frame.replicate_without_data(telescope_pointing=telescope_pointing)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -324,10 +341,14 @@ "outputs": [], "source": [ "plt.figure(figsize=(6, 6))\n", + "\n", + "\n", "disp = CameraDisplay(geom, image=image, cmap=\"gray_r\")\n", "\n", - "coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom.frame)\n", - "coord_in_another_frame = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=CameraFrame())\n", + "coord = SkyCoord(fov_lon=0.5 * u.deg, fov_lat=0.7 * u.deg, frame=geom.frame)\n", + "\n", + "coord_in_another_frame = SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom_camframe.frame)\n", + "\n", "disp.overlay_coordinate(coord, markersize=20, marker=\"*\")\n", "disp.overlay_coordinate(\n", " coord_in_another_frame, markersize=20, marker=\"*\", keep_old=True\n", @@ -438,9 +459,9 @@ ") as source:\n", " event = next(iter(source))\n", "\n", - "tel_id = list(event.r0.tel.keys())[0]\n", + "tel_id = list(event.tel.keys())[0]\n", "geom = source.subarray.tel[tel_id].camera.geometry\n", - "waveform = event.r0.tel[tel_id].waveform\n", + "waveform = event.tel[tel_id].r0.waveform\n", "n_chan, n_pix, n_samp = waveform.shape" ] }, @@ -534,25 +555,26 @@ "pixel 5 has value 2.44\n", "```" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, "language_info": { "codemirror_mode": { - "name": "ipython" + "name": "ipython", + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", - "nbconvert_exporter": "python" + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.16" } }, "nbformat": 4, diff --git a/docs/tutorials/ctapipe_overview.ipynb b/docs/tutorials/ctapipe_overview.ipynb index f7b020fb610..89f30b399e7 100644 --- a/docs/tutorials/ctapipe_overview.ipynb +++ b/docs/tutorials/ctapipe_overview.ipynb @@ -20,14 +20,14 @@ "source": [ "
\n", "\n", - "\"ctapipe\"/\n", + "\"ctapipe\"/\n", "\n", "\n", - "

Initially presented @ LST Analysis Bootcamp

\n", + "

Initially presented @ LST Analysis Bootcamp, update since to be compatible with latest ctapipe version

\n", "\n", "

Padova, 26.11.2018

\n", "\n", - "

Maximilian Nöthe (@maxnoe) & Kai A. Brügge (@mackaiver)

\n", + "

Maximilian Linhoff (@maxnoe) & Kai A. Brügge (@mackaiver)

\n", "\n", "
" ] @@ -126,8 +126,9 @@ "* ctapipe is developed as Open Source Software (BSD 3-Clause License) at \n", "\n", "* We use the \"Github-Workflow\": \n", - " * Few people (e.g. @kosack, @maxnoe) have write access to the main repository\n", - " * Contributors fork the main repository and work on branches\n", + " * Few people – the maintainers (@kosack, @maxnoe) – have write access to the main repository\n", + " * Regular contributors are part of the cta-dev team of the cta-observatory organization and work on branches\n", + " * External contributors can fork the main repository and also work on branches\n", " * Pull Requests are merged after Code Review and automatic execution of the test suite\n", "\n", "* Early developement stage ⇒ backwards-incompatible API changes might and will happen \n" @@ -149,8 +150,7 @@ "* Camera and Array plotting\n", "* Coordinate frames and transformations \n", "* Stereo-reconstruction using line intersections\n", - " \n", - " " + "* Machine-Learning Reconstruction of energy, particle type (gammaness) and mono direction using disp" ] }, { @@ -164,7 +164,7 @@ "source": [ "### What's still missing?\n", "\n", - "* Good integration with machine learning techniques\n", + "* Good integration with more advanced machine learning techniques\n", "* IRF calculation \n", "* Documentation, e.g. formal definitions of coordinate frames " ] @@ -199,7 +199,16 @@ } }, "source": [ - "## A simple hillas analysis" + "## A simple hillas analysis\n", + "\n", + "This section shows how to use functions and classes in ctapipe to do a \"custom\" analysis.\n", + "\n", + "ctapipe has command line tools to perform the \"standard\" analysis:\n", + "\n", + "- `ctapipe-process` to perform event-wise processing from DL0 (or R1 or even R0 for simtel array) to DL2 and store results in HDF5 files\n", + "- `ctapipe-merge` to merge the resulting HDF5 files together, e.g. to create datasets for training machine learning models\n", + "- `ctapipe-train-{energy-regressor,particle-classifier,disp-regressor}` to train machine learning models\n", + "- `ctapipe-apply-models` to apply those models to data sets" ] }, { @@ -242,7 +251,7 @@ "for event in source:\n", " print(\n", " \"Id: {}, E = {:1.3f}, Telescopes: {}\".format(\n", - " event.count, event.simulation.shower.energy, len(event.r0.tel)\n", + " event.count, event.simulation.shower.energy, len(event.tel)\n", " )\n", " )" ] @@ -253,7 +262,7 @@ "pycharm": {} }, "source": [ - "Each event is a `DataContainer` holding several `Field`s of data, which can be containers or just numbers.\n", + "Each event is a `SubarrayEventContainer` holding several `Field`s of data, which can be containers or just numbers.\n", "Let's look a one event:" ] }, @@ -268,15 +277,36 @@ "event" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All telescope-specific data is stored in the `TelescopeEventContainer` and the `SubarrayEventContainer` has a `Map` of telescope id to the telescope event:" + ] + }, { "cell_type": "code", "execution_count": null, - "metadata": { - "pycharm": {} - }, + "metadata": {}, "outputs": [], "source": [ - "source.subarray.camera_types" + "event.tel.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "event.tel[3]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The source also provides information about the subarray:" ] }, { @@ -287,7 +317,25 @@ }, "outputs": [], "source": [ - "len(event.r0.tel), len(event.r1.tel)" + "source.subarray.peek()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source.subarray.tel[3]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source.subarray.camera_types" ] }, { @@ -298,7 +346,9 @@ "source": [ "### Data calibration\n", "\n", - "The `CameraCalibrator` calibrates the event (obtaining the `dl1` images)." + "The `CameraCalibrator` calibrates the event from R1 to DL0 (performing data volume reduction) and from DL0 to DL1 (image extraction).\n", + "\n", + "The first step will be done by ACADA for actual observed data but for now is also implemented in ctapipe for prototyping and applying it to simulated data." ] }, { @@ -344,7 +394,7 @@ }, "outputs": [], "source": [ - "event.dl1.tel.keys()" + "event.tel.keys()" ] }, { @@ -355,7 +405,8 @@ }, "outputs": [], "source": [ - "tel_id = 130" + "tel_id = 130\n", + "tel_event = event.tel[tel_id]" ] }, { @@ -367,9 +418,19 @@ "outputs": [], "source": [ "geometry = source.subarray.tel[tel_id].camera.geometry\n", - "dl1 = event.dl1.tel[tel_id]\n", - "\n", - "geometry, dl1" + "geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "pycharm": {} + }, + "outputs": [], + "source": [ + "dl1 = tel_event.dl1\n", + "dl1" ] }, { @@ -395,8 +456,6 @@ "\n", "display = CameraDisplay(geometry)\n", "\n", - "# right now, there might be one image per gain channel.\n", - "# This will change as soon as\n", "display.image = dl1.image\n", "display.add_colorbar()" ] @@ -418,7 +477,7 @@ }, "outputs": [], "source": [ - "from ctapipe.image.cleaning import tailcuts_clean" + "from ctapipe.image.cleaning import TailcutsImageCleaner" ] }, { @@ -429,13 +488,24 @@ }, "outputs": [], "source": [ - "# unoptimized cleaning levels\n", - "cleaning_level = {\n", - " \"CHEC\": (2, 4, 2),\n", - " \"LSTCam\": (3.5, 7, 2),\n", - " \"FlashCam\": (3.5, 7, 2),\n", - " \"NectarCam\": (4, 8, 2),\n", - "}" + "# configure cleaning thresholds per telescope type:\n", + "\n", + "tailcuts_image_cleaner = TailcutsImageCleaner(\n", + " source.subarray,\n", + " picture_threshold_pe=[\n", + " (\"type\", \"LST_*\", 7.0),\n", + " (\"type\", \"*FlashCam\", 7.0),\n", + " (\"type\", \"*NectarCam\", 8.0),\n", + " (\"type\", \"SST_*\", 4.0),\n", + " ],\n", + " boundary_threshold_pe=[\n", + " (\"type\", \"LST_*\", 3.5),\n", + " (\"type\", \"*FlashCam\", 3.5),\n", + " (\"type\", \"*NectarCam\", 4.0),\n", + " (\"type\", \"SST_*\", 2.0),\n", + " ],\n", + " min_picture_neighbors=2,\n", + ")" ] }, { @@ -446,14 +516,10 @@ }, "outputs": [], "source": [ - "boundary, picture, min_neighbors = cleaning_level[geometry.name]\n", - "\n", - "clean = tailcuts_clean(\n", - " geometry,\n", + "clean = tailcuts_image_cleaner(\n", + " tel_id,\n", " dl1.image,\n", - " boundary_thresh=boundary,\n", - " picture_thresh=picture,\n", - " min_number_picture_neighbors=min_neighbors,\n", + " dl1.peak_time,\n", ")" ] }, @@ -540,7 +606,7 @@ "display.image = cleaned\n", "display.add_colorbar()\n", "\n", - "display.overlay_moments(hillas, color=\"xkcd:red\")" + "display.overlay_moments(hillas, color=\"xkcd:light blue\", n_sigma=2)" ] }, { @@ -710,9 +776,11 @@ "\n", " for event in source:\n", " energy = event.simulation.shower.energy\n", - " n_telescopes_r1 = len(event.r1.tel)\n", + " n_telescopes = len(event.dl0.trigger.tels_with_trigger)\n", " event_id = event.index.event_id\n", - " print(f\"Id: {event_id}, E = {energy:1.3f}, Telescopes (R1): {n_telescopes_r1}\")\n", + " print(\n", + " f\"Id: {event_id}, E = {energy:1.3f}, Telescopes (Triggered): {n_telescopes}\"\n", + " )\n", "\n", " calibrator(event)\n", " image_processor(event)\n", @@ -768,7 +836,7 @@ ")\n", "\n", "plt.hist(theta.to_value(u.deg) ** 2, bins=25, range=[0, 0.3])\n", - "plt.xlabel(r\"$\\theta² / deg²$\")\n", + "plt.xlabel(r\"$\\theta^2 \\,/\\, \\mathrm{deg}^2$\")\n", "None" ] }, @@ -795,11 +863,13 @@ "angle_offset = plotting_event.pointing.azimuth\n", "\n", "plotting_hillas = {\n", - " tel_id: dl1.parameters.hillas for tel_id, dl1 in plotting_event.dl1.tel.items()\n", + " tel_id: tel_event.dl1.parameters.hillas\n", + " for tel_id, tel_event in plotting_event.tel.items()\n", "}\n", "\n", "plotting_core = {\n", - " tel_id: dl1.parameters.core.psi for tel_id, dl1 in plotting_event.dl1.tel.items()\n", + " tel_id: tel_event.dl1.parameters.core.psi\n", + " for tel_id, tel_event in plotting_event.tel.items()\n", "}\n", "\n", "\n", @@ -907,7 +977,7 @@ "source": [ "from ctapipe.image import toymodel\n", "from ctapipe.instrument import SubarrayDescription\n", - "\n", + "from ctapipe.image.cleaning import tailcuts_clean\n", "\n", "geometry = loader.subarray.tel[1].camera.geometry" ] @@ -1136,14 +1206,14 @@ "pycharm": {} }, "source": [ - "**A lot less code, and a factor 3 speed improvement**" + "**A lot less code, and a factor 5 speed improvement**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Finally, current ctapipe implementation is using numba:" + "Finally, current ctapipe implementation is using numba, directly operating on the sparse array representation of the neighbor matrix:" ] }, { @@ -1172,7 +1242,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.13" + "version": "3.9.16" } }, "nbformat": 4, diff --git a/docs/tutorials/raw_data_exploration.ipynb b/docs/tutorials/raw_data_exploration.ipynb index 4a5595c4b28..ab3534b4aaf 100644 --- a/docs/tutorials/raw_data_exploration.ipynb +++ b/docs/tutorials/raw_data_exploration.ipynb @@ -76,7 +76,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "the event is just a class with a bunch of data items in it. You can see a more compact represntation via:" + "the event is just a class with a bunch of data items in it. You can see a more compact representation via:" ] }, { @@ -85,23 +85,32 @@ "metadata": {}, "outputs": [], "source": [ - "event.r0" + "event" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "printing the event structure, will currently print the value all items under it (so you get a lot of output if you print a high-level container):" + "printing the event structure, will currently print all items under it (so you get a lot of output if you print a high-level container):" ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ - "print(event.simulation.shower)" + "print(event.simulation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All telescope wise data is in the `.tel` `Map`, which maps telescope ids to `TelescopeEventContainer`:" ] }, { @@ -110,7 +119,7 @@ "metadata": {}, "outputs": [], "source": [ - "print(event.r0.tel.keys())" + "print(event.tel.keys())" ] }, { @@ -127,7 +136,7 @@ "outputs": [], "source": [ "event = next(event_iterator)\n", - "print(event.r0.tel.keys())" + "print(event.tel.keys())" ] }, { @@ -140,12 +149,14 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ - "teldata = event.r0.tel[26]\n", - "print(teldata)\n", - "teldata" + "tel_event = event.tel[26]\n", + "print(repr(tel_event))\n", + "print(tel_event)" ] }, { @@ -170,7 +181,7 @@ "metadata": {}, "outputs": [], "source": [ - "event.simulation.shower.energy.to('GeV')" + "event.simulation.shower.energy.to(\"GeV\")" ] }, { @@ -179,7 +190,7 @@ "metadata": {}, "outputs": [], "source": [ - "event.simulation.shower.energy.to('J')" + "event.simulation.shower.energy.to(\"J\")" ] }, { @@ -205,7 +216,7 @@ "metadata": {}, "source": [ "## Look for signal pixels in a camera\n", - "again, `event.r0.tel[x]` contains a data structure for the telescope data, with some fields like `waveform`.\n", + "again, `tel_event.r0` contains a data structure for the R0 telescope data, with some fields like `waveform`.\n", "\n", "Let's make a 2D plot of the sample data (sample vs pixel), so we can see if we see which pixels contain Cherenkov light signals:" ] @@ -216,7 +227,7 @@ "metadata": {}, "outputs": [], "source": [ - "plt.pcolormesh(teldata.waveform[0]) # note the [0] is for channel 0\n", + "plt.pcolormesh(tel_event.r0.waveform[0]) # note the [0] is for channel 0\n", "plt.colorbar()\n", "plt.xlabel(\"sample number\")\n", "plt.ylabel(\"Pixel_id\")" @@ -235,12 +246,14 @@ "metadata": {}, "outputs": [], "source": [ - "plt.pcolormesh(teldata.waveform[0])\n", + "plt.pcolormesh(tel_event.r0.waveform[0])\n", "plt.colorbar()\n", - "plt.ylim(700,750)\n", + "plt.ylim(700, 750)\n", "plt.xlabel(\"sample number\")\n", "plt.ylabel(\"pixel_id\")\n", - "print(\"waveform[0] is an array of shape (N_pix,N_slice) =\",teldata.waveform[0].shape)" + "print(\n", + " \"waveform[0] is an array of shape (N_pix,N_slice) =\", tel_event.r0.waveform[0].shape\n", + ")" ] }, { @@ -258,8 +271,8 @@ "metadata": {}, "outputs": [], "source": [ - "trace = teldata.waveform[0][719] \n", - "plt.plot(trace, drawstyle='steps')" + "trace = tel_event.r0.waveform[0, 719]\n", + "plt.plot(trace, drawstyle=\"steps\")" ] }, { @@ -277,8 +290,8 @@ "metadata": {}, "outputs": [], "source": [ - "for pix_id in range(718,723):\n", - " plt.plot(teldata.waveform[0][pix_id], label=\"pix {}\".format(pix_id), drawstyle='steps')\n", + "for pix_id in range(718, 723):\n", + " plt.plot(tel_event.r0.waveform[0, pix_id], label=f\"pix {pix_id}\", drawstyle=\"steps\")\n", "plt.legend()" ] }, @@ -300,10 +313,10 @@ "metadata": {}, "outputs": [], "source": [ - "for pix_id in range(718,723):\n", - " plt.plot(teldata.waveform[0][pix_id],'+-')\n", - "plt.fill_betweenx([0,1600],19,24,color='red',alpha=0.3, label='Ped window')\n", - "plt.fill_betweenx([0,1600],5,9,color='green',alpha=0.3, label='Signal window')\n", + "for pix_id in range(718, 723):\n", + " plt.plot(tel_event.r0.waveform[0, pix_id], \"+-\")\n", + "plt.fill_betweenx([0, 1600], 19, 24, color=\"red\", alpha=0.3, label=\"Ped window\")\n", + "plt.fill_betweenx([0, 1600], 5, 9, color=\"green\", alpha=0.3, label=\"Signal window\")\n", "plt.legend()" ] }, @@ -321,9 +334,9 @@ "metadata": {}, "outputs": [], "source": [ - "data = teldata.waveform[0]\n", + "data = tel_event.r0.waveform[0]\n", "peds = data[:, 19:24].mean(axis=1) # mean of samples 20 to 29 for all pixels\n", - "sums = data[:, 5:9].sum(axis=1)/(13-8) # simple sum integration" + "sums = data[:, 5:9].sum(axis=1) / (13 - 8) # simple sum integration" ] }, { @@ -332,7 +345,7 @@ "metadata": {}, "outputs": [], "source": [ - "phist = plt.hist(peds, bins=50, range=[0,150])\n", + "phist = plt.hist(peds, bins=50)\n", "plt.title(\"Pedestal Distribution of all pixels for a single event\")" ] }, @@ -368,8 +381,8 @@ "outputs": [], "source": [ "# we can also subtract the pedestals from the traces themselves, which would be needed to compare peaks properly\n", - "for ii in range(270,280):\n", - " plt.plot(data[ii] - peds[ii], drawstyle='steps', label=\"pix{}\".format(ii))\n", + "for ii in range(270, 280):\n", + " plt.plot(data[ii] - peds[ii], drawstyle=\"steps\", label=\"pix{}\".format(ii))\n", "plt.legend()" ] }, @@ -399,9 +412,9 @@ "metadata": {}, "outputs": [], "source": [ - "title=\"CT24, run {} event {} ped-sub\".format(event.index.obs_id,event.index.event_id)\n", - "disp = CameraDisplay(camgeom,title=title)\n", - "disp.image = sums - peds \n", + "title = \"CT24, run {} event {} ped-sub\".format(event.index.obs_id, event.index.event_id)\n", + "disp = CameraDisplay(camgeom, title=title)\n", + "disp.image = sums - peds\n", "disp.cmap = plt.cm.RdBu_r\n", "disp.add_colorbar()\n", "disp.set_limits_percent(95) # autoscale" @@ -426,12 +439,14 @@ }, "outputs": [], "source": [ - "for tel in event.r0.tel.keys():\n", + "for tel_id, tel_event in event.tel.items():\n", " plt.figure()\n", - " camgeom = source.subarray.tel[tel].camera.geometry\n", - " title=\"CT{}, run {} event {}\".format(tel,event.index.obs_id,event.index.event_id)\n", - " disp = CameraDisplay(camgeom,title=title)\n", - " disp.image = event.r0.tel[tel].waveform[0].sum(axis=1)\n", + " camgeom = source.subarray.tel[tel_id].camera.geometry\n", + "\n", + " title = f\"CT{tel_id}, run {event.index.obs_id} event {event.index.event_id}\"\n", + "\n", + " disp = CameraDisplay(camgeom, title=title)\n", + " disp.image = tel_event.r0.waveform[0].sum(axis=1)\n", " disp.cmap = plt.cm.RdBu_r\n", " disp.add_colorbar()\n", " disp.set_limits_percent(95)" @@ -460,19 +475,26 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ "pix_ids = np.arange(len(data))\n", - "has_signal = sums > 300\n", + "has_signal = sums > 250\n", "\n", - "widths = np.array([8,]) # peak widths to search for (let's fix it at 8 samples, about the width of the peak)\n", - "peaks = [signal.find_peaks_cwt(trace,widths) for trace in data[has_signal] ]\n", + "widths = np.array(\n", + " [\n", + " 8,\n", + " ]\n", + ") # peak widths to search for (let's fix it at 8 samples, about the width of the peak)\n", + "peaks = [signal.find_peaks_cwt(trace, widths) for trace in data[has_signal]]\n", "\n", - "for p,s in zip(pix_ids[has_signal],peaks):\n", - " print(\"pix{} has peaks at sample {}\".format(p,s))\n", - " plt.plot(data[p], drawstyle='steps-mid')\n", - " plt.scatter(np.array(s),data[p,s])" + "plt.figure()\n", + "for p, s in zip(pix_ids[has_signal], peaks):\n", + " print(\"pix{} has peaks at sample {}\".format(p, s))\n", + " plt.plot(data[p], drawstyle=\"steps-mid\")\n", + " plt.scatter(np.array(s), data[p, s])" ] }, { @@ -481,13 +503,6 @@ "source": [ "clearly the signal needs to be filtered first, or an appropriate wavelet used, but the idea is nice" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -506,7 +521,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.13" + "version": "3.9.16" } }, "nbformat": 4,