Modified Example Notebooks with new ori file

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)

cosipy/threeml/COSILike.py

@@ -61,9 +61,11 @@
         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 @@
         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 @@
                         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 @@
 
         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
 

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,

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)