Remove option of only calculating gh cuts with cut-opt tool, but keep…

… possibility to save an OptimizationResult with only gh cuts for now.

src/ctapipe/irf/optimize.py

@@ -188,7 +188,6 @@ def optimize_cuts(
         background: QTable,
         precuts: EventPreProcessor,
         clf_prefix: str,
-        point_like: bool,
     ) -> OptimizationResult:
         """
         Optimize G/H (and optionally theta) cuts
@@ -206,9 +205,6 @@ def optimize_cuts(
         clf_prefix: str
             Prefix of the output from the G/H classifier for which the
             cut will be optimized
-        point_like: bool
-            Whether a theta cut should be calculated (True) or only a
-            G/H cut (False)
         """
 
 
@@ -331,7 +327,6 @@ def optimize_cuts(
         background: QTable,
         precuts: EventPreProcessor,
         clf_prefix: str,
-        point_like: bool,
     ) -> OptimizationResult:
         reco_energy_bins = make_bins_per_decade(
             self.reco_energy_min.to(u.TeV),
@@ -343,18 +338,17 @@ def optimize_cuts(
             signal["reco_energy"],
             reco_energy_bins,
         )
-        if point_like:
-            gh_mask = evaluate_binned_cut(
-                signal["gh_score"],
-                signal["reco_energy"],
-                gh_cuts,
-                op=operator.ge,
-            )
-            theta_cuts = self.theta.calculate_theta_cut(
-                signal["theta"][gh_mask],
-                signal["reco_energy"][gh_mask],
-                reco_energy_bins,
-            )
+        gh_mask = evaluate_binned_cut(
+            signal["gh_score"],
+            signal["reco_energy"],
+            gh_cuts,
+            op=operator.ge,
+        )
+        theta_cuts = self.theta.calculate_theta_cut(
+            signal["theta"][gh_mask],
+            signal["reco_energy"][gh_mask],
+            reco_energy_bins,
+        )
 
         result = OptimizationResult(
             precuts=precuts,
@@ -365,7 +359,7 @@ def optimize_cuts(
             # A single set of cuts is calculated for the whole fov atm
             valid_offset_min=0 * u.deg,
             valid_offset_max=np.inf * u.deg,
-            theta_cuts=theta_cuts if point_like else None,
+            theta_cuts=theta_cuts,
         )
         return result
 
@@ -406,49 +400,35 @@ def optimize_cuts(
         background: QTable,
         precuts: EventPreProcessor,
         clf_prefix: str,
-        point_like: bool,
     ) -> OptimizationResult:
         reco_energy_bins = make_bins_per_decade(
             self.reco_energy_min.to(u.TeV),
             self.reco_energy_max.to(u.TeV),
             self.reco_energy_n_bins_per_decade,
         )
 
-        if point_like:
-            initial_gh_cuts = calculate_percentile_cut(
-                signal["gh_score"],
-                signal["reco_energy"],
-                bins=reco_energy_bins,
-                fill_value=0.0,
-                percentile=100 * (1 - self.initial_gh_cut_efficency),
-                min_events=10,
-                smoothing=1,
-            )
-            initial_gh_mask = evaluate_binned_cut(
-                signal["gh_score"],
-                signal["reco_energy"],
-                initial_gh_cuts,
-                op=operator.gt,
-            )
+        initial_gh_cuts = calculate_percentile_cut(
+            signal["gh_score"],
+            signal["reco_energy"],
+            bins=reco_energy_bins,
+            fill_value=0.0,
+            percentile=100 * (1 - self.initial_gh_cut_efficency),
+            min_events=10,
+            smoothing=1,
+        )
+        initial_gh_mask = evaluate_binned_cut(
+            signal["gh_score"],
+            signal["reco_energy"],
+            initial_gh_cuts,
+            op=operator.gt,
+        )
 
-            theta_cuts = self.theta.calculate_theta_cut(
-                signal["theta"][initial_gh_mask],
-                signal["reco_energy"][initial_gh_mask],
-                reco_energy_bins,
-            )
-            self.log.info("Optimizing G/H separation cut for best sensitivity")
-        else:
-            # Create a dummy theta cut since `pyirf.cut_optimization.optimize_gh_cut`
-            # needs a theta cut atm.
-            theta_cuts = QTable()
-            theta_cuts["low"] = reco_energy_bins[:-1]
-            theta_cuts["center"] = 0.5 * (reco_energy_bins[:-1] + reco_energy_bins[1:])
-            theta_cuts["high"] = reco_energy_bins[1:]
-            theta_cuts["cut"] = self.max_bkg_fov_offset
-            self.log.info(
-                "Optimizing G/H separation cut for best sensitivity "
-                "with `max_bkg_fov_offset` as theta cut."
-            )
+        theta_cuts = self.theta.calculate_theta_cut(
+            signal["theta"][initial_gh_mask],
+            signal["reco_energy"][initial_gh_mask],
+            reco_energy_bins,
+        )
+        self.log.info("Optimizing G/H separation cut for best sensitivity")
 
         gh_cut_efficiencies = np.arange(
             self.gh_cut_efficiency_step,
@@ -469,18 +449,17 @@ def optimize_cuts(
         valid_energy = self._get_valid_energy_range(opt_sens)
 
         # Re-calculate theta cut with optimized g/h cut
-        if point_like:
-            signal["selected_gh"] = evaluate_binned_cut(
-                signal["gh_score"],
-                signal["reco_energy"],
-                gh_cuts,
-                operator.ge,
-            )
-            theta_cuts_opt = self.theta.calculate_theta_cut(
-                signal[signal["selected_gh"]]["theta"],
-                signal[signal["selected_gh"]]["reco_energy"],
-                reco_energy_bins,
-            )
+        signal["selected_gh"] = evaluate_binned_cut(
+            signal["gh_score"],
+            signal["reco_energy"],
+            gh_cuts,
+            operator.ge,
+        )
+        theta_cuts_opt = self.theta.calculate_theta_cut(
+            signal[signal["selected_gh"]]["theta"],
+            signal[signal["selected_gh"]]["reco_energy"],
+            reco_energy_bins,
+        )
 
         result = OptimizationResult(
             precuts=precuts,
@@ -491,7 +470,7 @@ def optimize_cuts(
             # A single set of cuts is calculated for the whole fov atm
             valid_offset_min=self.min_bkg_fov_offset,
             valid_offset_max=self.max_bkg_fov_offset,
-            theta_cuts=theta_cuts_opt if point_like else None,
+            theta_cuts=theta_cuts_opt,
         )
         return result
 

src/ctapipe/irf/tests/test_optimize.py

@@ -116,7 +116,6 @@ def test_cut_optimizer(
         background=proton_events,
         precuts=gamma_loader.epp,  # identical precuts for all particle types
         clf_prefix="ExtraTreesClassifier",
-        point_like=True,
     )
     assert isinstance(result, OptimizationResult)
     assert result.clf_prefix == "ExtraTreesClassifier"

src/ctapipe/tools/compute_irf.py

@@ -13,7 +13,7 @@
 import astropy.units as u
 import numpy as np
 from astropy.io import fits
-from astropy.table import QTable, vstack
+from astropy.table import vstack
 from pyirf.cuts import evaluate_binned_cut
 from pyirf.io import create_rad_max_hdu
 
@@ -432,29 +432,19 @@ def _make_benchmark_hdus(self, hdus):
             )
         )
         if self.do_background:
-            if not self.point_like:
-                # Create a dummy theta cut since `pyirf.sensitivity.estimate_background`
-                # needs a theta cut atm.
-                self.log.info(
-                    "Using all signal events with `theta < fov_offset_max` "
-                    "to compute the sensitivity."
-                )
-                theta_cuts = QTable()
-                theta_cuts["center"] = 0.5 * (
-                    self.sensitivity.reco_energy_bins[:-1]
-                    + self.sensitivity.reco_energy_bins[1:]
+            if self.opt_result.theta_cuts is None:
+                raise ValueError(
+                    "Calculating the point-source sensitivity requires "
+                    f"theta cuts, but {self.cuts_file} does not contain any."
                 )
-                theta_cuts["cut"] = self.sensitivity.fov_offset_max
-            else:
-                theta_cuts = self.opt_result.theta_cuts
 
             hdus.append(
                 self.sensitivity.make_sensitivity_hdu(
                     signal_events=self.signal_events[self.signal_events["selected"]],
                     background_events=self.background_events[
                         self.background_events["selected_gh"]
                     ],
-                    theta_cut=theta_cuts,
+                    theta_cut=self.opt_result.theta_cuts,
                     gamma_spectrum=self.gamma_target_spectrum,
                 )
             )

src/ctapipe/tools/optimize_event_selection.py

@@ -4,7 +4,7 @@
 from astropy.table import vstack
 
 from ..core import Provenance, Tool, traits
-from ..core.traits import AstroQuantity, Bool, Integer, classes_with_traits, flag
+from ..core.traits import AstroQuantity, Integer, classes_with_traits
 from ..irf import EventLoader, Spectra
 from ..irf.optimize import CutOptimizerBase
 
@@ -94,14 +94,6 @@ class EventSelectionOptimizer(Tool):
         help="The cut optimization algorithm to be used.",
     ).tag(config=True)
 
-    point_like = Bool(
-        False,
-        help=(
-            "Compute a theta cut in addition to the G/H separation cut "
-            "for a point-like IRF."
-        ),
-    ).tag(config=True)
-
     aliases = {
         "gamma-file": "EventSelectionOptimizer.gamma_file",
         "proton-file": "EventSelectionOptimizer.proton_file",
@@ -110,15 +102,6 @@ class EventSelectionOptimizer(Tool):
         "chunk_size": "EventSelectionOptimizer.chunk_size",
     }
 
-    flags = {
-        **flag(
-            "point-like",
-            "EventSelectionOptimizer.point_like",
-            "Compute a theta cut and a G/H separation cut.",
-            "Compute only a G/H separation cut.",
-        )
-    }
-
     classes = [EventLoader] + classes_with_traits(CutOptimizerBase)
 
     def setup(self):
@@ -229,7 +212,6 @@ def start(self):
             else None,
             precuts=self.particles[0].epp,  # identical precuts for all particle types
             clf_prefix=self.particles[0].epp.gammaness_classifier,
-            point_like=self.point_like,
         )
         self.result = result
 

src/ctapipe/tools/tests/test_compute_irf.py

@@ -19,7 +19,6 @@ def dummy_cuts_file(
 ):
     from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer
 
-    # Do "point-like" cuts to have both g/h and theta cuts in the file
     output_path = irf_tmp_path / "dummy_cuts.fits"
     run_tool(
         EventSelectionOptimizer(),
@@ -30,7 +29,6 @@ def dummy_cuts_file(
             f"--electron-file={gamma_diffuse_full_reco_file}",
             f"--output={output_path}",
             f"--config={event_loader_config_path}",
-            "--point-like",
         ],
     )
     return output_path
@@ -181,24 +179,16 @@ def test_point_like_irf_no_theta_cut(
     gamma_diffuse_full_reco_file,
     proton_full_reco_file,
     event_loader_config_path,
+    dummy_cuts_file,
     tmp_path,
 ):
+    from ctapipe.irf import OptimizationResult
     from ctapipe.tools.compute_irf import IrfTool
-    from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer
 
     gh_cuts_path = tmp_path / "gh_cuts.fits"
-    # Without the "--point-like" flag only G/H cuts are produced.
-    run_tool(
-        EventSelectionOptimizer(),
-        argv=[
-            f"--gamma-file={gamma_diffuse_full_reco_file}",
-            f"--proton-file={proton_full_reco_file}",
-            # Use diffuse gammas weighted to electron spectrum as stand-in
-            f"--electron-file={gamma_diffuse_full_reco_file}",
-            f"--output={gh_cuts_path}",
-            f"--config={event_loader_config_path}",
-        ],
-    )
+    cuts = OptimizationResult.read(dummy_cuts_file)
+    cuts.theta_cuts = None
+    cuts.write(gh_cuts_path)
     assert gh_cuts_path.exists()
 
     output_path = tmp_path / "irf.fits.gz"

src/ctapipe/tools/tests/test_optimize_event_selection.py

@@ -9,13 +9,11 @@
 pytest.importorskip("pyirf")
 
 
-@pytest.mark.parametrize("point_like", (True, False))
 def test_cuts_optimization(
     gamma_diffuse_full_reco_file,
     proton_full_reco_file,
     event_loader_config_path,
     tmp_path,
-    point_like,
 ):
     from ctapipe.irf import (
         OptimizationResult,
@@ -33,9 +31,6 @@ def test_cuts_optimization(
         f"--output={output_path}",
         f"--config={event_loader_config_path}",
     ]
-    if point_like:
-        argv.append("--point-like")
-
     ret = run_tool(EventSelectionOptimizer(), argv=argv)
     assert ret == 0
 
@@ -47,16 +42,14 @@ def test_cuts_optimization(
     assert isinstance(result.gh_cuts, QTable)
     assert result.clf_prefix == "ExtraTreesClassifier"
     assert "cut" in result.gh_cuts.colnames
-    if point_like:
-        assert isinstance(result.theta_cuts, QTable)
-        assert "cut" in result.theta_cuts.colnames
+    assert isinstance(result.theta_cuts, QTable)
+    assert "cut" in result.theta_cuts.colnames
 
     for c in ["low", "center", "high"]:
         assert c in result.gh_cuts.colnames
         assert result.gh_cuts[c].unit == u.TeV
-        if point_like:
-            assert c in result.theta_cuts.colnames
-            assert result.theta_cuts[c].unit == u.TeV
+        assert c in result.theta_cuts.colnames
+        assert result.theta_cuts[c].unit == u.TeV
 
 
 def test_cuts_opt_no_electrons(