SAA Passage

cosipy/data_io/UnBinnedData.py

@@ -606,9 +606,9 @@
 
         return this_dict
 
-    def select_data(self, output_name=None, unbinned_data=None):
+    def select_data_time(self, output_name=None, unbinned_data=None):
 
-        """Applies cuts to unbinnned data dictionary. 
+        """Applies time cuts to unbinnned data dictionary. 
 
         Parameters
         ----------
@@ -617,10 +617,6 @@
         output_name : str, optional
             Prefix of output file (default is None, in which case no 
             file is saved).
-
-        Note
-        ----
-        Only cuts in time are allowed for now. 
         """
 
         logger.info("Making data selections...")
@@ -685,3 +681,106 @@
             self.write_unbinned_output(output_name)
 
         return
+
+    def find_bad_intervals(self, times, values, bad_value=0.0):
+
+        """Finds intervals where livetime is 0.0.
+
+        Parameters
+        ----------
+        times : array
+            Array of times from ori file.
+        values : array
+            Array of livetimes from ori file.        
+        bad_value : float or int
+            The value that defines a bad time interval. It must match 
+            exactly the value in the ori file, including the 
+            type (float or int). Default is 0.0. 
+
+        Returns
+        -------
+        bad_intervals : list of tuples
+            List of bad time intervals.
+        """
+
+        bad_intervals = []
+        start = None
+
+        for i in range(len(times)-1):
+            if values[i] == bad_value:
+                if start is None:
+                    start = times[i]
+            else:
+                if start is not None:
+                    bad_intervals.append((start.value, times[i].value))
+                    start = None
+
+        if start is not None:
+            bad_intervals.append((start.value, times[-1].value))
+
+        return bad_intervals
+
+    def filter_good_data(self, times, bad_intervals):
+
+        """Removes entries that fall within bad intervals.
+
+        Parameters
+        ----------
+        times : array
+            Array of photon event times.
+        bad_intervals : list
+            List of bad intervals defined by livetime = 0.0.
+
+        Returns
+        -------
+        filtered_index : list
+            List of indices of good events. 
+        """
+
+        filtered_index = []
+
+        # Get indices for good times.
+        # The inequality below corresponds to left end point in the ori file.
+        for i in range(len(times)-1):
+            if not any(start < times[i] < end for start, end in bad_intervals):
+                filtered_index.append(i)
+
+        return filtered_index
+
+    def cut_SAA_events(self, unbinned_data=None, output_name=None):
+
+        """Cuts events corresponding to SAA passage based on input ori file.
+
+        Parameters
+        ----------
+        unbinned_data : str, optional
+            Name of unbinned dictionary file.
+        output_name : str, optional
+            Prefix of output file (default is None, in which case no 
+            file is saved).
+        """
+
+        # Option to read in unbinned data file:
+        if unbinned_data:
+            self.cosi_dataset = self.get_dict(unbinned_data)
+
+        # Get ori info:
+        ori = SpacecraftFile.parse_from_file(self.ori_file)
+
+        # Get bad time intervals:
+        bti = self.find_bad_intervals(ori._time, ori.livetime)
+
+        # Get indices for good photons:
+        time_keep_index = self.filter_good_data(self.cosi_dataset['TimeTags'], bti)
+
+        # Apply cuts to dictionary:
+        for key in self.cosi_dataset:
+
+            self.cosi_dataset[key] = self.cosi_dataset[key][time_keep_index]
+
+        # Write unbinned data to file (either fits or hdf5):
+        if output_name != None:
+            logger.info("Saving file...")
+            self.write_unbinned_output(output_name)
+
+        return

cosipy/spacecraftfile/SpacecraftFile.py

@@ -22,7 +22,8 @@ class SpacecraftFile():
 
     def __init__(self, time, x_pointings = None, y_pointings = None, \
             z_pointings = None, earth_zenith = None, altitude = None,\
-            attitude = None, instrument = "COSI", frame = "galactic"):
+            attitude = None, livetime = None, instrument = "COSI", \
+            frame = "galactic"):
 
         """
         Handles the spacecraft orientation. Calculates the dwell time 
@@ -48,8 +49,12 @@ def __init__(self, time, x_pointings = None, y_pointings = None, \
         earth_zenith : astropy.coordinates.SkyCoord, optional
             The pointings (galactic system) of the Earth zenith (the 
             default is `None`, which implies no input for the earth pointings).
-	altitude : array, optional 
+	    altitude : array, optional 
             Altitude of the spacecraft in km.
+        livetime : array, optional 
+            Time in seconds the instrument is live for the corresponding 
+            energy bin (using left endpoints so that the last entry in 
+            the ori file is 0).
         attitude : numpy.ndarray, optional 
             The attitude of the spacecraft (the default is `None`, 
             which implies no input for the attitude of the spacecraft).
@@ -70,6 +75,10 @@ def __init__(self, time, x_pointings = None, y_pointings = None, \
         if not isinstance(altitude, (type(None))):
             self._altitude = np.array(altitude)
 
+        # livetime
+        if not isinstance(livetime, (type(None))):
+            self.livetime = np.array(livetime)
+
         # x pointings
         if isinstance(x_pointings, (SkyCoord, type(None))):
             self.x_pointings = x_pointings
@@ -131,19 +140,21 @@ def parse_from_file(cls, file):
             The SpacecraftFile object.
         """
 
-        orientation_file = np.loadtxt(file, usecols=(1, 2, 3, 4, 5, 6, 7, 8),delimiter=' ', skiprows=1, comments=("#", "EN"))
+        orientation_file = np.loadtxt(file, usecols=(1, 2, 3, 4, 5, 6, 7, 8, 9),delimiter=' ', skiprows=1, comments=("#", "EN"))
         time_stamps = orientation_file[:, 0]
         axis_1 = orientation_file[:, [2, 1]]
         axis_2 = orientation_file[:, [4, 3]]
         axis_3 = orientation_file[:, [7, 6]]
         altitude = np.array(orientation_file[:, 5]) 
-
+        livetime = np.array(orientation_file[:, 8])
+        livetime = livetime[:-1] # left end points, so remove last bin. 
+
         time = Time(time_stamps, format = "unix")
         xpointings = SkyCoord(l = axis_1[:,0]*u.deg, b = axis_1[:,1]*u.deg, frame = "galactic")
         zpointings = SkyCoord(l = axis_2[:,0]*u.deg, b = axis_2[:,1]*u.deg, frame = "galactic")
         earthpointings = SkyCoord(l = axis_3[:,0]*u.deg, b = axis_3[:,1]*u.deg, frame = "galactic")
 
-        return cls(time, x_pointings = xpointings, z_pointings = zpointings, earth_zenith = earthpointings, altitude = altitude)
+        return cls(time, x_pointings = xpointings, z_pointings = zpointings, earth_zenith = earthpointings, altitude = altitude, livetime=livetime)
 
     def get_time(self, time_array = None):
 
@@ -577,7 +588,8 @@ def get_scatt_map(self,
         earth_occ_index = src_angle.value >= max_angle
 
         # Define weights and set to 0 if blocked by Earth:
-        weight = np.diff(timestamps.gps)*u.s
+        weight = self.livetime*u.s
+
         if earth_occ == True:
             weight[earth_occ_index[:-1]] = 0        
 

cosipy/test_data/20280301_2s.ori

@@ -1,5 +1,5 @@
 Type OrientationsGalactic 
-OG 1835487300.0 68.44034002307066 44.61117227945379 -21.559659976929343 44.61117227945379 0.0 0.0 0.0
-OG 1835487301.0 68.38920658776064 44.642124027021296 -21.610793412239364 44.642124027021296 0.0 0.0 0.0
-OG 1835487302.0 68.3380787943012 44.67309722321497 -21.661921205698793 44.67309722321497 0.0 0.0 0.0
+OG 1835487300.0 68.44034002307066 44.61117227945379 -21.559659976929343 44.61117227945379 0.0 0.0 0.0 1.0
+OG 1835487301.0 68.38920658776064 44.642124027021296 -21.610793412239364 44.642124027021296 0.0 0.0 0.0 1.0 
+OG 1835487302.0 68.3380787943012 44.67309722321497 -21.661921205698793 44.67309722321497 0.0 0.0 0.0 0.0
 EN

cosipy/test_data/20280301_first_10sec.ori

@@ -1,13 +1,13 @@
 Type OrientationsGalactic 
-OG 1835478000.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0
-OG 1835478001.0 73.09517926980278 41.88225011209611 16.904820730197223 221.8822501120961 0.0 0.0 0.0
-OG 1835478002.0 73.04128380352786 41.90629597072256 16.95871619647214 221.90629597072257 0.0 0.0 0.0
-OG 1835478003.0 72.98739108131268 41.93035532675578 17.012608918687327 221.93035532675577 0.0 0.0 0.0
-OG 1835478004.0 72.9335011165853 41.954428243823145 17.066498883414702 221.95442824382317 0.0 0.0 0.0
-OG 1835478005.0 72.87961392277379 41.978514785552235 17.120386077226204 221.97851478555222 0.0 0.0 0.0
-OG 1835478006.0 72.82572951330626 42.002615015570285 17.174270486693747 222.0026150155703 0.0 0.0 0.0
-OG 1835478007.0 72.77184790161073 42.02672899750497 17.228152098389273 222.02672899750493 0.0 0.0 0.0
-OG 1835478008.0 72.7179691011153 42.05085679498347 17.282030898884702 222.05085679498347 0.0 0.0 0.0
-OG 1835478009.0 72.66409312524804 42.07499847163346 17.335906874751963 222.07499847163342 0.0 0.0 0.0
-OG 1835478010.0 72.61021998743702 42.09915409108222 17.38978001256298 222.09915409108223 0.0 0.0 0.0
+OG 1835478000.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0 1.0
+OG 1835478001.0 73.09517926980278 41.88225011209611 16.904820730197223 221.8822501120961 0.0 0.0 0.0 1.0 
+OG 1835478002.0 73.04128380352786 41.90629597072256 16.95871619647214 221.90629597072257 0.0 0.0 0.0 1.0
+OG 1835478003.0 72.98739108131268 41.93035532675578 17.012608918687327 221.93035532675577 0.0 0.0 0.0 1.0 
+OG 1835478004.0 72.9335011165853 41.954428243823145 17.066498883414702 221.95442824382317 0.0 0.0 0.0 1.0 
+OG 1835478005.0 72.87961392277379 41.978514785552235 17.120386077226204 221.97851478555222 0.0 0.0 0.0 1.0
+OG 1835478006.0 72.82572951330626 42.002615015570285 17.174270486693747 222.0026150155703 0.0 0.0 0.0 1.0
+OG 1835478007.0 72.77184790161073 42.02672899750497 17.228152098389273 222.02672899750493 0.0 0.0 0.0 1.0
+OG 1835478008.0 72.7179691011153 42.05085679498347 17.282030898884702 222.05085679498347 0.0 0.0 0.0 1.0 
+OG 1835478009.0 72.66409312524804 42.07499847163346 17.335906874751963 222.07499847163342 0.0 0.0 0.0 1.0
+OG 1835478010.0 72.61021998743702 42.09915409108222 17.38978001256298 222.09915409108223 0.0 0.0 0.0 0.0
 EN

cosipy/test_data/polarization_ori.ori

@@ -1,4 +1,4 @@
 Type OrientationsGalactic 
-OG 0.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0
-OG 15.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0
-EN
+OG 0.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0 15.0 
+OG 15.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0 0.0
+EN

docs/tutorials/DataIO/DataIO_example.ipynb

@@ -1075,7 +1075,7 @@
    ],
    "source": [
     "analysis = BinnedData(\"inputs_half_time.yaml\")\n",
-    "analysis.select_data(unbinned_data=\"combined_unbinned_data.hdf5\", output_name=\"selected_unbinned_data\")"
+    "analysis.select_data_time(unbinned_data=\"combined_unbinned_data.hdf5\", output_name=\"selected_unbinned_data\")"
    ]
   },
   {

docs/tutorials/background_estimation/continuum_estimation/BG_estimation_example.ipynb

@@ -9,7 +9,7 @@
    "source": [
     "# Continuum Background Estimation\n",
     "\n",
-    "This example shows how to use cosipy for estimating the continuum background. For more details on how the algorithm works, see the other notebook in this directory. In short, the method is based on the traditional on-off analysis, where the background for a source region on the sky is estimated by performing some kind of interpolation of the nearby surrounding region. The main difference with a Compton telescope is that we are now performing the on-off analysis in the Compton data space.\n",
+    "This example shows how to use cosipy for estimating the continuum background. In short, the method is based on the traditional on-off analysis, where the background for a source region on the sky is estimated by performing some kind of interpolation of the nearby surrounding region. The main difference with a Compton telescope is that we are now performing the on-off analysis in the Compton data space.\n",
     "\n",
     "In this particular example, we want to estimate the background for the Crab. We start with the full dataset. This contains the full background, which includes instrumental + astrophysical, where the latter consists of all astrophysical sources other than the Crab. Then, for each bin of Em and Phi, we mask the Crab in the PsiChi plane based on the point source response. Finally, we interpolate over the masked region using an inpainting method.   \n",
     "\n",
@@ -586,9 +586,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "COSI",
+   "display_name": "Python [conda env:COSIPY]",
    "language": "python",
-   "name": "cosi"
+   "name": "conda-env-COSIPY-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -600,7 +600,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.14"
+   "version": "3.10.15"
   }
  },
  "nbformat": 4,

tests/dataIO/test_binned_data.py

@@ -24,7 +24,7 @@ def test_binned_data(tmp_path):
     # Read in test file and select first 0.1 seconds to reduce number of photons:
     analysis.read_tra()
     analysis.tmax = 1835478000.1
-    analysis.select_data(output_name=tmp_path/"temp_unbinned")
+    analysis.select_data_time(output_name=tmp_path/"temp_unbinned")
 
     # Test default binning (which is in Galactic coordinates).
     # Use small number of bins to speed things up. 

tests/dataIO/test_unbinned_data_all.py

@@ -86,12 +86,15 @@ def test_unbinned_data_all(tmp_path):
     # Also test reading in the hdf5 file:
     analysis.unbinned_output = "hdf5"
     analysis.tmax = 1835478001.0
-    analysis.select_data(unbinned_data=tmp_path/"test_h5.hdf5",output_name=tmp_path/"temp_test_file")
+    analysis.select_data_time(unbinned_data=tmp_path/"test_h5.hdf5",output_name=tmp_path/"temp_test_file")
     assert np.amax(analysis.cosi_dataset['TimeTags']) <= 1835478001.0    
 
     # Test reading tra with no pointing info:
     analysis.data_file = os.path.join(test_data.path,\
             "GalacticScan.inc1.id1.crab10sec.extracted.testsample.nopointinginfo.tra.gz")
     analysis.read_tra()
 
+    # Test SAA cut:
+    analysis.cut_SAA_events(unbinned_data=tmp_path/"test_h5.hdf5")
+
     return