diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 248d26a8fac..2e235d95e8d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -75,7 +75,7 @@ jobs:
key: ctapipe-test-data
- name: Prepare mamba installation
- if: matrix.install-method == 'mamba'
+ if: matrix.install-method == 'mamba' && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
env:
PYTHON_VERSION: ${{ matrix.python-version }}
run: |
@@ -83,20 +83,20 @@ jobs:
sed -i -e "s/- python=.*/- python=$PYTHON_VERSION/g" environment.yml
- name: mamba setup
- if: matrix.install-method == 'mamba'
+ if: matrix.install-method == 'mamba' && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
uses: mamba-org/setup-micromamba@v1
with:
environment-file: environment.yml
cache-downloads: true
- name: Python setup
- if: matrix.install-method == 'pip'
+ if: matrix.install-method == 'pip' && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
uses: actions/setup-python@v4
with:
python-version: ${{ matrix.python-version }}
check-latest: true
- - if: ${{ matrix.install-method == 'pip' && runner.os == 'macOS' }}
+ - if: ${{ matrix.install-method == 'pip' && runner.os == 'macOS' }} && contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
name: Fix Python PATH on macOS
# See https://github.com/actions/setup-python/issues/132 and
# https://github.com/actions/setup-python/issues/132#issuecomment-779406058
@@ -107,6 +107,7 @@ jobs:
tee -a ~/.bash_profile <<<'export PATH="$pythonLocation/bin:$PATH"'
- name: Install dependencies
+ if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
run: |
python --version
pip install pytest-cov restructuredtext-lint pytest-xdist 'coverage!=6.3.0'
@@ -115,14 +116,17 @@ jobs:
pip freeze
- name: Static codechecks
+ if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
run: |
restructuredtext-lint README.rst
- name: ctapipe-info
+ if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
run: |
ctapipe-info --all
- name: Tests
+ if: contains(github.event.pull_request.labels.*.name, 'documentation-only') == false
run: |
cd $(mktemp -d)
pytest -n auto --dist loadscope \
diff --git a/README.rst b/README.rst
index b477415b5fe..dc241c4939e 100644
--- a/README.rst
+++ b/README.rst
@@ -71,11 +71,11 @@ or via::
pip install ctapipe
-**Note**: to install a specific version of ctapipe take look at the documentation `here `__.
+**Note**: to install a specific version of ctapipe take look at the documentation `here `__.
**Note**: ``mamba`` is a C++ reimplementation of conda and can be found `here `__.
Note this is *pre-alpha* software and is not yet stable enough for end-users (expect large API changes until the first stable 1.0 release).
Developers should follow the development install instructions found in the
-`documentation `__.
+`documentation `__.
diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index cb4282c5d4c..3fbb88a762c 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -772,6 +772,7 @@ class SimulationConfigContainer(Container):
Configuration parameters of the simulation
"""
+ run_number = Field(np.int32(-1), description="Original sim_telarray run number")
corsika_version = Field(nan, description="CORSIKA version * 1000")
simtel_version = Field(nan, description="sim_telarray version * 1000")
energy_range_min = Field(
@@ -920,12 +921,17 @@ class ReconstructedGeometryContainer(Container):
"reconstructed core position uncertainty along tilted frame Y axis",
unit=u.m,
)
- h_max = Field(nan * u.m, "reconstructed height of the shower maximum", unit=u.m)
+ h_max = Field(
+ nan * u.m,
+ "reconstructed vertical height above sea level of the shower maximum",
+ unit=u.m,
+ )
h_max_uncert = Field(nan * u.m, "uncertainty of h_max", unit=u.m)
+
is_valid = Field(
False,
(
- "direction validity flag. True if the shower direction"
+ "Geometry validity flag. True if the shower geometry"
"was properly reconstructed by the algorithm"
),
)
diff --git a/ctapipe/core/telescope_component.py b/ctapipe/core/telescope_component.py
index 2bac22e615e..6b49dae9227 100644
--- a/ctapipe/core/telescope_component.py
+++ b/ctapipe/core/telescope_component.py
@@ -217,8 +217,10 @@ def __getitem__(self, tel: Optional[Union[int, str]]):
try:
return self._value_for_tel_id[tel]
except KeyError:
+ if tel not in self._subarray.tel:
+ raise KeyError(f"No telescope with id {tel} in subarray")
raise KeyError(
- f"TelescopeParameterLookup: no "
+ "TelescopeParameterLookup: no "
f"parameter value was set for telescope with tel_id="
f"{tel}. Please set it explicitly, "
f"or by telescope type or '*'."
diff --git a/ctapipe/core/tests/test_telescope_component.py b/ctapipe/core/tests/test_telescope_component.py
index 7aae1de47f9..e7f389953c2 100644
--- a/ctapipe/core/tests/test_telescope_component.py
+++ b/ctapipe/core/tests/test_telescope_component.py
@@ -342,3 +342,16 @@ class Foo(TelescopeComponent):
f = Foo(mock_subarray, bar=[("type", "*", 1), ("id", 1, None)])
assert f.bar.tel[1] is None
+
+
+def test_telescope_parameter_nonexistent_telescope(mock_subarray):
+ class Foo(TelescopeComponent):
+ bar = IntTelescopeParameter(
+ default_value=None,
+ allow_none=True,
+ ).tag(config=True)
+
+ foo = Foo(subarray=mock_subarray)
+
+ with pytest.raises(KeyError, match="No telescope with id 0"):
+ foo.bar.tel[0]
diff --git a/ctapipe/image/cleaning.py b/ctapipe/image/cleaning.py
index a62876da286..189a5fb8b97 100644
--- a/ctapipe/image/cleaning.py
+++ b/ctapipe/image/cleaning.py
@@ -15,6 +15,7 @@
__all__ = [
"tailcuts_clean",
"dilate",
+ "time_clustering",
"mars_cleaning_1st_pass",
"fact_image_cleaning",
"apply_time_delta_cleaning",
@@ -27,6 +28,7 @@
from abc import abstractmethod
import numpy as np
+from sklearn.cluster import DBSCAN
from ..core import TelescopeComponent
from ..core.traits import (
@@ -111,6 +113,71 @@ def tailcuts_clean(
)
+def time_clustering(
+ geom,
+ image,
+ time,
+ minpts=5,
+ eps=1.0,
+ time_scale_ns=4.0,
+ space_scale_m=0.25,
+ hard_cut_pe=4,
+):
+ """
+ Clean an image by selecting pixels which pass a time clustering algorithm using DBSCAN.
+ Previously used for HESS [timecleaning]_.
+
+ As a neighbor-based image extractor algorithm can lead to biases in the time reconstruction of noise pixels,
+ specially those next to the shower, a cut in the minimum signal image is applied.
+
+ DBSCAN runs with the reconstructed times and pixel positions after scaling. Scaling is needed because eps
+ is not dimension dependent. If scaling is performed properly, eps can be set to 1. DBSCAN returns the
+ cluster IDs of each point. Pixels associated to cluster ID -1 are classified as noise.
+
+ Parameters
+ ----------
+ geom: `ctapipe.instrument.CameraGeometry`
+ Camera geometry information
+ image: array
+ pixel charge information
+ time: array
+ pixel timing information
+ minpts: int
+ Minimum number of points to consider a cluster
+ eps: float
+ Minimum distance in dbscan
+ time_scale_ns: float
+ Time scale in ns
+ space_scale: float
+ Space scale in m
+ hard_cut_pe: float
+ Hard cut in the number of signal pe
+
+ Returns
+ -------
+ A boolean mask of *clean* pixels.
+ """
+ precut_mask = image > hard_cut_pe
+
+ arr = np.zeros(len(image), dtype=float)
+ arr[~precut_mask] = -1
+
+ pix_x = geom.pix_x.value[precut_mask] / space_scale_m
+ pix_y = geom.pix_y.value[precut_mask] / space_scale_m
+
+ X = np.column_stack((time[precut_mask] / time_scale_ns, pix_x, pix_y))
+
+ labels = DBSCAN(eps=eps, min_samples=minpts).fit_predict(X)
+
+ # no_clusters = len(np.unique(labels))-1 # Could be used for gh separation
+
+ y = np.array(arr[(arr == 0)])
+ y[(labels == -1)] = -1
+ arr[arr == 0] = y
+ mask = arr == 0 # we keep these events
+ return mask
+
+
def mars_cleaning_1st_pass(
geom,
image,
@@ -532,6 +599,43 @@ def __call__(
)
+class TimeCleaner(ImageCleaner):
+ """
+ Clean images using the time clustering cleaning method
+ """
+
+ space_scale_m = FloatTelescopeParameter(
+ default_value=0.25, help="Pixel space scaling parameter in m"
+ ).tag(config=True)
+ time_scale_ns = FloatTelescopeParameter(
+ default_value=4.0, help="Time scale parameter in ns"
+ ).tag(config=True)
+ minpts = IntTelescopeParameter(
+ default_value=5, help="minimum number of points to form a cluster"
+ ).tag(config=True)
+ eps = FloatTelescopeParameter(
+ default_value=1.0, help="minimum distance in DBSCAN"
+ ).tag(config=True)
+ hard_cut_pe = FloatTelescopeParameter(
+ default_value=2.5, help="Hard cut in the number of pe"
+ ).tag(config=True)
+
+ def __call__(
+ self, tel_id: int, image: np.ndarray, arrival_times=None
+ ) -> np.ndarray:
+ """Apply HESS image cleaning. see ImageCleaner.__call__()"""
+ return time_clustering(
+ geom=self.subarray.tel[tel_id].camera.geometry,
+ image=image,
+ time=arrival_times,
+ eps=self.eps.tel[tel_id],
+ space_scale_m=self.space_scale_m.tel[tel_id],
+ time_scale_ns=self.time_scale_ns.tel[tel_id],
+ minpts=self.minpts.tel[tel_id],
+ hard_cut_pe=self.hard_cut_pe.tel[tel_id],
+ )
+
+
class MARSImageCleaner(TailcutsImageCleaner):
"""
1st-pass MARS-like Image cleaner (See `ctapipe.image.mars_cleaning_1st_pass`)
diff --git a/ctapipe/image/tests/test_cleaning.py b/ctapipe/image/tests/test_cleaning.py
index 91b8e85e47f..92b03f1416c 100644
--- a/ctapipe/image/tests/test_cleaning.py
+++ b/ctapipe/image/tests/test_cleaning.py
@@ -334,6 +334,73 @@ def test_apply_time_delta_cleaning(prod3_lst):
assert (test_mask == td_mask).all()
+def test_time_cleaning():
+ geom = CameraGeometry.from_name("LSTCam")
+ charge = np.zeros(geom.n_pixels, dtype=np.float64)
+ peak_time = np.zeros(geom.n_pixels, dtype=np.float64)
+
+ core_pixel = 10
+ core_neighbors = geom.neighbors[core_pixel]
+
+ charge[core_pixel], charge[core_neighbors] = 15, 5
+ peak_time[core_pixel], peak_time[core_neighbors] = (0.0, 0.0)
+
+ mask = cleaning.time_clustering(
+ geom,
+ charge,
+ peak_time,
+ eps=1.0,
+ space_scale_m=0.25,
+ time_scale_ns=4.0,
+ hard_cut_pe=10,
+ minpts=5,
+ )
+
+ assert np.count_nonzero(mask) == 0
+
+ mask = cleaning.time_clustering(
+ geom,
+ charge,
+ peak_time,
+ eps=1.0,
+ space_scale_m=0.25,
+ time_scale_ns=4.0,
+ hard_cut_pe=1,
+ minpts=5,
+ )
+
+ assert np.count_nonzero(mask) == 7
+
+ mask = cleaning.time_clustering(
+ geom,
+ charge,
+ peak_time,
+ eps=1.0,
+ space_scale_m=0.25,
+ time_scale_ns=4.0,
+ hard_cut_pe=1,
+ minpts=8,
+ )
+
+ assert np.count_nonzero(mask) == 0
+
+ charge[core_pixel], charge[core_neighbors] = 15, 5
+ peak_time[core_pixel], peak_time[core_neighbors] = (0.0, 30.0)
+
+ mask = cleaning.time_clustering(
+ geom,
+ charge,
+ peak_time,
+ eps=1.0,
+ space_scale_m=0.25,
+ time_scale_ns=4.0,
+ hard_cut_pe=1,
+ minpts=5,
+ )
+
+ assert np.count_nonzero(mask) == 6
+
+
def test_time_constrained_clean():
geom = CameraGeometry.from_name("LSTCam")
charge = np.zeros(geom.n_pixels, dtype=np.float64)
diff --git a/ctapipe/instrument/camera/geometry.py b/ctapipe/instrument/camera/geometry.py
index dd55b14425b..8d22dc6e2b6 100644
--- a/ctapipe/instrument/camera/geometry.py
+++ b/ctapipe/instrument/camera/geometry.py
@@ -800,7 +800,7 @@ def pixel_moment_matrix(self):
x = self.pix_x.value
y = self.pix_y.value
- return np.row_stack(
+ return np.vstack(
[
x,
y,
@@ -954,12 +954,13 @@ def position_to_pix_index(self, x, y):
pix_indices: Pixel index or array of pixel indices. Returns -1 if position falls
outside camera
"""
-
if not self._all_pixel_areas_equal:
logger.warning(
" Method not implemented for cameras with varying pixel sizes"
)
unit = x.unit
+ scalar = x.ndim == 0
+
points_searched = np.dstack([x.to_value(unit), y.to_value(unit)])
circum_rad = self._pixel_circumradius[0].to_value(unit)
kdtree = self._kdtree
@@ -969,8 +970,9 @@ def position_to_pix_index(self, x, y):
del dist
pix_indices = pix_indices.flatten()
+ invalid = np.iinfo(pix_indices.dtype).min
# 1. Mark all points outside pixel circumeference as lying outside camera
- pix_indices[pix_indices == self.n_pixels] = -1
+ pix_indices[pix_indices == self.n_pixels] = invalid
# 2. Accurate check for the remaing cases (within circumference, but still outside
# camera). It is first checked if any border pixel numbers are returned.
@@ -1006,17 +1008,9 @@ def position_to_pix_index(self, x, y):
)
del dist_check
if index_check != insidepix_index:
- pix_indices[index] = -1
-
- # print warning:
- for index in np.where(pix_indices == -1)[0]:
- logger.warning(
- " Coordinate ({} m, {} m) lies outside camera".format(
- points_searched[0][index, 0], points_searched[0][index, 1]
- )
- )
+ pix_indices[index] = invalid
- return pix_indices if len(pix_indices) > 1 else pix_indices[0]
+ return np.squeeze(pix_indices) if scalar else pix_indices
@staticmethod
def simtel_shape_to_type(pixel_shape):
diff --git a/ctapipe/instrument/camera/image_conversion.py b/ctapipe/instrument/camera/image_conversion.py
index c4131a77cc5..759375449ca 100644
--- a/ctapipe/instrument/camera/image_conversion.py
+++ b/ctapipe/instrument/camera/image_conversion.py
@@ -1,6 +1,5 @@
-import numpy as np
import astropy.units as u
-
+import numpy as np
__all__ = [
"unskew_hex_pixel_grid",
@@ -294,8 +293,8 @@ def get_orthogonal_grid_edges(pix_x, pix_y, scale_aspect=True):
# with the maximal extension of the axes and the size of the pixels,
# determine the number of bins in each direction
- n_bins_x = int(np.round_(np.abs(np.max(pix_x) - np.min(pix_x)) / d_x) + 2)
- n_bins_y = int(np.round_(np.abs(np.max(pix_y) - np.min(pix_y)) / d_y) + 2)
+ n_bins_x = int(np.round(np.abs(np.max(pix_x) - np.min(pix_x)) / d_x) + 2)
+ n_bins_y = int(np.round(np.abs(np.max(pix_y) - np.min(pix_y)) / d_y) + 2)
x_edges = np.linspace(pix_x.min(), pix_x.max(), n_bins_x)
y_edges = np.linspace(pix_y.min(), pix_y.max(), n_bins_y)
diff --git a/ctapipe/instrument/camera/tests/test_geometry.py b/ctapipe/instrument/camera/tests/test_geometry.py
index 10410e0efc6..8e3a7545d95 100644
--- a/ctapipe/instrument/camera/tests/test_geometry.py
+++ b/ctapipe/instrument/camera/tests/test_geometry.py
@@ -68,10 +68,20 @@ def test_load_lst_camera(prod5_lst):
def test_position_to_pix_index(prod5_lst):
"""test that we can lookup a pixel from a coordinate"""
+ geometry = prod5_lst.camera.geometry
+
x, y = (0.80 * u.m, 0.79 * u.m)
- pix_id = prod5_lst.camera.geometry.position_to_pix_index(x, y)
+
+ pix_id = geometry.position_to_pix_index(x, y)
+
assert pix_id == 1575
+ pix_ids = geometry.position_to_pix_index([0.8, 0.8] * u.m, [0.79, 0.79] * u.m)
+ np.testing.assert_array_equal(pix_ids, [1575, 1575])
+
+ assert len(geometry.position_to_pix_index([] * u.m, [] * u.m)) == 0
+ assert geometry.position_to_pix_index(5 * u.m, 5 * u.m) == np.iinfo(int).min
+
def test_find_neighbor_pixels():
"""test basic neighbor functionality"""
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 632277a566c..70e8e93f81a 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -49,7 +49,7 @@
from ..coordinates import CameraFrame, shower_impact_distance
from ..core import Map
from ..core.provenance import Provenance
-from ..core.traits import Bool, ComponentName, Float, Undefined, UseEnum
+from ..core.traits import Bool, ComponentName, Float, Integer, Undefined, UseEnum
from ..instrument import (
CameraDescription,
CameraGeometry,
@@ -483,6 +483,12 @@ class SimTelEventSource(EventSource):
),
).tag(config=True)
+ override_obs_id = Integer(
+ default_value=None,
+ allow_none=True,
+ help="Use the given obs_id instead of the run number from sim_telarray",
+ ).tag(config=True)
+
def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs):
"""
EventSource for simtelarray files using the pyeventio library.
@@ -701,11 +707,10 @@ def _make_dummy_location(self):
SimulationConfigContainer), but with the lat/lon set to (0,0), i.e.
on "Null Island"
"""
- obs_id = self.file_.header["run"]
return EarthLocation(
lon=0 * u.deg,
lat=0 * u.deg,
- height=self._simulation_config[obs_id].prod_site_alt,
+ height=self._simulation_config[self.obs_id].prod_site_alt,
)
@staticmethod
@@ -743,7 +748,7 @@ def _generate_events(self):
pseudo_event_id -= 1
event_id = pseudo_event_id
- obs_id = self.file_.header["run"]
+ obs_id = self.obs_id
trigger = self._fill_trigger_info(array_event)
if trigger.event_type == EventType.SUBARRAY:
@@ -1007,12 +1012,14 @@ def _parse_simulation_header(self):
point in time, this dictionary will always have
length 1.
"""
- assert len(self.obs_ids) == 1
- obs_id = self.obs_ids[0]
+ obs_id = self.obs_id
# With only one run, we can take the first entry:
mc_run_head = self.file_.mc_run_headers[-1]
simulation_config = SimulationConfigContainer(
+ run_number=self.file_.header["run"],
+ detector_prog_start=mc_run_head["detector_prog_start"],
+ detector_prog_id=mc_run_head["detector_prog_id"],
corsika_version=mc_run_head["shower_prog_vers"],
simtel_version=mc_run_head["detector_prog_vers"],
energy_range_min=mc_run_head["E_range"][0] * u.TeV,
@@ -1024,8 +1031,6 @@ def _parse_simulation_header(self):
spectral_index=mc_run_head["spectral_index"],
shower_prog_start=mc_run_head["shower_prog_start"],
shower_prog_id=mc_run_head["shower_prog_id"],
- detector_prog_start=mc_run_head["detector_prog_start"],
- detector_prog_id=mc_run_head["detector_prog_id"],
n_showers=mc_run_head["n_showers"],
shower_reuse=mc_run_head["n_use"],
max_alt=_clip_altitude_if_close(mc_run_head["alt_range"][1]) * u.rad,
@@ -1056,11 +1061,20 @@ def _fill_scheduling_and_observation_blocks(self):
"""
az, alt = self.file_.header["direction"]
- obs_id = self.file_.header["run"]
+
+ # this event source always contains only a single OB, so we can
+ # also assign a single obs_id
+ if self.override_obs_id is not None:
+ self.obs_id = self.override_obs_id
+ else:
+ self.obs_id = self.file_.header["run"]
+
+ # simulations at the moment do not have SBs, use OB id
+ self.sb_id = self.obs_id
sb_dict = {
- obs_id: SchedulingBlockContainer(
- sb_id=np.uint64(obs_id), # simulations have no SBs, use OB id
+ self.sb_id: SchedulingBlockContainer(
+ sb_id=np.uint64(self.sb_id),
sb_type=SchedulingBlockType.OBSERVATION,
producer_id="simulation",
observing_mode=ObservingMode.UNKNOWN,
@@ -1069,9 +1083,9 @@ def _fill_scheduling_and_observation_blocks(self):
}
ob_dict = {
- obs_id: ObservationBlockContainer(
- obs_id=np.uint64(obs_id),
- sb_id=np.uint64(obs_id), # see comment above
+ self.obs_id: ObservationBlockContainer(
+ obs_id=np.uint64(self.obs_id),
+ sb_id=np.uint64(self.sb_id),
producer_id="simulation",
state=ObservationBlockState.COMPLETED_SUCCEDED,
subarray_pointing_lat=alt * u.rad,
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index b4d3c0054ff..116a8a91f8a 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -604,3 +604,25 @@ def test_starting_grammage():
with SimTelEventSource(path, focal_length_choice="EQUIVALENT") as source:
e = next(iter(source))
assert e.simulation.shower.starting_grammage == 580 * u.g / u.cm**2
+
+
+@pytest.mark.parametrize("override_obs_id,expected_obs_id", [(None, 1), (5, 5)])
+def test_override_obs_id(override_obs_id, expected_obs_id, prod5_gamma_simtel_path):
+ """Test for the override_obs_id option"""
+ original_run_number = 1
+
+ with SimTelEventSource(
+ prod5_gamma_simtel_path, override_obs_id=override_obs_id
+ ) as s:
+ assert s.obs_id == expected_obs_id
+ assert s.obs_ids == [expected_obs_id]
+
+ assert s.simulation_config.keys() == {expected_obs_id}
+ assert s.observation_blocks.keys() == {expected_obs_id}
+ assert s.scheduling_blocks.keys() == {expected_obs_id}
+
+ # this should alway be the original run number
+ assert s.simulation_config[s.obs_id].run_number == original_run_number
+
+ for e in s:
+ assert e.index.obs_id == expected_obs_id
diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index d7288b58576..0906cf0ff4b 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -44,6 +44,7 @@
)
FOV_ANGULAR_DISTANCE_LIMIT_RAD = (45 * u.deg).to_value(u.rad)
+H_MAX_UPPER_LIMIT_M = 100_000
def _far_outside_fov(fov_lat, fov_lon):
@@ -260,7 +261,8 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
sky_pos = nom.transform_to(array_pointing.frame)
tilt = SkyCoord(x=core_x * u.m, y=core_y * u.m, z=0 * u.m, frame=tilted_frame)
grd = project_to_ground(tilt)
- x_max = self.reconstruct_xmax(
+
+ h_max = self.reconstruct_h_max(
nom.fov_lon,
nom.fov_lat,
tilt.x,
@@ -287,8 +289,8 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
is_valid=True,
alt_uncert=src_error.to(u.rad),
az_uncert=src_error.to(u.rad),
- h_max=x_max,
- h_max_uncert=u.Quantity(np.nan * x_max.unit),
+ h_max=h_max,
+ h_max_uncert=u.Quantity(np.nan * h_max.unit),
goodness_of_fit=np.nan,
prefix=self.__class__.__name__,
)
@@ -309,7 +311,7 @@ def reconstruct_nominal(self, hillas_parameters):
"""
if len(hillas_parameters) < 2:
- return None # Throw away events with < 2 images
+ return (np.nan, np.nan, np.nan, np.nan) # Throw away events with < 2 images
# Find all pairs of Hillas parameters
combos = itertools.combinations(list(hillas_parameters.values()), 2)
@@ -381,7 +383,8 @@ def reconstruct_tilted(self, hillas_parameters, tel_x, tel_y):
core uncertainty X
"""
if len(hillas_parameters) < 2:
- return None # Throw away events with < 2 images
+ return (np.nan, np.nan, np.nan, np.nan) # Throw away events with < 2 images
+
hill_list = list()
tx = list()
ty = list()
@@ -437,7 +440,7 @@ def reconstruct_tilted(self, hillas_parameters, tel_x, tel_y):
return x_pos, y_pos, np.sqrt(var_x), np.sqrt(var_y)
- def reconstruct_xmax(
+ def reconstruct_h_max(
self, source_x, source_y, core_x, core_y, hillas_parameters, tel_x, tel_y, zen
):
"""
@@ -495,25 +498,19 @@ def reconstruct_xmax(
np.array(ty),
)
weight = np.array(amp)
- mean_height = np.sum(height * weight) / np.sum(weight)
+ mean_distance = np.average(height, weights=weight)
# This value is height above telescope in the tilted system,
# we should convert to height above ground
- mean_height *= np.cos(zen)
+ mean_height = mean_distance * np.cos(zen.to_value(u.rad))
# Add on the height of the detector above sea level
- mean_height += 2100 # TODO: replace with instrument info
-
- if mean_height > 100000 or np.isnan(mean_height):
- mean_height = 100000
+ mean_height += self.subarray.reference_location.geodetic.height.to_value(u.m)
- mean_height *= u.m
- # Lookup this height in the depth tables, the convert Hmax to Xmax
- # x_max = self.thickness_profile(mean_height.to(u.km))
- # Convert to slant depth
- # x_max /= np.cos(zen)
+ if mean_height > H_MAX_UPPER_LIMIT_M:
+ mean_height = np.nan
- return mean_height
+ return u.Quantity(mean_height, u.m)
@staticmethod
def intersect_lines(xp1, yp1, phi1, xp2, yp2, phi2):
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index 1fbd305d7c2..0f82e2bc51d 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -19,6 +19,7 @@
altaz_to_righthanded_cartesian,
project_to_ground,
)
+from ..instrument import SubarrayDescription
from .reconstructor import (
HillasGeometryReconstructor,
InvalidWidthException,
@@ -106,7 +107,7 @@ class that reconstructs the direction of an atmospheric shower
"""
- def __init__(self, subarray, **kwargs):
+ def __init__(self, subarray: SubarrayDescription, **kwargs):
super().__init__(subarray=subarray, **kwargs)
_cam_radius_m = {
cam: cam.geometry.guess_radius().to_value(u.m)
@@ -181,7 +182,10 @@ def __call__(self, event):
_, lat, lon = cartesian_to_spherical(*direction)
# estimate max height of shower
- h_max = self.estimate_h_max(cog_cartesian, telescope_positions)
+ h_max = (
+ self.estimate_relative_h_max(cog_cartesian, telescope_positions)
+ + self.subarray.reference_location.geodetic.height
+ )
# az is clockwise, lon counter-clockwise, make sure it stays in [0, 2pi)
az = Longitude(-lon)
@@ -412,18 +416,22 @@ def estimate_core_position(array_pointing, psi, positions):
return core_pos_ground, core_pos_tilted
@staticmethod
- def estimate_h_max(cog_vectors, positions):
- """
- Estimate the max height by intersecting the lines of the cog directions of each telescope.
+ def estimate_relative_h_max(cog_vectors, positions):
+ """Estimate the relative (to the observatory) vertical height of
+ shower-max by intersecting the lines of the cog directions of each
+ telescope.
Returns
-------
- astropy.unit.Quantity
- the estimated max height
+ astropy.unit.Quantity:
+ the estimated height above observatory level (not sea level) of the
+ shower-max point
+
"""
positions = positions.cartesian.xyz.T.to_value(u.m)
- # not sure if its better to return the length of the vector of the z component
- return u.Quantity(
- np.linalg.norm(line_line_intersection_3d(cog_vectors, positions)),
+ shower_max = u.Quantity(
+ line_line_intersection_3d(cog_vectors, positions),
u.m,
)
+
+ return shower_max[2] # the z-coordinate only
diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py
index 976ecb6079e..1ae0bc78158 100644
--- a/ctapipe/reco/sklearn.py
+++ b/ctapipe/reco/sklearn.py
@@ -416,7 +416,7 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table
result,
ReconstructedEnergyContainer,
prefix=self.prefix,
- stereo=False,
+ add_tel_prefix=True,
)
return {ReconstructionProperty.ENERGY: result}
@@ -480,7 +480,10 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table
}
)
add_defaults_and_meta(
- result, ParticleClassificationContainer, prefix=self.prefix, stereo=False
+ result,
+ ParticleClassificationContainer,
+ prefix=self.prefix,
+ add_tel_prefix=True,
)
return {ReconstructionProperty.PARTICLE_TYPE: result}
@@ -770,7 +773,8 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table
disp_result,
DispContainer,
prefix=f"{self.prefix}_parameter",
- stereo=False,
+ # disp is always per telescope, so no need to add the prefix
+ add_tel_prefix=False,
)
psi = table["hillas_psi"].quantity.to_value(u.rad)
@@ -797,7 +801,7 @@ def predict_table(self, key, table: Table) -> Dict[ReconstructionProperty, Table
altaz_result,
ReconstructedGeometryContainer,
prefix=self.prefix,
- stereo=False,
+ add_tel_prefix=True,
)
return {
diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py
index 9fe7c57d4de..fd8ca438495 100644
--- a/ctapipe/reco/tests/test_HillasReconstructor.py
+++ b/ctapipe/reco/tests/test_HillasReconstructor.py
@@ -44,12 +44,9 @@ def test_h_max_results():
)
cog_cart = altaz_to_righthanded_cartesian(cog_alt, cog_az)
-
- # creating the fit class and setting the the great circle member
-
- # performing the direction fit with the minimisation algorithm
- # and a seed that is perpendicular to the up direction
- h_max_reco = HillasReconstructor.estimate_h_max(cog_cart, positions)
+ h_max_reco = HillasReconstructor.estimate_relative_h_max(
+ cog_vectors=cog_cart, positions=positions
+ )
# the results should be close to the direction straight up
assert u.isclose(h_max_reco, 1 * u.km)
@@ -210,7 +207,6 @@ def test_reconstruction_against_simulation_camera_frame(
],
)
def test_CameraFrame_against_TelescopeFrame(filename):
-
input_file = get_dataset_path(filename)
# "gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz"
# )
@@ -243,7 +239,6 @@ def test_CameraFrame_against_TelescopeFrame(filename):
reconstructed_events = 0
for event_telescope_frame in source:
-
calib(event_telescope_frame)
# make a copy of the calibrated event for the camera frame case
# later we clean and paramretrize the 2 events in the same way
@@ -281,7 +276,6 @@ def test_CameraFrame_against_TelescopeFrame(filename):
elif isinstance(cam, list):
assert cam == tel
else:
-
if cam == 0 or tel == 0:
kwargs["atol"] = 1e-6
assert np.isclose(cam, tel, **kwargs)
diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py
index 94492aba3d6..924699b74ed 100644
--- a/ctapipe/reco/tests/test_hillas_intersection.py
+++ b/ctapipe/reco/tests/test_hillas_intersection.py
@@ -30,8 +30,8 @@ def test_intersect():
def test_parallel():
"""
- Simple test to check the intersection of lines. Try to intersect positions at (0,0) and (0,1)
- with angles parallel and check the behaviour
+ ? Simple test to check the intersection of lines. Try to intersect positions at (0,0) and (0,1)
+ with angles parallel and check the behaviour
"""
x1 = 0
y1 = 0
@@ -83,7 +83,7 @@ def test_intersection_xmax_reco(example_subarray):
),
}
- x_max = hill_inter.reconstruct_xmax(
+ h_max = hill_inter.reconstruct_h_max(
source_x=nom_pos_reco.fov_lon,
source_y=nom_pos_reco.fov_lat,
core_x=0 * u.m,
@@ -93,7 +93,8 @@ def test_intersection_xmax_reco(example_subarray):
tel_y={1: (0 * u.m), 2: (150 * u.m)},
zen=zen_pointing,
)
- print(x_max)
+
+ assert h_max > 0 * u.m
def test_intersection_reco_impact_point_tilted(example_subarray):
diff --git a/ctapipe/reco/utils.py b/ctapipe/reco/utils.py
index bf26d948074..48e14445b17 100644
--- a/ctapipe/reco/utils.py
+++ b/ctapipe/reco/utils.py
@@ -1,4 +1,4 @@
-def add_defaults_and_meta(table, container, prefix=None, stereo=True):
+def add_defaults_and_meta(table, container, prefix=None, add_tel_prefix=False):
"""
Fill column descriptions and default values into table for container
@@ -10,21 +10,21 @@ def add_defaults_and_meta(table, container, prefix=None, stereo=True):
the container class to add columns and descriptions to the table
prefix : str
prefix for the column names
- stereo : bool
- If False, add a ``tel_`` prefix to the column names to signal it's
+ add_tel_prefix : bool
+ If True, add a ``tel_`` prefix to the column names to signal it's
telescope-wise quantity
"""
if prefix is None:
prefix = container.default_prefix
for name, field in container.fields.items():
- if not stereo and name == "telescopes":
+ if add_tel_prefix and name == "telescopes":
continue
- if stereo:
- colname = f"{prefix}_{name}"
- else:
+ if add_tel_prefix:
colname = f"{prefix}_tel_{name}"
+ else:
+ colname = f"{prefix}_{name}"
if colname not in table.colnames and field.default is not None:
table[colname] = field.default
diff --git a/ctapipe/tools/apply_models.py b/ctapipe/tools/apply_models.py
index 0a6285dcbab..2df691d8bf4 100644
--- a/ctapipe/tools/apply_models.py
+++ b/ctapipe/tools/apply_models.py
@@ -177,7 +177,13 @@ def _apply(self, reconstructor, tel_tables, start, stop):
for tel_id, table in tel_tables.items():
tel = self.loader.subarray.tel[tel_id]
- if tel not in reconstructor._models:
+ if len(table) == 0:
+ self.log.info("No events for telescope %d", tel_id)
+ continue
+
+ try:
+ predictions = reconstructor.predict_table(tel, table)
+ except KeyError:
self.log.warning(
"No model in %s for telescope type %s, skipping tel %d",
reconstructor,
@@ -186,11 +192,6 @@ def _apply(self, reconstructor, tel_tables, start, stop):
)
continue
- if len(table) == 0:
- self.log.info("No events for telescope %d", tel_id)
- continue
-
- predictions = reconstructor.predict_table(tel, table)
for prop, prediction_table in predictions.items():
# copy/overwrite columns into full feature table
new_columns = prediction_table.colnames
diff --git a/ctapipe/tools/quickstart.py b/ctapipe/tools/quickstart.py
index 317d5723ab2..319ccf316c1 100644
--- a/ctapipe/tools/quickstart.py
+++ b/ctapipe/tools/quickstart.py
@@ -17,6 +17,7 @@
"ml_preprocessing_config.yaml",
"train_energy_regressor.yaml",
"train_particle_classifier.yaml",
+ "train_disp_reconstructor.yaml",
]
README_TEXT = f"""
@@ -55,13 +56,16 @@
ctapipe-process --help-all
```
-## ctapipe-train-energy-regressor / ctapipe-train-particle-classifier configs
+## ctapipe-train-energy-regressor / ctapipe-train-particle-classifier / ctapipe-train-disp-reconstructor configs
Included here are also base configurations for training machine learning (ML)
-models for energy regression and gamma/hadron separation.
+models for energy regression, gamma/hadron separation and disp origin reconstruction.
+NOTE: As these files are used for unit tests, they are optimized for very fast training
+and will not result in well performing models.
- `train_energy_regressor.yaml`: configuration of energy regression model
- `train_particle_classifier.yaml`: configuration of particle classification model
+- `train_disp_reconstructor.yaml`: configuration of disp reconstruction models
This file was generated using ctapipe version {VERSION}
diff --git a/ctapipe/tools/tests/test_apply_models.py b/ctapipe/tools/tests/test_apply_models.py
index 13f8edb805a..32b0e4bcc86 100644
--- a/ctapipe/tools/tests/test_apply_models.py
+++ b/ctapipe/tools/tests/test_apply_models.py
@@ -103,7 +103,6 @@ def test_apply_all(
energy_regressor_path,
particle_classifier_path,
disp_reconstructor_path,
- dl2_shower_geometry_file_lapalma,
tmp_path,
):
from ctapipe.tools.apply_models import ApplyModels
@@ -211,6 +210,7 @@ def test_apply_all(
assert f"{prefix_disp}_tel_is_valid" in tel_events.colnames
assert f"{prefix_disp}_parameter_norm" in tel_events.colnames
assert f"{prefix_disp}_parameter_is_valid" in tel_events.colnames
+ assert f"{prefix_disp}_parameter_tel_is_valid" not in tel_events.colnames
# check that the "--no-dl1-parameters" option worked
assert "hillas_intensity" not in tel_events.colnames
diff --git a/ctapipe/tools/train_disp_reconstructor.py b/ctapipe/tools/train_disp_reconstructor.py
index 5953bda8038..80352a0bb29 100644
--- a/ctapipe/tools/train_disp_reconstructor.py
+++ b/ctapipe/tools/train_disp_reconstructor.py
@@ -1,6 +1,10 @@
+"""
+Tool for training the DispReconstructor
+"""
import astropy.units as u
import numpy as np
+from ctapipe.containers import CoordinateFrameType
from ctapipe.core import Tool
from ctapipe.core.traits import Bool, Int, IntTelescopeParameter, Path
from ctapipe.exceptions import TooFewEvents
@@ -8,6 +12,10 @@
from ctapipe.reco import CrossValidator, DispReconstructor
from ctapipe.reco.preprocessing import check_valid_rows, horizontal_to_telescope
+__all__ = [
+ "TrainDispReconstructor",
+]
+
class TrainDispReconstructor(Tool):
"""
@@ -15,7 +23,8 @@ class TrainDispReconstructor(Tool):
The tool first performs a cross validation to give an initial estimate
on the quality of the estimation and then finally trains two models
- (|disp| and sign(disp)) per telescope type on the full dataset.
+ (estimating ``norm(disp)`` and ``sign(disp)`` respectively) per
+ telescope type on the full dataset.
"""
name = "ctapipe-train-disp-reconstructor"
@@ -54,9 +63,9 @@ class TrainDispReconstructor(Tool):
project_disp = Bool(
default_value=False,
help=(
- "If true, true |disp| is the distance between shower cog and"
+ "If true, ``true_disp`` is the distance between shower cog and"
" the true source position along the reconstructed main shower axis."
- "If false, true |disp| is the distance between shower cog"
+ "If false, ``true_disp`` is the distance between shower cog"
" and the true source position."
),
).tag(config=True)
@@ -127,6 +136,13 @@ def _read_table(self, telescope_type):
f"No events after quality query for telescope type {telescope_type}"
)
+ if not np.all(
+ table["subarray_pointing_frame"] == CoordinateFrameType.ALTAZ.value
+ ):
+ raise ValueError(
+ "Pointing information for training data has to be provided in horizontal coordinates"
+ )
+
table = self.models.feature_generator(table, subarray=self.loader.subarray)
table[self.models.target] = self._get_true_disp(table)
diff --git a/docs/api-reference/tools/index.rst b/docs/api-reference/tools/index.rst
index 0d1c2dc8e6d..e9e9823e253 100644
--- a/docs/api-reference/tools/index.rst
+++ b/docs/api-reference/tools/index.rst
@@ -93,5 +93,8 @@ Reference/API
.. automodapi:: ctapipe.tools.train_particle_classifier
:no-inheritance-diagram:
+.. automodapi:: ctapipe.tools.train_disp_reconstructor
+ :no-inheritance-diagram:
+
.. automodapi:: ctapipe.tools.apply_models
:no-inheritance-diagram:
diff --git a/docs/bibliography.rst b/docs/bibliography.rst
index 3ae82af0d7e..2b37cc06369 100644
--- a/docs/bibliography.rst
+++ b/docs/bibliography.rst
@@ -28,6 +28,8 @@ References
methodologies for the Cherenkov Telescope Array",
Ph.D. thesis, Academic Year 2018/2019.
+.. [timecleaning] Steinmaßl, S. F., "Probing particle acceleration in stellar binary systems using gamma-ray observations", Ph.D. thesis, 2023.
+
.. [lts_regression] P. J. Rousseeuw, K. van Driessen, "Computing LTS Regression for Large Data Sets",
Data Mining and Knowledge Discovery, 12, 29-45, DOI: 10.1007/s10618-005-0024-4
diff --git a/docs/changes/2397.api.rst b/docs/changes/2397.api.rst
new file mode 100644
index 00000000000..45b64781bb8
--- /dev/null
+++ b/docs/changes/2397.api.rst
@@ -0,0 +1,8 @@
+``CameraGeometry.position_to_pix_index`` will now return the minimum integer value for invalid
+pixel coordinates instead of -1 due to the danger of using -1 as an index in python accessing
+the last element of a data array for invalid pixels.
+The function will now also no longer raise an error if the arguments are empty arrays and instead
+just return an empty index array.
+The function will also no longer log a warning in case of coordinates that do not match a camera pixel.
+The function is very low-level and if not finding a pixel at the tested position warrants a warning or
+is expected will depend on the calling code.
diff --git a/docs/changes/2401.feature.rst b/docs/changes/2401.feature.rst
new file mode 100644
index 00000000000..24cdf7bb76b
--- /dev/null
+++ b/docs/changes/2401.feature.rst
@@ -0,0 +1,2 @@
+Add an image cleaning algorithm that clusters points in a high-density region in a space with one temporal dimension and
+two spatial dimensions
\ No newline at end of file
diff --git a/docs/changes/2403.bugfix.rst b/docs/changes/2403.bugfix.rst
new file mode 100644
index 00000000000..042b14c59bd
--- /dev/null
+++ b/docs/changes/2403.bugfix.rst
@@ -0,0 +1,14 @@
+Fixed the definition of ``h_max``, which was both inconsistent between
+`~ctapipe.reco.HillasReconstructor` and `~ctapipe.reco.HillasIntersection`
+implementations, and was also incorrect since it was measured from the
+observatory elevation rather than from sea level.
+
+The value of ``h_max`` is now defined as the height above sea level of the
+shower-max point (in meters), not the distance to that point. Therefore it is
+not corrected for the zenith angle of the shower. This is consistent with the
+options currently used for *CORSIKA*, where the *SLANT* option is set to false,
+meaning heights are actual heights not distances from the impact point, and
+``x_max`` is a *depth*, not a *slant depth*. Note that this definition may be
+inconsistent with other observatories where slant-depths are used, and also note
+that the slant depth or distance to shower max are the more useful quantities
+for shower physics.
diff --git a/docs/changes/2406.maintenance.rst b/docs/changes/2406.maintenance.rst
new file mode 100644
index 00000000000..2d4cf50ceac
--- /dev/null
+++ b/docs/changes/2406.maintenance.rst
@@ -0,0 +1 @@
+Updated some numpy calls to not use depreciated functions.
diff --git a/docs/changes/2411.features.rst b/docs/changes/2411.features.rst
new file mode 100644
index 00000000000..d81a88b977f
--- /dev/null
+++ b/docs/changes/2411.features.rst
@@ -0,0 +1,3 @@
+Add option ``override_obs_id`` to ``SimTelEventSource`` which allows
+assigning new, unique ``obs_ids`` in case productions re-use CORSIKA run
+numbers.
diff --git a/docs/changes/2414.bugfix.rst b/docs/changes/2414.bugfix.rst
new file mode 100644
index 00000000000..0e6c4243c33
--- /dev/null
+++ b/docs/changes/2414.bugfix.rst
@@ -0,0 +1,2 @@
+Add the example config for ctapipe-train-disp-reconstructor
+to the list of configs generated by ctapipe-quickstart.
diff --git a/docs/changes/2418.bugfix.rst b/docs/changes/2418.bugfix.rst
new file mode 100644
index 00000000000..5af5d6d3661
--- /dev/null
+++ b/docs/changes/2418.bugfix.rst
@@ -0,0 +1 @@
+Do not use a hidden attribute of ``SKLearnReconstructor`` in ``ctapipe-apply-models``.
diff --git a/docs/changes/2420.bugfix.rst b/docs/changes/2420.bugfix.rst
new file mode 100644
index 00000000000..aa21eafa4fb
--- /dev/null
+++ b/docs/changes/2420.bugfix.rst
@@ -0,0 +1 @@
+Add docstring for ``ctapipe-train-disp-reconstructor``.
diff --git a/docs/changes/2429.feature.rst b/docs/changes/2429.feature.rst
new file mode 100644
index 00000000000..6dc4cec824d
--- /dev/null
+++ b/docs/changes/2429.feature.rst
@@ -0,0 +1,3 @@
+In case no configuration is found for a telescope in ``TelescopeParameter``,
+it is now checked whether the telescope exists at all to provide a better
+error message.
diff --git a/docs/changes/2431.bugfix.rst b/docs/changes/2431.bugfix.rst
new file mode 100644
index 00000000000..d5da9c18960
--- /dev/null
+++ b/docs/changes/2431.bugfix.rst
@@ -0,0 +1,2 @@
+Check that the array pointing is given in horizontal coordinates
+before training a ``DispReconstructor``.
diff --git a/docs/changes/2440.bugfix.rst b/docs/changes/2440.bugfix.rst
new file mode 100644
index 00000000000..435d4288113
--- /dev/null
+++ b/docs/changes/2440.bugfix.rst
@@ -0,0 +1 @@
+Fix additional, unwanted columns being written into disp prediction output.
diff --git a/docs/conf.py b/docs/conf.py
index ba59f8e41f0..4bcb2131e40 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -131,8 +131,14 @@ def setup(app):
("py:class", "traitlets.config.application.Application"),
("py:class", "traitlets.utils.sentinel.Sentinel"),
("py:class", "traitlets.traitlets.ObserveHandler"),
+ ("py:class", "dict[K, V]"),
+ ("py:class", "G"),
+ ("py:class", "t.Sequence"),
+ ("py:class", "StrDict"),
+ ("py:class", "ClassesType"),
("py:obj", "traitlets.traitlets.G"),
("py:obj", "traitlets.traitlets.S"),
+ ("py:obj", "traitlets.traitlets.T"),
("py:class", "traitlets.traitlets.T"),
("py:class", "re.Pattern[t.Any]"),
("py:class", "Sentinel"),
@@ -288,6 +294,7 @@ def setup(app):
"version_match": version_match,
"json_url": json_url,
},
+ "navigation_with_keys": False,
"use_edit_page_button": True,
"icon_links_label": "Quick Links",
"icon_links": [
diff --git a/docs/developer-guide/ceps/proposed/cep-003-remove-image-parameters-in-camera-frame.rst b/docs/developer-guide/ceps/proposed/cep-003-remove-image-parameters-in-camera-frame.rst
new file mode 100644
index 00000000000..463df4ac2a5
--- /dev/null
+++ b/docs/developer-guide/ceps/proposed/cep-003-remove-image-parameters-in-camera-frame.rst
@@ -0,0 +1,52 @@
+.. _cep-003:
+
+
+************************************************************
+CEP 3 - Dropping support for image parameters in CameraFrame
+************************************************************
+
+* Status: draft
+* Discussion: NA
+* Date accepted: NA
+* Last revised: 2023-09-22
+* Author: Maximilian Linhoff
+* Created: 2023-09-22
+
+Abstract
+========
+
+Currently, ctapipe supports computing all image parameters in two variants:
+
+* Using a ``CameraGeometry`` where pixel coordinates are expressed in ``CameraFrame``, i.e.
+ in length units (most commonly meters) on the camera focal plane.
+* Using a ``CameraGeometry`` where pixel coordinates are expressed in ``TelescopeFrame``, i.e.
+ in angular units (most commonly degree) on sky.
+
+We propose to drop support for the first, to simplify code in multiple places and reduce
+possibility for confusing the two similar variants of the image parameters.
+
+The overhead of supporting both ``TelescopeFrame`` and ``CameraFrame`` representations
+of the image parameters is quite significant, as it e.g. requires dealing with both
+possible definitions in all Hillas-style dl2 reconstructors.
+
+Advantages of computation in TelescopeFrame
+===========================================
+
+Computing the image parameters in ``TelescopeFrame`` – angular units on the sky –
+has the following advantages:
+
+* Parameters are easier to compare across different telescope types.
+* Pointing corrections can directly be applied in the conversion from ``CameraFrame``
+ to ``TelescopeFrame`` and are then automatically included in the image parameters,
+ which is much more straight forward than trying to correct image parameters that
+ are affected to different degrees after they have been computed.
+* Conversion from ``CameraFrame`` to ``TelescopeFrame`` will include any necessary
+ special handling of the curved cameras of dual mirror telescopes.
+
+
+Previous discussions
+====================
+* Issue discussing the removal of the camera frame image parameters: `#2061 `_
+* Original issue for introducing the computation of image parameters in telescope frame: `#1090 `_
+* Pull Request implementing image parameters in telescope frame, also setting it as the default: `#1591 `_
+* Adapting the reconstructors to also work with image parameters in telescope frame: `#1408 `_