Repository: cta-observatory/ctapipe (python)

Issue:
Wrong sign in disp reconstruction
Hi there!

I'm attempting to create IRFs of a single SST using the `ctapipe 0.24.0` CLI.

After disp reconstruction, many point-source gamma-rays are well-reconstructed, but many are not, clustered in a ring around the true source direction. The resulting angular resolution is thus very poor.

I don't know how best to diagnose this. My first thoughts are that it maybe has something to do with sign reconstruction, or perhaps coordinate transformations in a dual-mirror telescope?

![Image](https://github.com/user-attachments/assets/ba016b8e-f5af-426c-baba-4709a8ff9f09)

![Image](https://github.com/user-attachments/assets/407de11d-5e96-4e5c-96fe-b4953f60bf62)

Here is my `train_disp_reconstructor.yaml`

```
TrainDispReconstructor:
  random_seed: 0
  n_events:

  CrossValidator:
    n_cross_validations: 5
    rng_seed: 1337

  DispReconstructor:
    n_jobs: -1

    norm_cls: RandomForestRegressor
    log_target: True
    sign_cls: RandomForestClassifier

    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:
        - ["area", "hillas_width * hillas_length"]

    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
      - peak_time_std
      - peak_time_skewness
      - peak_time_kurtosis
      - intensity_max
      - intensity_mean
      - intensity_std
      - intensity_skewness
      - intensity_kurtosis
      - area
```

Pull request description:
Fix wrong proba column being used for disp sign
fixes #2737

Patch:
diff --git a/docs/changes/2738.bugfix.rst b/docs/changes/2738.bugfix.rst
new file mode 100644
index 00000000..565aeebf
--- /dev/null
+++ b/docs/changes/2738.bugfix.rst
@@ -0,0 +1,1 @@
+Fix inverted sign of ``DispReconstructor`` prediction.
diff --git a/src/ctapipe/reco/sklearn.py b/src/ctapipe/reco/sklearn.py
index 835339f3..acc77e96 100644
--- a/src/ctapipe/reco/sklearn.py
+++ b/src/ctapipe/reco/sklearn.py
@@ -703,7 +703,7 @@ class DispReconstructor(Reconstructor):
             else:
                 prediction[valid] = valid_norms
 
-            sign_proba = self._models[key][1].predict_proba(X)[:, 0]
+            sign_proba = self._models[key][1].predict_proba(X)[:, 1]
             # proba is [0 and 1] where 0 => very certain -1, 1 => very certain 1
             # and 0.5 means random guessing either. So we transform to a score
             # where 0 means "guessing" and 1 means "very certain"


Test patch:
diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py
index f4555044..57f0a618 100644
--- a/src/ctapipe/conftest.py
+++ b/src/ctapipe/conftest.py
@@ -642,6 +642,7 @@ 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"
+    cv_out_file = model_tmp_path / "cv_disp_reconstructor.h5"
     with FileLock(out_file.with_suffix(out_file.suffix + ".lock")):
         if out_file.is_file():
             return out_file
@@ -653,12 +654,13 @@ def disp_reconstructor_path(model_tmp_path, gamma_train_clf):
             argv=[
                 f"--input={gamma_train_clf}",
                 f"--output={out_file}",
+                f"--cv-output={cv_out_file}",
                 f"--config={config}",
                 "--log-level=INFO",
             ],
         )
         assert ret == 0
-        return out_file
+        return out_file, cv_out_file
 
 
 @pytest.fixture(scope="session")
diff --git a/src/ctapipe/tools/tests/test_apply_models.py b/src/ctapipe/tools/tests/test_apply_models.py
index 6958a3b5..40ace9a5 100644
--- a/src/ctapipe/tools/tests/test_apply_models.py
+++ b/src/ctapipe/tools/tests/test_apply_models.py
@@ -110,6 +110,8 @@ def test_apply_all(
 ):
     from ctapipe.tools.apply_models import ApplyModels
 
+    disp_reconstructor_path, _ = disp_reconstructor_path
+
     input_path = get_dataset_path("gamma_diffuse_dl2_train_small.dl2.h5")
     output_path = tmp_path / "particle-and-energy-and-disp.dl2.h5"
 
diff --git a/src/ctapipe/tools/tests/test_process_ml.py b/src/ctapipe/tools/tests/test_process_ml.py
index 6a25f503..f9445696 100644
--- a/src/ctapipe/tools/tests/test_process_ml.py
+++ b/src/ctapipe/tools/tests/test_process_ml.py
@@ -122,6 +122,8 @@ def test_process_apply_disp(
 ):
     from ctapipe.tools.process import ProcessorTool
 
+    disp_reconstructor_path, _ = disp_reconstructor_path
+
     output = tmp_path / "gamma_prod5.dl2_disp.h5"
 
     config_path = tmp_path / "config.json"
diff --git a/src/ctapipe/tools/tests/test_train.py b/src/ctapipe/tools/tests/test_train.py
index acdb6092..a5485820 100644
--- a/src/ctapipe/tools/tests/test_train.py
+++ b/src/ctapipe/tools/tests/test_train.py
@@ -3,6 +3,7 @@ import pytest
 
 from ctapipe.core import ToolConfigurationError, run_tool
 from ctapipe.exceptions import TooFewEvents
+from ctapipe.io import read_table
 from ctapipe.utils.datasets import resource_file
 
 
@@ -21,7 +22,15 @@ def test_train_particle_classifier(particle_classifier_path):
 def test_train_disp_reconstructor(disp_reconstructor_path):
     from ctapipe.reco import DispReconstructor
 
-    DispReconstructor.read(disp_reconstructor_path)
+    model_path, cv_path = disp_reconstructor_path
+
+    DispReconstructor.read(model_path)
+
+    cv_table = read_table(cv_path, "/cv_predictions/LST_LST_LSTCam")
+    disp = cv_table["disp_parameter"]
+    true_disp = cv_table["truth"]
+    accuracy = np.count_nonzero(np.sign(disp) == np.sign(true_disp)) / len(disp)
+    assert accuracy > 0.75
 
 
 def test_too_few_events(tmp_path, dl2_shower_geometry_file):