Reworked how initial theta cuts are calculated, changed some logging …

…printouts
cta-observatory · Sep 19, 2023 · 2d5ce22 · 2d5ce22
1 parent b6eb6ff
commit 2d5ce22
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 12 deletions.
diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py
@@ -74,18 +74,26 @@ def reco_energy_bins(self):
         return reco_energy
 
     def optimise_gh_cut(self, signal, background, alpha, max_bg_radius):
-        INITIAL_GH_CUT = np.quantile(
-            signal["gh_score"], (1 - self.initial_gh_cut_efficency)
-        )
-        self.log.info(
-            f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts"
+        initial_gh_cuts = calculate_percentile_cut(
+            signal["gh_score"],
+            signal["reco_energy"],
+            bins=self.reco_energy_bins(),
+            fill_value=0.0,
+            percentile=100 * (1 - self.initial_gh_cut_efficency),
+            min_events=25,
+            smoothing=1,
         )
 
-        mask_theta_cuts = signal["gh_score"] >= INITIAL_GH_CUT
+        initial_gh_mask = evaluate_binned_cut(
+            signal["gh_score"],
+            signal["reco_energy"],
+            initial_gh_cuts,
+            op=operator.gt,
+        )
 
         theta_cuts = calculate_percentile_cut(
-            signal["theta"][mask_theta_cuts],
-            signal["reco_energy"][mask_theta_cuts],
+            signal["theta"][initial_gh_mask],
+            signal["reco_energy"][initial_gh_mask],
             bins=self.reco_energy_bins(),
             min_value=self.theta_min_angle * u.deg,
             max_value=self.theta_max_angle * u.deg,
@@ -325,7 +333,7 @@ class DataBinning(Component):
 
     fov_offset_max = Float(
         help="Maximum value for FoV offset bins in degrees",
-        default_value=2.0,
+        default_value=5.0,
     ).tag(config=True)
 
     fov_offset_n_edges = Integer(

diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py
@@ -233,11 +233,11 @@ def load_preselected_events(self):
                 len(reduced_events["electron"]),
             )
         )
-        self.log.debug(self.epp.to_table())
         select_fov = (
             reduced_events["gamma"]["true_source_fov_offset"]
             <= self.bins.fov_offset_max * u.deg
         )
+        # TODO: verify that this fov cut on only gamma is ok
         self.signal_events = reduced_events["gamma"][select_fov]
         self.background_events = vstack(
             [reduced_events["proton"], reduced_events["electron"]]
@@ -319,15 +319,15 @@ def finish(self):
             fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"),
         ]
 
+        self.log.debug("True Energy bins", self.true_energy_bins)
+        self.log.debug("FoV offset bins", self.fov_offset_bins)
         for label, mask in masks.items():
             effective_area = effective_area_per_energy_and_fov(
                 self.signal_events[mask],
                 self.sim_info,
                 true_energy_bins=self.true_energy_bins,
                 fov_offset_bins=self.fov_offset_bins,
             )
-            self.log.debug(self.true_energy_bins)
-            self.log.debug(self.fov_offset_bins)
             hdus.append(
                 create_aeff2d_hdu(
                     effective_area[..., np.newaxis],  # +1 dimension for FOV offset