Repository: cta-observatory/ctapipe
Pull request #2138 (2022-11-25 15:48:45)

Disp reconstruction
This is a continuation of #2064 after #1767 was merged.

For now I only implemented a simple weighted average to combine single telescope predictions into event predictions. The estimation of the error of this average seems to underestimate the error a bit, based on a small MC study I did. 
I will look into other ways to combine the telescope predictions in the future (including other error estimates), but for now this is better than no error at all, in my opinion. 

codecov[bot]: # [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
Base: **92.74**% // Head: **92.78**% // Increases project coverage by **`+0.03%`** :tada:
> Coverage data is based on head [(`596811b`)](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) compared to base [(`03d320c`)](https://codecov.io/gh/cta-observatory/ctapipe/commit/03d320cb9534e66ada513e2486794a262c54e9ed?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).
> Patch coverage: 94.29% of modified lines in pull request are covered.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #2138      +/-   ##
==========================================
+ Coverage   92.74%   92.78%   +0.03%     
==========================================
  Files         214      216       +2     
  Lines       17911    18385     +474     
==========================================
+ Hits        16612    17058     +446     
- Misses       1299     1327      +28     
```


| [Impacted Files](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | Coverage Δ | |
|---|---|---|
| [ctapipe/tools/process.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS90b29scy9wcm9jZXNzLnB5) | `87.87% <ø> (ø)` | |
| [ctapipe/tools/train\_energy\_regressor.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS90b29scy90cmFpbl9lbmVyZ3lfcmVncmVzc29yLnB5) | `85.29% <ø> (ø)` | |
| [ctapipe/tools/train\_particle\_classifier.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS90b29scy90cmFpbl9wYXJ0aWNsZV9jbGFzc2lmaWVyLnB5) | `86.58% <ø> (ø)` | |
| [ctapipe/tools/train\_disp\_reconstructor.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS90b29scy90cmFpbl9kaXNwX3JlY29uc3RydWN0b3IucHk=) | `87.05% <87.05%> (ø)` | |
| [ctapipe/conftest.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jb25mdGVzdC5weQ==) | `94.42% <91.66%> (-0.17%)` | :arrow_down: |
| [ctapipe/tools/apply\_models.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS90b29scy9hcHBseV9tb2RlbHMucHk=) | `95.04% <92.30%> (-1.34%)` | :arrow_down: |
| [ctapipe/reco/sklearn.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9yZWNvL3NrbGVhcm4ucHk=) | `94.60% <93.41%> (-0.89%)` | :arrow_down: |
| [ctapipe/reco/stereo\_combination.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9yZWNvL3N0ZXJlb19jb21iaW5hdGlvbi5weQ==) | `95.31% <96.15%> (+0.27%)` | :arrow_up: |
| [ctapipe/containers.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jb250YWluZXJzLnB5) | `100.00% <100.00%> (ø)` | |
| [ctapipe/coordinates/disp.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jb29yZGluYXRlcy9kaXNwLnB5) | `100.00% <100.00%> (ø)` | |
| ... and [7 more](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2138?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).

LukasBeiske: I seem to have missed that commit [71672d8](https://github.com/cta-observatory/ctapipe/pull/2138/commits/71672d849274133f0de27471da145e2f49e66815) breaks `apply-models` for the energy regressor (haven't tested the particle classifier yet) which is also why some tests are failing atm. 
However, I don't understand which change causes this error, especially since I did not actually change anything at the point in the code where the error occurs.
```
Traceback (most recent call last):
  File "/scratch/lbeiske/envs/cta/lib/python3.8/site-packages/tables/table.py", line 2192, in append
    wbufRA = np.rec.array(rows, dtype=self._v_dtype)
  File "/scratch/lbeiske/envs/cta/lib/python3.8/site-packages/numpy/core/records.py", line 1078, in array
    new = obj.view(dtype)
ValueError: When changing to a larger dtype, its size must be a divisor of the total size in bytes of the last axis of the array.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/net/nfshome/home/lbeiske/ctapipe/ctapipe/core/tool.py", line 358, in run
    self.start()
  File "/net/nfshome/home/lbeiske/ctapipe/ctapipe/tools/apply_models.py", line 178, in start
    self._apply(reconstructor)
  File "/net/nfshome/home/lbeiske/ctapipe/ctapipe/tools/apply_models.py", line 272, in _apply
    write_table(
  File "/net/nfshome/home/lbeiske/ctapipe/ctapipe/io/astropy_helpers.py", line 211, in write_table
    h5_table.append(table.as_array())
  File "/scratch/lbeiske/envs/cta/lib/python3.8/site-packages/tables/table.py", line 2194, in append
    raise ValueError("rows parameter cannot be converted into a "
ValueError: rows parameter cannot be converted into a recarray object compliant with table '/dl2/event/telescope/energy/RandomForestRegressor/tel_001 (Table(18158,)fletcher32, shuffle, blosc:zstd(5)) '''. The error was: <When changing to a larger dtype, its size must be a divisor of the total size in bytes of the last axis of the array.>
```
nbiederbeck: We found the issue in https://github.com/cta-observatory/ctapipe/blob/f789773f94bc5f33b33258dab13b5d4bd6c467dd/ctapipe/io/tableloader.py#L326-L329

```python
(Pdb) table['obs_id'].dtype
dtype('int32')
(Pdb) observation_table['obs_id'].dtype
dtype('uint64')
```

The join makes a float64 out of this.
maxnoe: @LukasBeiske Could you show some performance plots if you think this is ready?
LukasBeiske: Sure. [cv_disp.pdf](https://github.com/cta-observatory/ctapipe/files/10186763/cv_disp.pdf) these are the performance plots for the model I am currently using. 

maxnoe: Why is core-y so important? That's confusing me a bit.

Could you show the *single* LST performance? I.e. no stereo features for LST?
LukasBeiske: For mono I used this config:
```yaml
TrainDispReconstructor:
  CrossValidator:
    n_cross_validations: 5

  random_seed: 42

  DispReconstructor:
    norm_cls: RandomForestRegressor
    norm_config:
      n_estimators: 100
      max_features: "sqrt"
      min_samples_leaf: 0.00001
      n_jobs: 40

    log_target: True

    sign_cls: RandomForestClassifier
    sign_config:
      n_estimators: 100
      max_features: "sqrt"
      min_samples_split: 0.00001
      n_jobs: 40

    features:
      - hillas_intensity
      - hillas_length
      - hillas_width
      - hillas_skewness
      - hillas_kurtosis
      - timing_slope
      - timing_deviation
      - peak_time_std
      - concentration_cog
      - leakage_intensity_width_2
      - morphology_n_pixels
      - scaled_length
      - scaled_width
      - RandomForestRegressor_energy
      - log_abs_timing_slope

    QualityQuery:
      quality_criteria:
        - ["enough intensity", "hillas_intensity > 50"]
        - ["Positive width", "hillas_width > 0"]
        - ["enough pixels", "morphology_n_pixels > 3"]
        - ["not clipped", "leakage_intensity_width_2 < 0.5"]

    FeatureGenerator:
      features:
        - ["log_abs_timing_slope", "log(abs(timing_slope))"]
```
And the performance plots still look pretty good imo (sign performs slightly worse): [cv_disp_mono.pdf](https://github.com/cta-observatory/ctapipe/files/10187475/cv_disp_mono.pdf)

maxnoe: @LukasBeiske Please rebase to fix CI build

Patch:
diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index d31c4f05..ebf56d83 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -45,6 +45,7 @@ __all__ = [
     "ReconstructedContainer",
     "ReconstructedEnergyContainer",
     "ReconstructedGeometryContainer",
+    "DispContainer",
     "SimulatedCameraContainer",
     "SimulatedShowerContainer",
     "SimulatedShowerDistribution",
@@ -754,6 +755,11 @@ class ReconstructedGeometryContainer(Container):
     alt_uncert = Field(nan * u.deg, "reconstructed altitude uncertainty", unit=u.deg)
     az = Field(nan * u.deg, "reconstructed azimuth", unit=u.deg)
     az_uncert = Field(nan * u.deg, "reconstructed azimuth uncertainty", unit=u.deg)
+    ang_distance_uncert = Field(
+        nan * u.deg,
+        "uncertainty radius of reconstructed altitude-azimuth position",
+        unit=u.deg,
+    )
     core_x = Field(
         nan * u.m, "reconstructed x coordinate of the core position", unit=u.m
     )
@@ -843,6 +849,17 @@ class ParticleClassificationContainer(Container):
     telescopes = Field(None, "Telescopes used if stereo, or None if Mono")
 
 
+class DispContainer(Container):
+    """
+    Standard output of disp reconstruction algorithms for origin reconstruction
+    """
+
+    default_prefix = "disp_parameter"
+
+    norm = Field(nan * u.deg, "reconstructed value for disp", unit=u.deg)
+    is_valid = Field(False, "true if the predictions are valid")
+
+
 class ReconstructedContainer(Container):
     """Reconstructed shower info from multiple algorithms"""
 
@@ -875,6 +892,10 @@ class TelescopeReconstructedContainer(ReconstructedContainer):
         default_factory=partial(Map, TelescopeImpactParameterContainer),
         description="map of algorithm to impact parameter info",
     )
+    disp = Field(
+        default_factory=partial(Map, DispContainer),
+        description="map of algorithm to reconstructed disp parameters",
+    )
 
 
 class DL2Container(Container):
diff --git a/ctapipe/reco/__init__.py b/ctapipe/reco/__init__.py
index e6f48adc..54bb2a17 100644
--- a/ctapipe/reco/__init__.py
+++ b/ctapipe/reco/__init__.py
@@ -7,7 +7,12 @@ from .hillas_reconstructor import HillasReconstructor
 from .impact import ImPACTReconstructor
 from .reconstructor import GeometryReconstructor, Reconstructor
 from .shower_processor import ShowerProcessor
-from .sklearn import CrossValidator, EnergyRegressor, ParticleClassifier
+from .sklearn import (
+    CrossValidator,
+    DispReconstructor,
+    EnergyRegressor,
+    ParticleClassifier,
+)
 from .stereo_combination import StereoCombiner, StereoMeanCombiner
 
 __all__ = [
@@ -19,6 +24,7 @@ __all__ = [
     "HillasIntersection",
     "EnergyRegressor",
     "ParticleClassifier",
+    "DispReconstructor",
     "StereoCombiner",
     "StereoMeanCombiner",
     "CrossValidator",
diff --git a/ctapipe/reco/preprocessing.py b/ctapipe/reco/preprocessing.py
index a60ae45c..d5adcf20 100644
--- a/ctapipe/reco/preprocessing.py
+++ b/ctapipe/reco/preprocessing.py
@@ -6,11 +6,17 @@ These functions are mainly taken from https://github.com/fact-project/aict-tools
 and adapted to work with astropy Tables instead of pandas dataframes
 """
 import logging
+import warnings
+from typing import List
 
+import astropy.units as u
 import numpy as np
+from astropy.coordinates import AltAz
 from astropy.table import QTable, Table
 from numpy.lib.recfunctions import structured_to_unstructured
 
+from ctapipe.coordinates import MissingFrameAttributeWarning, TelescopeFrame
+
 from ..containers import ArrayEventContainer
 
 LOG = logging.getLogger(__name__)
@@ -20,6 +26,9 @@ __all__ = [
     "check_valid_rows",
     "collect_features",
     "table_to_float",
+    "table_to_X",
+    "horizontal_to_telescope",
+    "telescope_to_horizontal",
 ]
 
 
@@ -48,6 +57,17 @@ def check_valid_rows(table: Table, warn=True, log=LOG) -> np.ndarray:
     return valid
 
 
+def table_to_X(table: Table, features: List[str], log=LOG):
+    """
+    Extract features as numpy ndarray to be given to sklearn from input table
+    dropping all events for which one or more training features are nan
+    """
+    feature_table = table[features]
+    valid = check_valid_rows(feature_table, log=log)
+    X = table_to_float(feature_table[valid])
+    return X, valid
+
+
 def collect_features(
     event: ArrayEventContainer, tel_id: int, subarray_table=None
 ) -> Table:
@@ -102,3 +122,33 @@ def collect_features(
         features.update(subarray_table.loc[tel_id])
 
     return Table({k: [v] for k, v in features.items()})
+
+
+@u.quantity_input(alt=u.deg, az=u.deg, pointing_alt=u.deg, pointing_az=u.deg)
+def horizontal_to_telescope(alt, az, pointing_alt, pointing_az):
+    """Transform coordinates from horizontal coordinates into TelescopeFrame"""
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", MissingFrameAttributeWarning)
+
+        horizontal_coord = AltAz(alt=alt, az=az)
+        pointing = AltAz(alt=pointing_alt, az=pointing_az)
+        tel_frame = TelescopeFrame(telescope_pointing=pointing)
+
+        tel_coord = horizontal_coord.transform_to(tel_frame)
+
+    return tel_coord.fov_lon.to(u.deg), tel_coord.fov_lat.to(u.deg)
+
+
+@u.quantity_input(lon=u.deg, lat=u.deg, pointing_alt=u.deg, pointing_az=u.deg)
+def telescope_to_horizontal(lon, lat, pointing_alt, pointing_az):
+    """Transform coordinates from TelescopeFrame into horizontal coordinates"""
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", MissingFrameAttributeWarning)
+
+        pointing = AltAz(alt=pointing_alt, az=pointing_az)
+        tel_coord = TelescopeFrame(
+            fov_lon=lon, fov_lat=lat, telescope_pointing=pointing
+        )
+        horizontal_coord = tel_coord.transform_to(AltAz())
+
+    return horizontal_coord.alt.to(u.deg), horizontal_coord.az.to(u.deg)
diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py
index 5fd447ad..42ff43d5 100644
--- a/ctapipe/reco/sklearn.py
+++ b/ctapipe/reco/sklearn.py
@@ -4,22 +4,27 @@ Component Wrappers around sklearn models
 import pathlib
 from abc import abstractmethod
 from collections import defaultdict
+from typing import Tuple
 
 import astropy.units as u
 import joblib
 import numpy as np
+from astropy.coordinates import AltAz
 from astropy.table import QTable, Table, vstack
 from astropy.utils.decorators import lazyproperty
-from sklearn.metrics import r2_score, roc_auc_score
+from sklearn.metrics import accuracy_score, r2_score, roc_auc_score
 from sklearn.model_selection import KFold, StratifiedKFold
 from sklearn.utils import all_estimators
 from tqdm import tqdm
 
 from ..containers import (
     ArrayEventContainer,
+    DispContainer,
     ParticleClassificationContainer,
     ReconstructedEnergyContainer,
+    ReconstructedGeometryContainer,
 )
+from ..coordinates import TelescopeFrame
 from ..core import Component, FeatureGenerator, Provenance, QualityQuery
 from ..core.traits import (
     Bool,
@@ -34,7 +39,7 @@ from ..core.traits import (
     Unicode,
 )
 from ..io import write_table
-from .preprocessing import check_valid_rows, collect_features, table_to_float
+from .preprocessing import collect_features, table_to_X, telescope_to_horizontal
 from .reconstructor import Reconstructor
 from .stereo_combination import StereoCombiner
 from .utils import add_defaults_and_meta
@@ -42,9 +47,10 @@ from .utils import add_defaults_and_meta
 __all__ = [
     "SKLearnReconstructor",
     "SKLearnRegressionReconstructor",
-    "SKLearnClassficationReconstructor",
+    "SKLearnClassificationReconstructor",
     "EnergyRegressor",
     "ParticleClassifier",
+    "DispReconstructor",
     "CrossValidator",
 ]
 
@@ -127,7 +133,6 @@ class SKLearnReconstructor(Reconstructor):
                 self.prefix = self.model_cls
 
             super().__init__(subarray, **kwargs)
-            self.subarray = subarray
             self.feature_generator = FeatureGenerator(parent=self)
             self.quality_query = MLQualityQuery(parent=self)
 
@@ -221,15 +226,6 @@ class SKLearnReconstructor(Reconstructor):
     def _new_model(self):
         return SUPPORTED_MODELS[self.model_cls](**self.model_config)
 
-    def _table_to_X(self, table):
-        """
-        Extract features as numpy ndarray to be given to sklearn from input table
-        """
-        feature_table = table[self.features]
-        valid = check_valid_rows(feature_table, log=self.log)
-        X = table_to_float(feature_table[valid])
-        return X, valid
-
     def _table_to_y(self, table, mask=None):
         """
         Extract target values as numpy array from input table
@@ -246,7 +242,7 @@ class SKLearnReconstructor(Reconstructor):
         """
         self._models[key] = self._new_model()
 
-        X, valid = self._table_to_X(table)
+        X, valid = table_to_X(table, self.features, self.log)
         self.unit = table[self.target].unit
         y = self._table_to_y(table, mask=valid)
         self._models[key].fit(X, y)
@@ -281,7 +277,7 @@ class SKLearnRegressionReconstructor(SKLearnReconstructor):
                 f"No model available for key {key},"
                 f" available models: {self._models.keys()}"
             )
-        X, valid = self._table_to_X(table)
+        X, valid = table_to_X(table, self.features, self.log)
         n_outputs = getattr(self._models[key], "n_outputs_", 1)
 
         if n_outputs > 1:
@@ -314,7 +310,7 @@ class SKLearnRegressionReconstructor(SKLearnReconstructor):
         return y
 
 
-class SKLearnClassficationReconstructor(SKLearnReconstructor):
+class SKLearnClassificationReconstructor(SKLearnReconstructor):
     """
     Base class for classification tasks
     """
@@ -323,7 +319,7 @@ class SKLearnClassficationReconstructor(SKLearnReconstructor):
         SUPPORTED_CLASSIFIERS.keys(),
         default_value=None,
         allow_none=True,
-        help="Which scikit-learn regression model to use.",
+        help="Which scikit-learn classification model to use.",
     ).tag(config=True)
 
     invalid_class = Integer(
@@ -346,7 +342,7 @@ class SKLearnClassficationReconstructor(SKLearnReconstructor):
                 f" available models: {self._models.keys()}"
             )
 
-        X, valid = self._table_to_X(table)
+        X, valid = table_to_X(table, self.features, self.log)
         n_outputs = getattr(self._models[key], "n_outputs_", 1)
 
         if n_outputs > 1:
@@ -367,7 +363,7 @@ class SKLearnClassficationReconstructor(SKLearnReconstructor):
                 f" available models: {self._models.keys()}"
             )
 
-        X, valid = self._table_to_X(table)
+        X, valid = table_to_X(table, self.features, self.log)
 
         n_classes = getattr(self._models[key], "n_classes_", 2)
         n_rows = len(table)
@@ -456,7 +452,7 @@ class EnergyRegressor(SKLearnRegressionReconstructor):
         return result
 
 
-class ParticleClassifier(SKLearnClassficationReconstructor):
+class ParticleClassifier(SKLearnClassificationReconstructor):
     """
     Predict dl2 particle classification
     """
@@ -519,6 +515,320 @@ class ParticleClassifier(SKLearnClassficationReconstructor):
         return result
 
 
+class DispReconstructor(Reconstructor):
+    """
+    Predict absolute value and sign for disp origin reconstruction for each telescope.
+    """
+
+    property = "geometry"
+    prefix = Unicode(default_value="disp", allow_none=False).tag(config=True)
+    features = List(Unicode(), help="Features to use for both models").tag(config=True)
+    target = "true_norm"
+    log_target = Bool(
+        default_value=False,
+        help="If True, the model is trained to predict the natural logarithm of the absolute value.",
+    ).tag(config=True)
+
+    norm_config = Dict({}, help="kwargs for the sklearn regressor").tag(config=True)
+    norm_cls = Enum(
+        SUPPORTED_REGRESSORS.keys(),
+        default_value=None,
+        allow_none=True,
+        help="Which scikit-learn regression model to use.",
+    ).tag(config=True)
+
+    sign_config = Dict({}, help="kwargs for the sklearn classifier").tag(config=True)
+    sign_cls = Enum(
+        SUPPORTED_CLASSIFIERS.keys(),
+        default_value=None,
+        allow_none=True,
+        help="Which scikit-learn classification model to use.",
+    ).tag(config=True)
+
+    stereo_combiner_cls = ComponentName(
+        StereoCombiner,
+        default_value="StereoMeanCombiner",
+        help="Which stereo combination method to use",
+    ).tag(config=True)
+
+    load_path = Path(
+        default_value=None,
+        allow_none=True,
+        help="If given, load serialized model from this path",
+    ).tag(config=True)
+
+    def __init__(self, subarray=None, models=None, **kwargs):
+        # Run the Component __init__ first to handle the configuration
+        # and make `self.load_path` available
+        Component.__init__(self, **kwargs)
+
+        if self.load_path is None:
+            if self.norm_cls is None or self.sign_cls is None:
+                raise TraitError(
+                    "Must provide `norm_cls` and `sign_cls` if not loading from file"
+                )
+
+            if subarray is None:
+                raise TypeError(
+                    "__init__() missing 1 required positional argument: 'subarray'"
+                )
+
+            super().__init__(subarray, **kwargs)
+            self.quality_query = MLQualityQuery(parent=self)
+            self.feature_generator = FeatureGenerator(parent=self)
+
+            # to verify settings
+            self._new_models()
+
+            self._models = {} if models is None else models
+            self.unit = None
+            self.stereo_combiner = StereoCombiner.from_name(
+                self.stereo_combiner_cls,
+                prefix=self.prefix,
+                property=self.property,
+                parent=self,
+            )
+        else:
+            loaded = self.read(self.load_path)
+            if (
+                subarray is not None
+                and loaded.subarray.telescope_types != subarray.telescope_types
+            ):
+                self.log.warning(
+                    "Supplied subarray has different telescopes than subarray loaded from file"
+                )
+            self.__dict__.update(loaded.__dict__)
+            self.subarray = subarray
+
+    def _new_models(self):
+        norm_regressor = SUPPORTED_REGRESSORS[self.norm_cls](**self.norm_config)
+        sign_classifier = SUPPORTED_CLASSIFIERS[self.sign_cls](**self.sign_config)
+        return norm_regressor, sign_classifier
+
+    def _table_to_y(self, table, mask=None):
+        """
+        Extract target values as numpy array from input table
+        """
+        # make sure we use the unit that was used during training
+        if self.unit is not None:
+            norm = table[mask][self.target].quantity.to_value(self.unit)
+        else:
+            norm = np.array(table[self.target][mask])
+
+        abs_norm = np.abs(norm)
+        sign_norm = np.sign(norm)
+
+        if self.log_target:
+            abs_norm = np.log(abs_norm)
+
+        return abs_norm, sign_norm
+
+    def fit(self, key, table):
+        """
+        Create and fit new models for ``key`` using the data in ``table``.
+        """
+        self._models[key] = self._new_models()
+
+        X, valid = table_to_X(table, self.features, self.log)
+        self.unit = table[self.target].unit
+        norm, sign = self._table_to_y(table, mask=valid)
+        self._models[key][0].fit(X, norm)
+        self._models[key][1].fit(X, sign)
+
+    def write(self, path, overwrite=False):
+        path = pathlib.Path(path)
+
+        if path.exists() and not overwrite:
+            raise IOError(f"Path {path} exists and overwrite=False")
+
+        with path.open("wb") as f:
+            Provenance().add_output_file(path, role="ml-models")
+            joblib.dump(self, f, compress=True)
+
+    @classmethod
+    def read(cls, path, **kwargs):
+        with open(path, "rb") as f:
+            instance = joblib.load(f)
+
+        for attr, value in kwargs.items():
+            setattr(instance, attr, value)
+
+        if not isinstance(instance, cls):
+            raise TypeError(
+                f"{path} did not contain an instance of {cls}, got {instance}"
+            )
+
+        Provenance().add_input_file(path, role="ml-models")
+        return instance
+
+    @lazyproperty
+    def instrument_table(self):
+        return self.subarray.to_table("joined")
+
+    def _predict(self, key, table):
+        if key not in self._models:
+            raise KeyError(
+                f"No model available for key {key},"
+                f" available models: {self._models.keys()}"
+            )
+        X, valid = table_to_X(table, self.features, self.log)
+        prediction = np.full(len(table), np.nan)
+
+        if np.any(valid):
+            valid_norms = self._models[key][0].predict(X)
+
+            if self.log_target:
+                prediction[valid] = np.exp(valid_norms)
+            else:
+                prediction[valid] = valid_norms
+
+            prediction[valid] *= self._models[key][1].predict(X)
+
+        if self.unit is not None:
+            prediction = u.Quantity(prediction, self.unit, copy=False)
+
+        return prediction, valid
+
+    def __call__(self, event: ArrayEventContainer) -> None:
+        """Event-wise prediction for the EventSource-Loop.
+
+        Fills the event.dl2.tel[tel_id].disp[prefix] container
+        and event.dl2.tel[tel_id].geometry[prefix] container.
+
+        Parameters
+        ----------
+        event: ArrayEventContainer
+        """
+        for tel_id in event.trigger.tels_with_trigger:
+            table = collect_features(event, tel_id, self.instrument_table)
+            table = self.feature_generator(table)
+
+            passes_quality_checks = self.quality_query.get_table_mask(table)[0]
+
+            if passes_quality_checks:
+                disp, valid = self._predict(self.subarray.tel[tel_id], table)
+                disp_container = DispContainer(
+                    norm=disp[0],
+                    is_valid=valid[0],
+                )
+
+                if valid:
+                    hillas = event.dl1.tel[tel_id].parameters.hillas
+                    psi = hillas.psi.to_value(u.rad)
+
+                    fov_lon = hillas.fov_lon + disp[0] * np.cos(psi)
+                    fov_lat = hillas.fov_lat + disp[0] * np.sin(psi)
+                    altaz = TelescopeFrame(
+                        fov_lon=fov_lon,
+                        fov_lat=fov_lat,
+                        telescope_pointing=AltAz(
+                            alt=event.pointing.tel[tel_id].altitude,
+                            az=event.pointing.tel[tel_id].azimuth,
+                        ),
+                    ).transform_to(AltAz())
+
+                    altaz_container = ReconstructedGeometryContainer(
+                        alt=altaz.alt, az=altaz.az, is_valid=True
+                    )
+
+                else:
+                    altaz_container = ReconstructedGeometryContainer(
+                        alt=u.Quantity(np.nan, u.deg, copy=False),
+                        az=u.Quantity(np.nan, u.deg, copy=False),
+                        is_valid=False,
+                    )
+            else:
+                disp_container = DispContainer(
+                    norm=u.Quantity(np.nan, self.unit),
+                    is_valid=False,
+                )
+                altaz_container = ReconstructedGeometryContainer(
+                    alt=u.Quantity(np.nan, u.deg, copy=False),
+                    az=u.Quantity(np.nan, u.deg, copy=False),
+                    is_valid=False,
+                )
+
+            disp_container.prefix = f"{self.prefix}_parameter"
+            altaz_container.prefix = f"{self.prefix}_tel"
+            event.dl2.tel[tel_id].disp[self.prefix] = disp_container
+            event.dl2.tel[tel_id].geometry[self.prefix] = altaz_container
+
+        self.stereo_combiner(event)
+
+    def predict_table(self, key, table: Table) -> Tuple[Table, Table]:
+        """Predict on a table of events
+
+        Parameters
+        ----------
+        table : `~astropy.table.Table`
+            Table of features
+
+        Returns
+        -------
+        disp_table : `~astropy.table.Table`
+            Table with disp predictions, matches the corresponding
+            container definition
+        altaz_table : `~astropy.table.Table`
+            Table with resulting predictions of horizontal coordinates
+        """
+        table = self.feature_generator(table)
+
+        n_rows = len(table)
+        disp = u.Quantity(np.full(n_rows, np.nan), self.unit, copy=False)
+        is_valid = np.full(n_rows, False)
+
+        valid = self.quality_query.get_table_mask(table)
+        disp[valid], is_valid[valid] = self._predict(key, table[valid])
+
+        disp_result = Table(
+            {
+                f"{self.prefix}_parameter_norm": disp,
+                f"{self.prefix}_parameter_is_valid": is_valid,
+            }
+        )
+        add_defaults_and_meta(
+            disp_result,
+            DispContainer,
+            prefix=f"{self.prefix}_parameter",
+            stereo=False,
+        )
+
+        psi = table["hillas_psi"].quantity.to_value(u.rad)
+        fov_lon = table["hillas_fov_lon"].quantity + disp * np.cos(psi)
+        fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi)
+
+        # FIXME: Assume constant and parallel pointing for each run
+        self.log.warning("Assuming constant and parallel pointing for each run")
+        alt, az = telescope_to_horizontal(
+            lon=fov_lon,
+            lat=fov_lat,
+            pointing_alt=table["subarray_pointing_lat"],
+            pointing_az=table["subarray_pointing_lon"],
+        )
+
+        altaz_result = Table(
+            {
+                f"{self.prefix}_tel_alt": alt,
+                f"{self.prefix}_tel_az": az,
+                f"{self.prefix}_tel_is_valid": is_valid,
+            }
+        )
+        add_defaults_and_meta(
+            altaz_result,
+            ReconstructedGeometryContainer,
+            prefix=self.prefix,
+            stereo=False,
+        )
+
+        return disp_result, altaz_result
+
+    def __getstate__(self):
+        state = self.__dict__.copy()
+        state["_trait_values"]["parent"] = None
+        state["_trait_notifiers"] = {}
+        return state
+
+
 class CrossValidator(Component):
     """Class to train sklearn based reconstructors in a cross validation"""
 
@@ -546,12 +856,15 @@ class CrossValidator(Component):
         self.model_component = model_component
         self.rng = np.random.default_rng(self.rng_seed)
 
-        if isinstance(self.model_component, SKLearnClassficationReconstructor):
+        if isinstance(self.model_component, SKLearnClassificationReconstructor):
             self.cross_validate = self._cross_validate_classification
             self.split_data = StratifiedKFold
         elif isinstance(self.model_component, SKLearnRegressionReconstructor):
             self.cross_validate = self._cross_validate_regressor
             self.split_data = KFold
+        elif isinstance(self.model_component, DispReconstructor):
+            self.cross_validate = self._cross_validate_disp
+            self.split_data = StratifiedKFold
         else:
             raise KeyError(
                 "Unsupported Model of type %s supplied", self.model_component
@@ -577,7 +890,11 @@ class CrossValidator(Component):
             random_state=self.rng.integers(0, 2**31 - 1),
         )
 
-        cv_it = kfold.split(table, table[self.model_component.target])
+        if isinstance(self.model_component, DispReconstructor):
+            cv_it = kfold.split(table, np.sign(table[self.model_component.target]))
+        else:
+            cv_it = kfold.split(table, table[self.model_component.target])
+
         for fold, (train_indices, test_indices) in enumerate(
             tqdm(
                 cv_it,
@@ -587,6 +904,7 @@ class CrossValidator(Component):
         ):
             train = table[train_indices]
             test = table[test_indices]
+
             cv_prediction, truth, metrics = self.cross_validate(
                 telescope_type, train, test
             )
@@ -597,6 +915,7 @@ class CrossValidator(Component):
                         "tel_type": [str(telescope_type)] * len(truth),
                         "prediction": cv_prediction,
                         "truth": truth,
+                        "true_energy": test["true_energy"],
                     }
                 )
             )
@@ -635,6 +954,15 @@ class CrossValidator(Component):
         roc_auc = roc_auc_score(truth, prediction)
         return prediction, truth, {"ROC AUC": roc_auc}
 
+    def _cross_validate_disp(self, telescope_type, train, test):
+        models = self.model_component
+        models.fit(telescope_type, train)
+        prediction, _ = models._predict(telescope_type, test)
+        truth = test[models.target]
+        r2 = r2_score(np.abs(truth), np.abs(prediction))
+        accuracy = accuracy_score(np.sign(truth), np.sign(prediction))
+        return prediction, truth, {"R^2": r2, "accuracy": accuracy}
+
     def write(self, overwrite=False):
         if self.output_path.exists() and not overwrite:
             raise IOError(f"Path {self.output_path} exists and overwrite=False")
@@ -642,5 +970,5 @@ class CrossValidator(Component):
         Provenance().add_output_file(self.output_path, role="ml-cross-validation")
         for tel_type, results in self.cv_predictions.items():
             write_table(
-                results, self.output_path, f"cv_predictions_{tel_type}", overwrite=True
+                results, self.output_path, f"/cv_predictions_{tel_type}", overwrite=True
             )
diff --git a/ctapipe/reco/stereo_combination.py b/ctapipe/reco/stereo_combination.py
index 392711bc..6789ea33 100644
--- a/ctapipe/reco/stereo_combination.py
+++ b/ctapipe/reco/stereo_combination.py
@@ -2,6 +2,7 @@ from abc import abstractmethod
 
 import astropy.units as u
 import numpy as np
+from astropy.coordinates import AltAz, CartesianRepresentation, SphericalRepresentation
 from astropy.table import Table
 
 from ctapipe.core import Component, Container
@@ -11,12 +12,14 @@ from ..containers import (
     ArrayEventContainer,
     ParticleClassificationContainer,
     ReconstructedEnergyContainer,
+    ReconstructedGeometryContainer,
 )
 from .utils import add_defaults_and_meta
 
 _containers = {
     "energy": ReconstructedEnergyContainer,
     "classification": ParticleClassificationContainer,
+    "geometry": ReconstructedGeometryContainer,
 }
 
 __all__ = [
@@ -102,10 +105,10 @@ class StereoMeanCombiner(StereoCombiner):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
-        if self.property not in {"energy", "classification"}:
+        if self.property not in {"energy", "classification", "geometry"}:
             raise NotImplementedError(
                 f"Cannot combine {self.property}."
-                "Implemented are: energy, classification"
+                "Implemented are: energy, classification, geometry"
             )
 
     def _calculate_weights(self, data):
@@ -205,6 +208,55 @@ class StereoMeanCombiner(StereoCombiner):
         )
         event.dl2.stereo.classification[self.prefix] = container
 
+    def _combine_disp(self, event):
+        ids = []
+        alt_values = []
+        az_values = []
+        weights = []
+
+        for tel_id, dl2 in event.dl2.tel.items():
+            mono = dl2.geometry[self.prefix]
+            if mono.is_valid:
+                alt_values.append(mono.alt)
+                az_values.append(mono.az)
+                dl1 = event.dl1.tel[tel_id].parameters
+                weights.append(self._calculate_weights(dl1) if dl1 else 1)
+                ids.append(tel_id)
+
+        if len(alt_values) > 0:  # by construction len(alt_values) == len(az_values)
+            coord = AltAz(alt=alt_values, az=az_values)
+            # https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution#Mean_direction
+            mono_x, mono_y, mono_z = coord.cartesian.get_xyz()
+            stereo_x = np.average(mono_x, weights=weights)
+            stereo_y = np.average(mono_y, weights=weights)
+            stereo_z = np.average(mono_z, weights=weights)
+
+            mean_cartesian = CartesianRepresentation(x=stereo_x, y=stereo_y, z=stereo_z)
+            mean_spherical = mean_cartesian.represent_as(SphericalRepresentation)
+            mean_altaz = AltAz(mean_spherical)
+
+            # https://en.wikipedia.org/wiki/Directional_statistics#Measures_of_location_and_spread
+            r = mean_spherical.distance.to_value()
+            std = np.sqrt(-2 * np.log(r))
+
+            valid = True
+        else:
+            mean_altaz = AltAz(
+                alt=u.Quantity(np.nan, u.deg, copy=False),
+                az=u.Quantity(np.nan, u.deg, copy=False),
+            )
+            std = np.nan
+            valid = False
+
+        event.dl2.stereo.geometry[self.prefix] = ReconstructedGeometryContainer(
+            alt=mean_altaz.alt,
+            az=mean_altaz.az,
+            ang_distance_uncert=u.Quantity(np.rad2deg(std), u.deg, copy=False),
+            telescopes=ids,
+            is_valid=valid,
+            prefix=self.prefix,
+        )
+
     def __call__(self, event: ArrayEventContainer) -> None:
         """
         Calculate the mean prediction for a single array event.
@@ -215,6 +267,9 @@ class StereoMeanCombiner(StereoCombiner):
         elif self.property == "classification":
             self._combine_classification(event)
 
+        elif self.property == "geometry":
+            self._combine_disp(event)
+
     def predict_table(self, mono_predictions: Table) -> Table:
         """
         Calculates the (array-)event-wise mean.
@@ -294,6 +349,55 @@ class StereoMeanCombiner(StereoCombiner):
             )
             stereo_table[f"{self.prefix}_is_valid"] = np.isfinite(stereo_energy)
             stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan
+
+        elif self.property == "geometry":
+            if len(valid_predictions) > 0:
+                coord = AltAz(
+                    alt=valid_predictions[f"{prefix}_alt"],
+                    az=valid_predictions[f"{prefix}_az"],
+                )
+                # https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution#Mean_direction
+                mono_x, mono_y, mono_z = coord.cartesian.get_xyz()
+
+                stereo_x = _weighted_mean_ufunc(
+                    mono_x, weights, n_array_events, indices[valid]
+                )
+                stereo_y = _weighted_mean_ufunc(
+                    mono_y, weights, n_array_events, indices[valid]
+                )
+                stereo_z = _weighted_mean_ufunc(
+                    mono_z, weights, n_array_events, indices[valid]
+                )
+
+                mean_cartesian = CartesianRepresentation(
+                    x=stereo_x, y=stereo_y, z=stereo_z
+                )
+                mean_spherical = mean_cartesian.represent_as(SphericalRepresentation)
+                mean_altaz = AltAz(mean_spherical)
+
+                # https://en.wikipedia.org/wiki/Directional_statistics#Measures_of_location_and_spread
+                r = mean_spherical.distance.to_value()
+                std = np.sqrt(-2 * np.log(r))
+            else:
+                mean_altaz = AltAz(
+                    alt=u.Quantity(np.full(n_array_events, np.nan), u.deg, copy=False),
+                    az=u.Quantity(np.full(n_array_events, np.nan), u.deg, copy=False),
+                )
+                std = np.full(n_array_events, np.nan)
+
+            stereo_table[f"{self.prefix}_alt"] = mean_altaz.alt.to(u.deg)
+            stereo_table[f"{self.prefix}_az"] = mean_altaz.az.to(u.deg)
+
+            stereo_table[f"{self.prefix}_ang_distance_uncert"] = u.Quantity(
+                np.rad2deg(std), u.deg, copy=False
+            )
+
+            stereo_table[f"{self.prefix}_is_valid"] = np.logical_and(
+                np.isfinite(stereo_table[f"{self.prefix}_alt"]),
+                np.isfinite(stereo_table[f"{self.prefix}_az"]),
+            )
+            stereo_table[f"{self.prefix}_goodness_of_fit"] = np.nan
+
         else:
             raise NotImplementedError()
 
diff --git a/ctapipe/resources/train_disp_reconstructor.yaml b/ctapipe/resources/train_disp_reconstructor.yaml
new file mode 100644
index 00000000..de10603b
--- /dev/null
+++ b/ctapipe/resources/train_disp_reconstructor.yaml
@@ -0,0 +1,68 @@
+# ======================================================================
+# ctapipe-train-disp-reconstructor config file.
+# version: VERSION
+#
+# Configuration for training machine-learning models
+# ======================================================================
+TrainDispReconstructor:
+  CrossValidator:
+    n_cross_validations: 5
+
+  DispReconstructor:
+    norm_cls: ExtraTreesRegressor
+    norm_config:
+      n_estimators: 10
+      max_depth: 10
+      n_jobs: -1
+
+    log_target: True
+
+    sign_cls: ExtraTreesClassifier
+    sign_config:
+      n_estimators: 10
+      max_depth: 10
+      n_jobs: -1
+
+    features:
+      - hillas_intensity
+      - hillas_r
+      - hillas_length
+      - hillas_width
+      - hillas_skewness
+      - hillas_kurtosis
+      - timing_slope
+      - timing_deviation
+      - leakage_intensity_width_1
+      - leakage_intensity_width_2
+      - leakage_pixels_width_1
+      - leakage_pixels_width_2
+      - concentration_cog
+      - concentration_core
+      - morphology_n_islands
+      - morphology_n_large_islands
+      - morphology_n_medium_islands
+      - morphology_n_pixels
+      - morphology_n_small_islands
+      - peak_time_std
+      - peak_time_skewness
+      - peak_time_kurtosis
+      - intensity_max
+      - intensity_mean
+      - intensity_std
+      - intensity_skewness
+      - intensity_kurtosis
+      - area
+      - ExtraTreesRegressor_energy
+      - ExtraTreesRegressor_tel_energy
+
+    QualityQuery:
+      quality_criteria:
+        - ["HillasValid", "HillasReconstructor_is_valid"]
+        - ["enough intensity", "hillas_intensity > 50"]
+        - ["Positive width", "hillas_width > 0"]
+        - ["enough pixels", "morphology_n_pixels > 3"]
+        - ["not clipped", "leakage_intensity_width_2 < 0.5"]
+
+    FeatureGenerator:
+      features:
+        - ["area", "hillas_width * hillas_length"]
\ No newline at end of file
diff --git a/ctapipe/tools/apply_models.py b/ctapipe/tools/apply_models.py
index c56bc4a2..ccb1501d 100644
--- a/ctapipe/tools/apply_models.py
+++ b/ctapipe/tools/apply_models.py
@@ -14,7 +14,12 @@ from ctapipe.io import TableLoader, write_table
 from ctapipe.io.astropy_helpers import read_table
 from ctapipe.io.tableio import TelListToMaskTransform
 from ctapipe.io.tableloader import _join_subarray_events
-from ctapipe.reco import EnergyRegressor, ParticleClassifier, StereoCombiner
+from ctapipe.reco import (
+    DispReconstructor,
+    EnergyRegressor,
+    ParticleClassifier,
+    StereoCombiner,
+)
 
 __all__ = [
     "ApplyModels",
@@ -76,6 +81,14 @@ class ApplyModels(Tool):
         help="Input path for the trained ParticleClassifier",
     ).tag(config=True)
 
+    disp_reconstructor_path = Path(
+        default_value=None,
+        allow_none=True,
+        exists=True,
+        directory_ok=False,
+        help="Input path for the trained DispReconstructor",
+    ).tag(config=True)
+
     chunk_size = Integer(
         default_value=100000,
         allow_none=True,
@@ -86,6 +99,7 @@ class ApplyModels(Tool):
         ("i", "input"): "ApplyModels.input_url",
         "energy-regressor": "ApplyModels.energy_regressor_path",
         "particle-classifier": "ApplyModels.particle_classifier_path",
+        "disp-reconstructor": "ApplyModels.disp_reconstructor_path",
         ("o", "output"): "ApplyModels.output_path",
         "chunk-size": "ApplyModels.chunk_size",
     }
@@ -107,6 +121,7 @@ class ApplyModels(Tool):
         TableLoader,
         EnergyRegressor,
         ParticleClassifier,
+        DispReconstructor,
         StereoCombiner,
     ]
 
@@ -129,6 +144,7 @@ class ApplyModels(Tool):
             load_instrument=True,
             load_dl1_images=False,
             load_simulated=False,
+            load_observation_info=True,
         )
 
         self._reconstructors = []
@@ -147,6 +163,13 @@ class ApplyModels(Tool):
                     parent=self,
                 )
             )
+        if self.disp_reconstructor_path is not None:
+            self._reconstructors.append(
+                DispReconstructor.read(
+                    self.disp_reconstructor_path,
+                    parent=self,
+                )
+            )
 
     def start(self):
         """Apply models to input tables"""
@@ -177,7 +200,6 @@ class ApplyModels(Tool):
 
             for tel_id, table in chunk.items():
                 tel = self.loader.subarray.tel[tel_id]
-
                 if tel not in reconstructor._models:
                     self.log.warning(
                         "No model in %s for telescope type %s, skipping tel %d",
@@ -195,18 +217,47 @@ class ApplyModels(Tool):
                     [c for c in table.colnames if c.startswith(prefix)]
                 )
 
-                predictions = reconstructor.predict_table(tel, table)
-                table = hstack(
-                    [table, predictions],
-                    join_type="exact",
-                    metadata_conflicts="ignore",
-                )
-                write_table(
-                    table[["obs_id", "event_id", "tel_id"] + predictions.colnames],
-                    self.output_path,
-                    f"/dl2/event/telescope/{property}/{prefix}/tel_{tel_id:03d}",
-                    append=True,
-                )
+                if isinstance(reconstructor, DispReconstructor):
+                    disp_predictions, altaz_predictions = reconstructor.predict_table(
+                        tel, table
+                    )
+                    table = hstack(
+                        [table, altaz_predictions, disp_predictions],
+                        join_type="exact",
+                        metadata_conflicts="ignore",
+                    )
+                    # tables should follow the container structure
+                    write_table(
+                        table[
+                            ["obs_id", "event_id", "tel_id"]
+                            + altaz_predictions.colnames
+                        ],
+                        self.output_path,
+                        f"/dl2/event/telescope/geometry/{prefix}/tel_{tel_id:03d}",
+                        append=True,
+                    )
+                    write_table(
+                        table[
+                            ["obs_id", "event_id", "tel_id"] + disp_predictions.colnames
+                        ],
+                        self.output_path,
+                        f"/dl2/event/telescope/disp/{prefix}/tel_{tel_id:03d}",
+                        append=True,
+                    )
+                else:
+                    predictions = reconstructor.predict_table(tel, table)
+                    table = hstack(
+                        [table, predictions],
+                        join_type="exact",
+                        metadata_conflicts="ignore",
+                    )
+                    write_table(
+                        table[["obs_id", "event_id", "tel_id"] + predictions.colnames],
+                        self.output_path,
+                        f"/dl2/event/telescope/{property}/{prefix}/tel_{tel_id:03d}",
+                        append=True,
+                    )
+
                 tel_tables.append(table)
 
             if len(tel_tables) == 0:
diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py
index b7ebbee3..5ceaa584 100644
--- a/ctapipe/tools/process.py
+++ b/ctapipe/tools/process.py
@@ -76,6 +76,7 @@ class ProcessorTool(Tool):
         ("m", "max-events"): "EventSource.max_events",
         "energy-regressor": "ShowerProcessor.EnergyRegressor.load_path",
         "particle-classifier": "ShowerProcessor.ParticleClassifier.load_path",
+        "disp-reconstructor": "ShowerProcessor.DispReconstructor.load_path",
         "image-cleaner-type": "ImageProcessor.image_cleaner_type",
     }
 
@@ -186,6 +187,7 @@ class ProcessorTool(Tool):
         reco_aliases = {
             "--energy-regressor": "EnergyRegressor",
             "--particle-classifier": "ParticleClassifier",
+            "--disp-reconstructor": "DispReconstructor",
         }
         for alias, name in reco_aliases.items():
             has_alias = any(arg.startswith(alias) for arg in self.argv)
diff --git a/ctapipe/tools/train_disp_reconstructor.py b/ctapipe/tools/train_disp_reconstructor.py
new file mode 100644
index 00000000..10ffa7c6
--- /dev/null
+++ b/ctapipe/tools/train_disp_reconstructor.py
@@ -0,0 +1,211 @@
+import astropy.units as u
+import numpy as np
+
+from ctapipe.core import Tool
+from ctapipe.core.traits import Bool, Int, IntTelescopeParameter, Path, TraitError, flag
+from ctapipe.io import TableLoader
+from ctapipe.reco import CrossValidator, DispReconstructor
+from ctapipe.reco.preprocessing import check_valid_rows, horizontal_to_telescope
+
+
+class TrainDispReconstructor(Tool):
+    """
+    Tool to train a `~ctapipe.reco.DispReconstructor` on dl1b/dl2 data.
+
+    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.
+    """
+
+    name = "ctapipe-train-disp-reconstructor"
+    description = __doc__
+
+    examples = """
+    ctapipe-train-disp-reconstructor \\
+        --config train_disp_reconstructor.yaml \\
+        --input gamma.dl2.h5 \\
+        --output disp_models.pkl
+    """
+
+    output_path = Path(
+        default_value=None,
+        allow_none=False,
+        directory_ok=False,
+        help=(
+            "Ouput path for the trained reconstructor."
+            " At the moment, pickle is the only supported format.",
+        ),
+    ).tag(config=True)
+
+    n_events = IntTelescopeParameter(
+        default_value=None,
+        allow_none=True,
+        help=(
+            "Number of events for training the models."
+            " If not given, all available events will be used."
+        ),
+    ).tag(config=True)
+
+    overwrite = Bool(help="overwrite existing output files").tag(config=True)
+
+    random_seed = Int(
+        default_value=0, help="Random seed for sampling and cross validation"
+    ).tag(config=True)
+
+    project_disp = Bool(
+        default_value=False,
+        help=(
+            "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"
+            " and the true source position."
+        ),
+    ).tag(config=True)
+
+    flags = {
+        **flag(
+            "overwrite",
+            "TrainDispReconstructor.overwrite",
+            "Overwrite output existing output files",
+            "Don't overwrite existing output files",
+        ),
+    }
+
+    aliases = {
+        ("i", "input"): "TableLoader.input_url",
+        ("o", "output"): "TrainDispReconstructor.output_path",
+        "n-events": "TrainDispReconstructor.n_events",
+        "cv-output": "CrossValidator.output_path",
+    }
+
+    classes = [TableLoader, DispReconstructor, CrossValidator]
+
+    def setup(self):
+        """
+        Initialize components from config.
+        """
+        self.loader = TableLoader(
+            parent=self,
+            load_dl1_images=False,
+            load_dl1_parameters=True,
+            load_dl2=True,
+            load_simulated=True,
+            load_instrument=True,
+            load_observation_info=True,
+        )
+        self.n_events.attach_subarray(self.loader.subarray)
+
+        self.models = DispReconstructor(self.loader.subarray, parent=self)
+        self.cross_validate = CrossValidator(parent=self, model_component=self.models)
+        self.rng = np.random.default_rng(self.random_seed)
+
+        if self.output_path.suffix != ".pkl":
+            self.log.warning(
+                "Expected .pkl extension for output_path, got %s",
+                self.output_path.suffix,
+            )
+
+        if self.output_path.exists() and not self.overwrite:
+            raise TraitError(
+                f"output_path '{self.output_path}' exists and overwrite=False"
+            )
+
+    def start(self):
+        """
+        Train models per telescope type using a cross-validation.
+        """
+        types = self.loader.subarray.telescope_types
+        self.log.info("Inputfile: %s", self.loader.input_url)
+
+        self.log.info("Training models for %d types", len(types))
+        for tel_type in types:
+            self.log.info("Loading events for %s", tel_type)
+            table = self._read_table(tel_type)
+
+            self.log.info("Train models on %s events", len(table))
+            self.cross_validate(tel_type, table)
+
+            self.log.info("Performing final fit for %s", tel_type)
+            self.models.fit(tel_type, table)
+            self.log.info("done")
+
+    def _read_table(self, telescope_type):
+        table = self.loader.read_telescope_events([telescope_type])
+        self.log.info("Events read from input: %d", len(table))
+
+        mask = self.models.quality_query.get_table_mask(table)
+        table = table[mask]
+        self.log.info("Events after applying quality query: %d", len(table))
+
+        table = self.models.feature_generator(table)
+
+        table[self.models.target] = self._get_true_disp(table)
+
+        # Add true energy for energy-dependent performance plots
+        columns = self.models.features + [self.models.target, "true_energy"]
+        table = table[columns]
+
+        valid = check_valid_rows(table)
+        if np.any(~valid):
+            self.log.warning("Dropping non-predicable events.")
+            table = table[valid]
+
+        n_events = self.n_events.tel[telescope_type]
+        if n_events is not None:
+            if n_events > len(table):
+                self.log.warning(
+                    "Number of events in table (%d) is less than requested number of events %d",
+                    len(table),
+                    n_events,
+                )
+            else:
+                self.log.info("Sampling %d events", n_events)
+                idx = self.rng.choice(len(table), n_events, replace=False)
+                idx.sort()
+                table = table[idx]
+
+        return table
+
+    def _get_true_disp(self, table):
+        fov_lon, fov_lat = horizontal_to_telescope(
+            alt=table["true_alt"],
+            az=table["true_az"],
+            pointing_alt=table["subarray_pointing_lat"],
+            pointing_az=table["subarray_pointing_lon"],
+        )
+
+        # numpy's trigonometric functions need radians
+        psi = table["hillas_psi"].quantity.to_value(u.rad)
+        cog_lon = table["hillas_fov_lon"].quantity
+        cog_lat = table["hillas_fov_lat"].quantity
+
+        delta_lon = fov_lon - cog_lon
+        delta_lat = fov_lat - cog_lat
+
+        true_disp = np.cos(psi) * delta_lon + np.sin(psi) * delta_lat
+        true_sign = np.sign(true_disp)
+
+        if self.project_disp:
+            true_norm = np.abs(true_disp)
+        else:
+            true_norm = np.sqrt((fov_lon - cog_lon) ** 2 + (fov_lat - cog_lat) ** 2)
+
+        return true_norm * true_sign
+
+    def finish(self):
+        """
+        Write-out trained models and cross-validation results.
+        """
+        self.log.info("Writing output")
+        self.models.write(self.output_path, overwrite=self.overwrite)
+        if self.cross_validate.output_path:
+            self.cross_validate.write(overwrite=self.overwrite)
+        self.loader.close()
+
+
+def main():
+    TrainDispReconstructor().run()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/ctapipe/tools/train_energy_regressor.py b/ctapipe/tools/train_energy_regressor.py
index d8159ac3..a4ca36a5 100644
--- a/ctapipe/tools/train_energy_regressor.py
+++ b/ctapipe/tools/train_energy_regressor.py
@@ -30,7 +30,7 @@ class TrainEnergyRegressor(Tool):
     ctapipe-train-energy-regressor \\
         --config train_energy_regressor.yaml \\
         --input gamma.dl2.h5 \\
-        --ouput energy_regressor.pkl
+        --output energy_regressor.pkl
     """
 
     output_path = Path(
diff --git a/ctapipe/tools/train_particle_classifier.py b/ctapipe/tools/train_particle_classifier.py
index 509cc31d..46f73dd8 100644
--- a/ctapipe/tools/train_particle_classifier.py
+++ b/ctapipe/tools/train_particle_classifier.py
@@ -87,7 +87,7 @@ class TrainParticleClassifier(Tool):
     flags = {
         **flag(
             "overwrite",
-            "TrainEnergyRegressor.overwrite",
+            "TrainParticleClassifier.overwrite",
             "Overwrite output existing output files",
             "Don't overwrite existing output files",
         ),
@@ -189,7 +189,8 @@ class TrainParticleClassifier(Tool):
 
         table = self.classifier.feature_generator(table)
 
-        columns = self.classifier.features + [self.classifier.target]
+        # Add true energy for energy-dependent performance plots
+        columns = self.classifier.features + [self.classifier.target, "true_energy"]
         table = table[columns]
 
         valid = check_valid_rows(table)
diff --git a/docs/changes/2138.feature.rst b/docs/changes/2138.feature.rst
new file mode 100644
index 00000000..557ed3d7
--- /dev/null
+++ b/docs/changes/2138.feature.rst
@@ -0,0 +1,5 @@
+This adds a disp reconstruction including ``ctapipe-train-disp-reconstructor`` for training
+two machine learning models and application of these models in ``ctapipe-apply-models``
+and ``ctapipe-process``.
+``DispContainer`` is added for storing the disp prediction.
+Stereo predictions are done by taking a weighted average of the single telescope predictions.
diff --git a/setup.cfg b/setup.cfg
index 3bcec0d1..bb217f6a 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -96,6 +96,7 @@ console_scripts =
     ctapipe-quickstart = ctapipe.tools.quickstart:main
     ctapipe-train-energy-regressor = ctapipe.tools.train_energy_regressor:main
     ctapipe-train-particle-classifier = ctapipe.tools.train_particle_classifier:main
+    ctapipe-train-disp-reconstructor = ctapipe.tools.train_disp_reconstructor:main
     ctapipe-apply-models = ctapipe.tools.apply_models:main
 
 


Test patch:
diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index 18396dba..dec57bff 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -541,9 +541,7 @@ def proton_train_clf(model_tmp_path, energy_regressor_path):
 
 
 @pytest.fixture(scope="session")
-def particle_classifier_path(
-    model_tmp_path, energy_regressor_path, gamma_train_clf, proton_train_clf
-):
+def particle_classifier_path(model_tmp_path, gamma_train_clf, proton_train_clf):
     from ctapipe.tools.train_particle_classifier import TrainParticleClassifier
 
     out_file = model_tmp_path / "particle_classifier.pkl"
@@ -565,3 +563,27 @@ def particle_classifier_path(
         )
         assert ret == 0
         return out_file
+
+
+@pytest.fixture(scope="session")
+def disp_reconstructor_path(model_tmp_path, gamma_train_clf):
+    from ctapipe.tools.train_disp_reconstructor import TrainDispReconstructor
+
+    out_file = model_tmp_path / "disp_reconstructor.pkl"
+    with FileLock(out_file.with_suffix(out_file.suffix + ".lock")):
+        if out_file.is_file():
+            return out_file
+
+        config = resource_file("train_disp_reconstructor.yaml")
+
+        ret = run_tool(
+            TrainDispReconstructor(),
+            argv=[
+                f"--input={gamma_train_clf}",
+                f"--output={out_file}",
+                f"--config={config}",
+                "--log-level=INFO",
+            ],
+        )
+        assert ret == 0
+        return out_file
diff --git a/ctapipe/reco/tests/test_stereo_combination.py b/ctapipe/reco/tests/test_stereo_combination.py
index 2cd8aec4..2baa151f 100644
--- a/ctapipe/reco/tests/test_stereo_combination.py
+++ b/ctapipe/reco/tests/test_stereo_combination.py
@@ -2,7 +2,7 @@ import astropy.units as u
 import numpy as np
 import pytest
 from astropy.table import Table
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_allclose, assert_array_equal
 
 from ctapipe.containers import (
     ArrayEventContainer,
@@ -11,6 +11,7 @@ from ctapipe.containers import (
     ParticleClassificationContainer,
     ReconstructedContainer,
     ReconstructedEnergyContainer,
+    ReconstructedGeometryContainer,
 )
 from ctapipe.reco.stereo_combination import StereoMeanCombiner
 
@@ -47,6 +48,16 @@ def mono_table():
                 True,
                 True,
             ],
+            "disp_tel_alt": [58.5, 58, 62.5, 72, 74.5, 81] * u.deg,
+            "disp_tel_az": [12.5, 15, 13, 21, 20, 14.5] * u.deg,
+            "disp_tel_is_valid": [
+                True,
+                False,
+                True,
+                True,
+                True,
+                True,
+            ],
         }
     )
 
@@ -105,6 +116,39 @@ def test_predict_mean_classification(mono_table):
     assert_array_equal(tel_ids[2], [1])
 
 
+def test_predict_mean_disp(mono_table):
+    combine = StereoMeanCombiner(
+        prefix="disp",
+        property="geometry",
+    )
+    stereo = combine.predict_table(mono_table)
+
+    for name, field in ReconstructedGeometryContainer.fields.items():
+        colname = f"disp_{name}"
+        assert colname in stereo.colnames
+        assert stereo[colname].description == field.description
+
+    assert "obs_id" in stereo.colnames
+    assert "event_id" in stereo.colnames
+
+    assert_array_equal(stereo["obs_id"], np.array([1, 1, 2]))
+    assert_array_equal(stereo["event_id"], np.array([1, 2, 1]))
+    assert_allclose(
+        stereo["disp_alt"].quantity,
+        [60.5002328, 73.2505989, 81] * u.deg,
+        atol=1e-7,
+    )
+    assert_allclose(
+        stereo["disp_az"].quantity,
+        [12.7345693, 20.5362510, 14.5] * u.deg,
+        atol=1e-7,
+    )
+    tel_ids = stereo["disp_telescopes"]
+    assert_array_equal(tel_ids[0], [1, 3])
+    assert_array_equal(tel_ids[1], [5, 7])
+    assert_array_equal(tel_ids[2], [1])
+
+
 @pytest.mark.parametrize("weights", ["konrad", "intensity", "none"])
 def test_mean_prediction_single_event(weights):
     event = ArrayEventContainer()
@@ -125,6 +169,11 @@ def test_mean_prediction_single_event(weights):
         classification={
             "dummy": ParticleClassificationContainer(prediction=1.0, is_valid=True)
         },
+        geometry={
+            "dummy": ReconstructedGeometryContainer(
+                alt=60 * u.deg, az=15 * u.deg, is_valid=True
+            )
+        },
     )
     event.dl2.tel[125] = ReconstructedContainer(
         energy={
@@ -133,6 +182,11 @@ def test_mean_prediction_single_event(weights):
         classification={
             "dummy": ParticleClassificationContainer(prediction=0.0, is_valid=True)
         },
+        geometry={
+            "dummy": ReconstructedGeometryContainer(
+                alt=50 * u.deg, az=30 * u.deg, is_valid=True
+            )
+        },
     )
     event.dl2.tel[130] = ReconstructedContainer(
         energy={
@@ -141,6 +195,11 @@ def test_mean_prediction_single_event(weights):
         classification={
             "dummy": ParticleClassificationContainer(prediction=0.8, is_valid=True)
         },
+        geometry={
+            "dummy": ReconstructedGeometryContainer(
+                alt=45 * u.deg, az=280 * u.deg, is_valid=True
+            )
+        },
     )
 
     combine_energy = StereoMeanCombiner(
@@ -153,10 +212,20 @@ def test_mean_prediction_single_event(weights):
         property="classification",
         weights=weights,
     )
+    combine_geometry = StereoMeanCombiner(
+        prefix="dummy",
+        property="geometry",
+        weights=weights,
+    )
     combine_energy(event)
     combine_classification(event)
+    combine_geometry(event)
     if weights == "none":
         assert u.isclose(event.dl2.stereo.energy["dummy"].energy, (70 / 3) * u.GeV)
+        assert u.isclose(event.dl2.stereo.geometry["dummy"].alt, 63.0738383 * u.deg)
+        assert u.isclose(event.dl2.stereo.geometry["dummy"].az, 348.0716693 * u.deg)
     elif weights == "intensity":
         assert u.isclose(event.dl2.stereo.energy["dummy"].energy, 30 * u.GeV)
+        assert u.isclose(event.dl2.stereo.geometry["dummy"].alt, 60.9748605 * u.deg)
+        assert u.isclose(event.dl2.stereo.geometry["dummy"].az, 316.0365515 * u.deg)
     assert event.dl2.stereo.classification["dummy"].prediction == 0.6
diff --git a/ctapipe/tools/tests/test_apply_models.py b/ctapipe/tools/tests/test_apply_models.py
index 2735f91a..a905efa1 100644
--- a/ctapipe/tools/tests/test_apply_models.py
+++ b/ctapipe/tools/tests/test_apply_models.py
@@ -4,6 +4,7 @@ from ctapipe.containers import (
     EventIndexContainer,
     ParticleClassificationContainer,
     ReconstructedEnergyContainer,
+    ReconstructedGeometryContainer,
 )
 from ctapipe.core import run_tool
 from ctapipe.io import TableLoader, read_table
@@ -68,16 +69,17 @@ def test_apply_energy_regressor(
     check_equal_array_event_order(trigger, energy)
 
 
-def test_apply_both(
+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
 
     input_path = dl2_shower_geometry_file_lapalma
-    output_path = tmp_path / "particle-and-energy.dl2.h5"
+    output_path = tmp_path / "particle-and-energy-and-disp.dl2.h5"
 
     ret = run_tool(
         ApplyModels(),
@@ -86,37 +88,64 @@ def test_apply_both(
             f"--output={output_path}",
             f"--energy-regressor={energy_regressor_path}",
             f"--particle-classifier={particle_classifier_path}",
+            f"--disp-reconstructor={disp_reconstructor_path}",
             "--StereoMeanCombiner.weights=konrad",
             "--chunk-size=5",  # small chunksize so we test multiple chunks for the test file
         ],
+        raises=True,
     )
     assert ret == 0
 
-    prefix = "ExtraTreesClassifier"
-    table = read_table(output_path, f"/dl2/event/subarray/classification/{prefix}")
+    prefix_clf = "ExtraTreesClassifier"
+    table = read_table(output_path, f"/dl2/event/subarray/classification/{prefix_clf}")
     for col in "obs_id", "event_id":
         assert table[col].description == EventIndexContainer.fields[col].description
 
     for name, field in ParticleClassificationContainer.fields.items():
-        colname = f"ExtraTreesClassifier_{name}"
+        colname = f"{prefix_clf}_{name}"
+        assert colname in table.colnames
+        assert table[colname].description == field.description
+
+    prefix_disp = "disp"
+    table = read_table(output_path, f"/dl2/event/subarray/geometry/{prefix_disp}")
+    for col in "obs_id", "event_id":
+        assert table[col].description == EventIndexContainer.fields[col].description
+
+    for name, field in ReconstructedGeometryContainer.fields.items():
+        colname = f"{prefix_disp}_{name}"
         assert colname in table.colnames
         assert table[colname].description == field.description
 
     loader = TableLoader(output_path, load_dl2=True)
     events = loader.read_subarray_events()
-    assert f"{prefix}_prediction" in events.colnames
-    assert f"{prefix}_telescopes" in events.colnames
-    assert f"{prefix}_is_valid" in events.colnames
-    assert f"{prefix}_goodness_of_fit" in events.colnames
+    assert f"{prefix_clf}_prediction" in events.colnames
+    assert f"{prefix_clf}_telescopes" in events.colnames
+    assert f"{prefix_clf}_is_valid" in events.colnames
+    assert f"{prefix_clf}_goodness_of_fit" in events.colnames
+    assert f"{prefix_disp}_alt" in events.colnames
+    assert f"{prefix_disp}_az" in events.colnames
+    assert f"{prefix_disp}_ang_distance_uncert" in events.colnames
+    assert f"{prefix_disp}_is_valid" in events.colnames
+    assert f"{prefix_disp}_goodness_of_fit" in events.colnames
 
     events = loader.read_telescope_events()
-    assert f"{prefix}_prediction" in events.colnames
-    assert f"{prefix}_telescopes" in events.colnames
-    assert f"{prefix}_is_valid" in events.colnames
-    assert f"{prefix}_goodness_of_fit" in events.colnames
-
-    assert f"{prefix}_tel_prediction" in events.colnames
-    assert f"{prefix}_tel_is_valid" in events.colnames
+    assert f"{prefix_clf}_prediction" in events.colnames
+    assert f"{prefix_clf}_telescopes" in events.colnames
+    assert f"{prefix_clf}_is_valid" in events.colnames
+    assert f"{prefix_clf}_goodness_of_fit" in events.colnames
+    assert f"{prefix_disp}_alt" in events.colnames
+    assert f"{prefix_disp}_az" in events.colnames
+    assert f"{prefix_disp}_ang_distance_uncert" in events.colnames
+    assert f"{prefix_disp}_is_valid" in events.colnames
+    assert f"{prefix_disp}_goodness_of_fit" in events.colnames
+
+    assert f"{prefix_clf}_tel_prediction" in events.colnames
+    assert f"{prefix_clf}_tel_is_valid" in events.colnames
+    assert f"{prefix_disp}_tel_alt" in events.colnames
+    assert f"{prefix_disp}_tel_az" in events.colnames
+    assert f"{prefix_disp}_tel_is_valid" in events.colnames
+    assert f"{prefix_disp}_parameter_norm" in events.colnames
+    assert f"{prefix_disp}_parameter_is_valid" in events.colnames
 
     assert "ExtraTreesRegressor_energy" in events.colnames
 
@@ -124,6 +153,8 @@ def test_apply_both(
 
     trigger = read_table(output_path, "/dl1/event/subarray/trigger")
     particle_clf = read_table(
-        output_path, "/dl2/event/subarray/classification/ExtraTreesClassifier"
+        output_path, f"/dl2/event/subarray/classification/{prefix_clf}"
     )
     check_equal_array_event_order(trigger, particle_clf)
+    disp_reco = read_table(output_path, f"/dl2/event/subarray/geometry/{prefix_disp}")
+    check_equal_array_event_order(trigger, disp_reco)
diff --git a/ctapipe/tools/tests/test_process_ml.py b/ctapipe/tools/tests/test_process_ml.py
index 3b011fd8..6e999bc7 100644
--- a/ctapipe/tools/tests/test_process_ml.py
+++ b/ctapipe/tools/tests/test_process_ml.py
@@ -115,3 +115,66 @@ def test_process_apply_classification(
     assert "ExtraTreesClassifier_prediction" in events.colnames
     np.testing.assert_array_equal(events["obs_id"], trigger["obs_id"])
     np.testing.assert_array_equal(events["event_id"], trigger["event_id"])
+
+
+def test_process_apply_disp(
+    tmp_path,
+    energy_regressor_path,
+    disp_reconstructor_path,
+    prod5_gamma_lapalma_simtel_path,
+):
+    from ctapipe.tools.process import ProcessorTool
+
+    output = tmp_path / "gamma_prod5.dl2_disp.h5"
+
+    config_path = tmp_path / "config.json"
+
+    input_url = prod5_gamma_lapalma_simtel_path
+
+    allowed_tels = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 19, 35]
+    config = {
+        "ProcessorTool": {
+            "EventSource": {
+                "allowed_tels": allowed_tels,
+            },
+            "ShowerProcessor": {
+                "reconstructor_types": [
+                    "HillasReconstructor",
+                    "EnergyRegressor",
+                    "DispReconstructor",
+                ]
+            },
+        }
+    }
+
+    with config_path.open("w") as f:
+        json.dump(config, f)
+
+    argv = [
+        f"--input={input_url}",
+        f"--output={output}",
+        "--write-images",
+        "--write-showers",
+        f"--energy-regressor={energy_regressor_path}",
+        f"--disp-reconstructor={disp_reconstructor_path}",
+        f"--config={config_path}",
+    ]
+    assert run_tool(ProcessorTool(), argv=argv, cwd=tmp_path, raises=True) == 0
+
+    tel_events = read_table(output, "/dl2/event/telescope/disp/disp/tel_004")
+    assert "disp_parameter_norm" in tel_events.colnames
+    assert "disp_parameter_is_valid" in tel_events.colnames
+
+    tel_events = read_table(output, "/dl2/event/telescope/geometry/disp/tel_004")
+    assert "disp_tel_alt" in tel_events.colnames
+    assert "disp_tel_az" in tel_events.colnames
+    assert "disp_tel_is_valid" in tel_events.colnames
+
+    events = read_table(output, "/dl2/event/subarray/geometry/disp")
+    trigger = read_table(output, "/dl1/event/subarray/trigger")
+    assert len(events) == len(trigger)
+    assert "disp_alt" in events.colnames
+    assert "disp_az" in events.colnames
+    assert "disp_is_valid" in events.colnames
+    np.testing.assert_array_equal(events["obs_id"], trigger["obs_id"])
+    np.testing.assert_array_equal(events["event_id"], trigger["event_id"])
diff --git a/ctapipe/tools/tests/test_train.py b/ctapipe/tools/tests/test_train.py
index 3d2e5f1e..ac1b28b9 100644
--- a/ctapipe/tools/tests/test_train.py
+++ b/ctapipe/tools/tests/test_train.py
@@ -16,6 +16,12 @@ def test_train_particle_classifier(particle_classifier_path):
     ParticleClassifier.read(particle_classifier_path)
 
 
+def test_train_disp_reconstructor(disp_reconstructor_path):
+    from ctapipe.reco import DispReconstructor
+
+    DispReconstructor.read(disp_reconstructor_path)
+
+
 def test_too_few_events(tmp_path, dl2_shower_geometry_file):
     from ctapipe.tools.train_energy_regressor import TrainEnergyRegressor
 
@@ -37,6 +43,7 @@ def test_too_few_events(tmp_path, dl2_shower_geometry_file):
 
 
 def test_cross_validation_results(tmp_path, gamma_train_clf, proton_train_clf):
+    from ctapipe.tools.train_disp_reconstructor import TrainDispReconstructor
     from ctapipe.tools.train_energy_regressor import TrainEnergyRegressor
     from ctapipe.tools.train_particle_classifier import TrainParticleClassifier
 
@@ -76,3 +83,21 @@ def test_cross_validation_results(tmp_path, gamma_train_clf, proton_train_clf):
     )
     assert ret == 0
     assert classifier_cv_out_file.exists()
+
+    tool = TrainDispReconstructor()
+    config = resource_file("train_disp_reconstructor.yaml")
+    out_file = tmp_path / "disp_reconstructor_.pkl"
+    disp_cv_out_file = tmp_path / "disp_cv_results.h5"
+
+    ret = run_tool(
+        tool,
+        argv=[
+            f"--input={gamma_train_clf}",
+            f"--output={out_file}",
+            f"--config={config}",
+            f"--cv-output={disp_cv_out_file}",
+            "--log-level=INFO",
+        ],
+    )
+    assert ret == 0
+    assert disp_cv_out_file.exists()