Merge pull request #264 from ckarwin/develop-Earth_Occultation

Modified Example Notebooks with new ori file
cositools · Nov 14, 2024 · 5ea0243 · 5ea0243
2 parents f71e55d + 5b693b5
commit 5ea0243
Show file tree

Hide file tree

Showing 7 changed files with 316 additions and 420 deletions.
diff --git a/cosipy/spacecraftfile/SpacecraftFile.py b/cosipy/spacecraftfile/SpacecraftFile.py
@@ -518,6 +518,7 @@ def get_scatt_map(self,
                        scheme = 'ring',
                        coordsys = 'galactic',
                        r_earth = 6378.0,
+                       earth_occ = True
                        ):
 
         """
@@ -537,7 +538,10 @@ def get_scatt_map(self,
             The coordinate system used in the scatt map (the default is "galactic).
         r_earth : float, optional
             Earth radius in km (default is 6378 km).
-        
+        earth_occ : bool, optional
+            Option to include Earth occultation in scatt map calculation.
+            Default is True. 
+
         Returns
         -------
         h_ori : cosipy.spacecraftfile.scatt_map.SpacecraftAttitudeMap
@@ -574,7 +578,8 @@ def get_scatt_map(self,
 
         # Define weights and set to 0 if blocked by Earth:
         weight = np.diff(timestamps.gps)*u.s
-        weight[earth_occ_index[:-1]] = 0        
+        if earth_occ == True:
+            weight[earth_occ_index[:-1]] = 0        
 
         # Fill histogram:
         h_ori.fill(x, y, weight = weight)

diff --git a/cosipy/threeml/COSILike.py b/cosipy/threeml/COSILike.py
@@ -61,9 +61,11 @@ class COSILike(PluginPrototype):
         attached to them
     precomputed_psr_file : str, optional
         Full path to precomputed point source response in Galactic coordinates
+    earth_occ : bool, optional
+        Option to include Earth occultation in fit (default is True).
     """
     def __init__(self, name, dr, data, bkg, sc_orientation, 
-                 nuisance_param=None, coordsys=None, precomputed_psr_file=None, **kwargs):
+                 nuisance_param=None, coordsys=None, precomputed_psr_file=None, earth_occ=True, **kwargs):
 
         # create the hash for the nuisance parameters. We have none for now.
         self._nuisance_parameters = collections.OrderedDict()
@@ -78,6 +80,7 @@ def __init__(self, name, dr, data, bkg, sc_orientation,
         self._data = data
         self._bkg = bkg
         self._sc_orientation = sc_orientation
+        self.earth_occ = earth_occ
 
         try:
             if data.axes["PsiChi"].coordsys.name != bkg.axes["PsiChi"].coordsys.name:
@@ -194,7 +197,7 @@ def set_model(self, model):
                         dwell_time_map = self._get_dwell_time_map(coord)
                         self._psr[name] = self._dr.get_point_source_response(exposure_map=dwell_time_map)
                     elif self._coordsys == 'galactic':
-                        scatt_map = self._get_scatt_map()
+                        scatt_map = self._get_scatt_map(coord)
                         self._psr[name] = self._dr.get_point_source_response(coord=coord, scatt_map=scatt_map)
                     else:
                         raise RuntimeError("Unknown coordinate system")
@@ -325,16 +328,22 @@ def _get_dwell_time_map(self, coord):
 
         return dwell_time_map
 
-    def _get_scatt_map(self):
+    def _get_scatt_map(self, coord):
         """
         Get the spacecraft attitude map of the source in the inertial (spacecraft) frame.
         
+        Parameters
+        ----------
+        coord : astropy.coordinates.SkyCoord
+            The coordinates of the target object.
+
         Returns
         -------
         scatt_map : cosipy.spacecraftfile.scatt_map.SpacecraftAttitudeMap
         """
 
-        scatt_map = self._sc_orientation.get_scatt_map(nside = self._dr.nside * 2, coordsys = 'galactic')
+        scatt_map = self._sc_orientation.get_scatt_map(coord, nside = self._dr.nside * 2, \
+                coordsys = 'galactic', earth_occ = self.earth_occ)
 
         return scatt_map
 

diff --git a/docs/tutorials/spectral_fits/continuum_fit/crab/SpectralFit_Crab.ipynb b/docs/tutorials/spectral_fits/continuum_fit/crab/SpectralFit_Crab.ipynb
diff --git a/docs/tutorials/spectral_fits/continuum_fit/grb/SpectralFit_GRB.ipynb b/docs/tutorials/spectral_fits/continuum_fit/grb/SpectralFit_GRB.ipynb
diff --git a/docs/tutorials/spectral_fits/extended_source_fit/diffuse_511_spectral_fit.ipynb b/docs/tutorials/spectral_fits/extended_source_fit/diffuse_511_spectral_fit.ipynb
@@ -23,7 +23,7 @@
     "This tutotrial also walks through all the steps needed when performing a spectral fit, starting with the unbinned data, i.e. creating the combined data set, and binning the data. \n",
     "\n",
     "For the first two examples, you will need the following files (available on wasabi):<br>\n",
-    "**20280301_3_month.ori <br>\n",
+    "**20280301_3_month_with_orbital_info.ori <br>\n",
     "cosmic_photons_3months_unbinned_data.fits.gz <br>\n",
     "511_Testing_3months.fits.gz <br>\n",
     "SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5 <br>\n",
@@ -94,7 +94,7 @@
     "# ori file:\n",
     "# wasabi path: COSI-SMEX/DC2/Data/Orientation/20280301_3_month.ori\n",
     "# File size: 684 MB\n",
-    "fetch_wasabi_file('COSI-SMEX/DC2/Data/Orientation/20280301_3_month.ori')"
+    "fetch_wasabi_file('COSI-SMEX/DC2/Data/Orientation/20280301_3_month_with_orbital_info.ori')"
    ]
   },
   {
@@ -770,7 +770,7 @@
    "source": [
     "response_file = \"SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5\"\n",
     "response = FullDetectorResponse.open(response_file)\n",
-    "ori = SpacecraftFile.parse_from_file(\"20280301_3_month.ori\")\n",
+    "ori = SpacecraftFile.parse_from_file(\"20280301_3_month_with_orbital_info.ori\")\n",
     "psr_file = \"psr_gal_511_DC2.h5\""
    ]
   },
@@ -3157,7 +3157,7 @@
     "# if not previously loaded in example 1, load the response, ori, and psr: \n",
     "response_file = \"SMEXv12.511keV.HEALPixO4.binnedimaging.imagingresponse.nonsparse_nside16.area.h5\"\n",
     "response = FullDetectorResponse.open(response_file)\n",
-    "ori = SpacecraftFile.parse_from_file(\"20280301_3_month.ori\")\n",
+    "ori = SpacecraftFile.parse_from_file(\"20280301_3_month_with_orbital_info.ori\")\n",
     "psr_file = \"psr_gal_511_DC2.h5\""
    ]
   },
@@ -4299,9 +4299,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [conda env:COSI]",
+   "display_name": "Python [conda env:analysis]",
    "language": "python",
-   "name": "conda-env-COSI-py"
+   "name": "conda-env-analysis-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -4313,7 +4313,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.15"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,

diff --git a/tests/spectral_fits/__init__.py b/tests/spectral_fits/__init__.py
diff --git a/tests/spectral_fits/test_spectral_fitting.py → tests/threeml/test_spectral_fitting.py b/tests/spectral_fits/test_spectral_fitting.py → tests/threeml/test_spectral_fitting.py
@@ -4,6 +4,7 @@
 import numpy as np
 from threeML import Band, PointSource, Model, JointLikelihood, DataList
 from astromodels import Parameter
+from astropy.coordinates import SkyCoord
 
 data_path = test_data.path
 
@@ -70,3 +71,7 @@ def test_point_source_spectral_fit():
                        [1.0743623124061388, -1.1000643881813548, -2.299033632814098, 449.99790270666415, 1.0], atol=[0.1, 0.1, 0.1, 1.0, 0.1])
 
     assert np.allclose([cosi.get_log_like()], [337.17196587486285], atol=[1.0])
+
+    # Test scatt map method:
+    coord = SkyCoord(l=184.56*u.deg,b=-5.78*u.deg,frame="galactic")
+    cosi._get_scatt_map(coord)