Merge pull request #61 from VERITAS-Observatory/on_exclusion

feat(background_models.py): allowing for user defined default exclusion

gammapy_tools/make_background/background_models.py

@@ -17,6 +17,7 @@
 # from here
 from ..utils import ExclusionFinder
 
+from typing import Optional
 
 class BackgroundModelEstimator:
     """
@@ -33,6 +34,7 @@ def __init__(
         excluded_sources: list = [],
         smooth_sigma: float = 1,
         njobs: int = 0,
+        default_exclusion: Optional[float] = None,
     ) -> None:
         """BackgroundModelEstimator class for estimating the background from runs
 
@@ -45,6 +47,8 @@ def __init__(
             excluded_sources (list)             - list of sources to be excluded
                                                  (astropy.coordinates.Skycoords)
             njobs (int)                         - Number of jobs to run
+            default_exclusion (float)           - Default exclusion radius for sources. Defaults to None
+                                                 which will use the exclusion finder to determine the exclusion radius
 
         Returns
         ----------
@@ -85,7 +89,7 @@ def __init__(
         # )
         # self.hawc = SourceCatalog3HWC("$GAMMAPY_DATA/catalogs/3HWC.ecsv")
 
-        self.exclusion = ExclusionFinder()
+        self.exclusion = ExclusionFinder(default_exclusion = default_exclusion)
 
     @staticmethod
     def _make_bkg2d(energy: MapAxis, offset: MapAxis, unit: str) -> Background2D:

gammapy_tools/make_background/make_background.py

@@ -131,12 +131,18 @@ def get_background_for_run(parms: tuple[float, dict]) -> tuple[str, list]:
         if "smooth_sigma" in config["background_selection"]:
             smooth_sigma = config["background_selection"]["smooth_sigma"]
 
+        default_exclusion = None 
+        if "sky_map" in config:
+            if "on_exclusion_region" in config["sky_map"]:
+                default_exclusion = config["sky_map"]["on_exclusion_region"]
+
         estimator = BackgroundModelEstimator(
             energy,
             offset,
             smooth=config["background_selection"]["smooth"],
             smooth_sigma=smooth_sigma,
             njobs=config["config"]["njobs"],
+            default_exclusion=default_exclusion,
         )
 
         estimator.run(observations)
@@ -213,8 +219,16 @@ def generate_background_from_run(parms: tuple[int, dict]) -> str:
             name="offset",
         )
 
+        default_exclusion = None 
+        if "sky_map" in config:
+            if "on_exclusion_region" in config["sky_map"]:
+                default_exclusion = config["sky_map"]["on_exclusion_region"]
+
+
         estimator = BackgroundModelEstimator(
-            energy, offset, smooth=config["background_selection"]["smooth"]
+            energy, offset, 
+            smooth=config["background_selection"]["smooth"],
+            default_exclusion=default_exclusion,
         )
         estimator.run(observations)
         if "BACKGROUND" in hdul:

gammapy_tools/utils/exclusion_finder.py

@@ -4,17 +4,24 @@
 from gammapy.catalog import SourceCatalog3HWC, SourceCatalogGammaCat
 
 
+from typing import Optional
+
 class ExclusionFinder:
     """
     Class for finding exclusion regions for the analysis
 
     """
 
-    def __init__(self):
+    def __init__(self, default_exclusion: Optional[float] = 0.35):
         """Initilization function
 
         Load in the star catalog, GammaCat and 3HWC
+
+        Parameters
+        ----------
+            default_exclusion (float)          - Default theta cut to use for exclusion region
         """
+        self.default_exclusion = default_exclusion
 
         # ToDo wrap all this into package info
         this_dir, _ = path.split(__file__)
@@ -141,8 +148,9 @@ def find_exclusion(
         ra: float,
         dec: float,
         theta: float = 2.0,
-        theta_cut: float = 0.35,
+        theta_cut: Optional[float] = None,
         mag_cut: float = 7,
+        star_theta_cut: float = 0.35,
     ) -> tuple[list[tuple[float, float, float]], list[str]]:
         """Find sources to exclude
 
@@ -153,10 +161,11 @@ def find_exclusion(
             theta (float)                       - Search radius (in degrees)
                                                   Defaults to 2.0
             theta_cut (float)                   - Default theta cut to use for exclusion region
-                                                  Defaults to 0.35
+                                                  Defaults to 0.35 or the default_exclusion
             mag_cut (float)                     - B Magnitude cut on star brightness
                                                   Defaults to mag 7
-
+            star_theta_cut (float)              - Theta cut for stars
+                                                  Defaults to 0.35
 
         Returns
         ----------
@@ -167,10 +176,13 @@ def find_exclusion(
         regions = []
         sources_excluded = []
 
+        if theta_cut is None:
+            theta_cut = self.default_exclusion
+
         # Get the stars in the region
         stars = self.find_stars(ra, dec, theta, mag_cut)
         for star in stars:
-            regions.append((star["ra"], star["dec"], theta_cut))
+            regions.append((star["ra"], star["dec"], star_theta_cut))
             sources_excluded.append(f"star_{star['id']:0.0f}")
 
         # Find the HAWC sources