Use galactic coordinates with TOAST2 destriper

.gitignore

@@ -1,3 +1,29 @@
+### PyCharm ###
+# Covers JetBrains IDEs: IntelliJ, RubyMine, PhpStorm, AppCode, PyCharm, CLion, Android Studio, WebStorm and Rider
+# Reference: https://intellij-support.jetbrains.com/hc/en-us/articles/206544839
+
+# User-specific stuff
+.idea/**/workspace.xml
+.idea/**/tasks.xml
+.idea/**/usage.statistics.xml
+.idea/**/dictionaries
+.idea/**/shelf
+
+# AWS User-specific
+.idea/**/aws.xml
+
+# Generated files
+.idea/**/contentModel.xml
+
+# Sensitive or high-churn files
+.idea/**/dataSources/
+.idea/**/dataSources.ids
+.idea/**/dataSources.local.xml
+.idea/**/sqlDataSources.xml
+.idea/**/dynamic.xml
+.idea/**/uiDesigner.xml
+.idea/**/dbnavigator.xml
+
 *~
 
 # Byte-compiled / optimized / DLL files

CHANGELOG.md

@@ -8,6 +8,8 @@
 
     -   Each `Simulation` object creates random number generators (field `Simulation.random`), in a way that is safe even for MPI applications
 
+-   Support the production of maps in Galactic coordinates through the TOAST2 map-maker
+
 -   Make `make_bin_map` compute pixel indices instead of requiring them as input [#176](https://github.com/litebird/litebird_sim/pull/176)
 
 -   Use a more robust algorithm to compute pointings [#175](https://github.com/litebird/litebird_sim/pull/175)

litebird_sim/destriper/__init__.py

@@ -15,6 +15,8 @@
 from toast.tod.interval import Interval
 import toast.mpi
 
+from litebird_sim.coordinates import CoordinateSystem, rotate_coordinates_e2g
+
 toast.mpi.use_mpi = lbs.MPI_ENABLED
 
 
@@ -26,6 +28,10 @@ class DestriperParameters:
 
     - ``nside``: the NSIDE parameter used to create the maps
 
+    - ``coordinate_system``: an instance of the :class:`.CoordinateSystem` enum.
+      It specifies if the map must be created in ecliptic (default) or
+      galactic coordinates.
+
     - ``nnz``: number of components per pixel. The default is 3 (I/Q/U).
 
     - ``baseline_length``: number of consecutive samples in a 1/f noise
@@ -64,6 +70,7 @@ class DestriperParameters:
     """
 
     nside: int = 512
+    coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic
     nnz: int = 3
     baseline_length: int = 100
     iter_max: int = 100
@@ -104,12 +111,20 @@ class DestriperResult:
     npp: Any = None
     invnpp: Any = None
     rcond: Any = None
+    coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic
 
 
 class _Toast2FakeCache:
     "This class simulates a TOAST2 cache"
 
-    def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside):
+    def __init__(
+        self,
+        spin2ecliptic_quats,
+        obs,
+        bore2spin_quat,
+        nside,
+        coordinates: CoordinateSystem,
+    ):
         self.obs = obs
 
         self.keydict = {"timestamps": obs.get_times()}
@@ -141,13 +156,27 @@ def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside):
                 logging.warning(
                     "converting TODs for %s from %s to float64",
                     obs.name[i],
-                    str(pointings[i].dtype),
+                    str(obs.tod[i].dtype),
                 )
                 self.keydict[f"signal_{det}"] = np.array(obs.tod[i], dtype=np.float64)
 
-            self.keydict[f"pixels_{det}"] = healpix_base.ang2pix(curpnt[:, 0:2])
+            theta_phi = curpnt[:, 0:2]
+            polangle = curpnt[:, 2]
+
+            if coordinates == CoordinateSystem.Galactic:
+                theta_phi, polangle = rotate_coordinates_e2g(
+                    pointings_ecl=theta_phi, pol_angle_ecl=polangle
+                )
+            elif coordinates == CoordinateSystem.Ecliptic:
+                pass  # Do nothing, "theta_phi" and "polangle" are ok
+            else:
+                assert ValueError(
+                    "unable to handle coordinate system {coordinates} in `destripe`"
+                )
+
+            self.keydict[f"pixels_{det}"] = healpix_base.ang2pix(theta_phi)
             self.keydict[f"weights_{det}"] = np.stack(
-                (np.ones(nsamples), np.cos(2 * curpnt[:, 2]), np.sin(2 * curpnt[:, 2]))
+                (np.ones(nsamples), np.cos(2 * polangle), np.sin(2 * polangle))
             ).transpose()
 
     def keys(self):
@@ -172,10 +201,19 @@ def destroy(self, name):
 class _Toast2FakeTod:
     "This class simulates a TOAST2 TOD"
 
-    def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside):
+    def __init__(
+        self,
+        spin2ecliptic_quats,
+        obs,
+        bore2spin_quat,
+        nside,
+        coordinates: CoordinateSystem,
+    ):
         self.obs = obs
         self.local_samples = (0, obs.tod[0].size)
-        self.cache = _Toast2FakeCache(spin2ecliptic_quats, obs, bore2spin_quat, nside)
+        self.cache = _Toast2FakeCache(
+            spin2ecliptic_quats, obs, bore2spin_quat, nside, coordinates
+        )
 
     def local_intervals(self, _):
         return [
@@ -215,9 +253,20 @@ def local_times(self):
 class _Toast2FakeData:
     "This class simulates a TOAST2 Data class"
 
-    def __init__(self, spin2ecliptic_quats, obs, bore2spin_quat, nside):
+    def __init__(
+        self,
+        spin2ecliptic_quats,
+        obs,
+        bore2spin_quat,
+        nside,
+        coordinates: CoordinateSystem,
+    ):
         self.obs = [
-            {"tod": _Toast2FakeTod(spin2ecliptic_quats, x, bore2spin_quat, nside)}
+            {
+                "tod": _Toast2FakeTod(
+                    spin2ecliptic_quats, x, bore2spin_quat, nside, coordinates
+                )
+            }
             for x in obs
         ]
         self.bore2spin_quat = bore2spin_quat
@@ -285,6 +334,7 @@ def destripe_observations(
         obs=observations,
         bore2spin_quat=bore2spin_quat,
         nside=params.nside,
+        coordinates=params.coordinate_system,
     )
     mapmaker = OpMapMaker(
         nside=params.nside,
@@ -306,6 +356,8 @@ def destripe_observations(
 
     result = DestriperResult()
 
+    result.coordinate_system = params.coordinate_system
+
     if params.return_hit_map:
         result.hit_map = healpy.read_map(
             base_path / (params.output_file_prefix + "hits.fits"),

litebird_sim/scan_map.py

@@ -140,7 +140,7 @@ def scan_map_in_observations(
             elif all(item in maps.keys() for item in cur_obs.channel):
                 input_names = cur_obs.channel
             else:
-                print(
+                raise ValueError(
                     "The dictionary maps does not contain all the relevant"
                     + "keys, please check the list of detectors and channels"
                 )