Tod effect driver

CHANGELOG.md

@@ -1,4 +1,5 @@
 # HEAD
+-   Extract a shared driver for the systematic-effect modules (`for_each_observation`, `for_each_observation_with_pointings`, `normalize_observations`) and rename `pointings_in_obs.py` to `observation_utilities.py`, which now collects all the helpers that operate on a collection of `Observation` objects. Noise, gain drift, non-linearity, HWP differential emission and the CMB dipole now share one place for observation normalization, per-detector RNG resolution and TOD-component indirection instead of re-implementing it each. No change to public effect interfaces.
 
 -   Stop supporting unused MPI processes and drop `MPI_COMM_GRID` class and its attributes. Now throwing an error when `comm_size != n_blocks_det * n_blocks_time`. It fixes [#364](https://github.com/litebird/litebird_sim/issues/364), while incorporating selective changes from [#372](https://github.com/litebird/litebird_sim/pull/372) and [#427](https://github.com/litebird/litebird_sim/pull/427). Reference PR: [#539](https://github.com/litebird/litebird_sim/pull/539).
 

docs/source/h_maps.rst

@@ -38,7 +38,7 @@ These maps capture how the detector orientation is distributed across the sky
 during the scanning strategy. Setting :math:`m=0` is equivalent to the definition of h maps in McCallum et al (2021) without HWP modulation.
 
 
-The maps can be generated by using the :func:`make_h_maps` function with a list of observations, or by using the observations of a simulation through the method :func:`.make_h_maps` of :class:`simulation`
+The maps can be generated by using the :func:`~litebird_sim.mapmaking.h_maps.make_h_maps` function with a list of observations, or by using the observations of a simulation through the :meth:`.Simulation.make_h_maps` method.
 
 Example
 -----------
@@ -100,7 +100,7 @@ Note: for the couple (0,0) the resulting map is not the one expected from the de
 Output: ``HMapsResult``
 -----------------------
 
-:func:`make_h_maps` returns a :class:`HMapsResult` object, which contains:
+:func:`~litebird_sim.mapmaking.h_maps.make_h_maps` returns a :class:`~litebird_sim.mapmaking.h_maps.HMapsResult` object, which contains:
 
 - ``h_maps``: a dictionary indexed by detector name, then by ``(n, m)``
   tuple, each entry being a :class:`.h_map_Re_and_Im` object.
@@ -149,7 +149,7 @@ Maps are saved in **HDF5** format, one file per detector:
        │   └── Im
        └── ...
 
-To reload maps from disk, use :func:`.load_h_maps_from_file`::
+To reload maps from disk, use :func:`~litebird_sim.mapmaking.h_maps.load_h_maps_from_file`::
 
    from litebird_sim.mapmaking.h_maps import load_h_maps_from_file
 
@@ -165,7 +165,7 @@ and all time samples.
 
 MPI support
 -----------
-The current implementation of :func:`make_h_maps` only allows to distribute observation by detector, i.e. each MPI process computes the h maps for a subset of detectors, but using all time samples.
+The current implementation of :func:`~litebird_sim.mapmaking.h_maps.make_h_maps` only allows to distribute observation by detector, i.e. each MPI process computes the h maps for a subset of detectors, but using all time samples.
 Example:
 
 .. code-block:: python

docs/source/index.rst

@@ -20,6 +20,21 @@ This is the User’s Manual of the LiteBIRD Simulation Framework.
    part6.rst
    appendix.rst
 
+How to cite
+===========
+
+If you use the LiteBIRD Simulation Framework in your work, please cite
+:cite:`2025:litebird:tomasi`:
+
+   M. Tomasi *et al.* (LiteBIRD Collaboration), "A simulation framework for the
+   LiteBIRD instruments," *JCAP* **11** (2025) 040,
+   doi:`10.1088/1475-7516/2025/11/040
+   <https://doi.org/10.1088/1475-7516/2025/11/040>`_,
+   `arXiv:2507.04918 <https://arxiv.org/abs/2507.04918>`_.
+
+The corresponding BibTeX entry (key ``2025:litebird:tomasi``) is available in
+:doc:`the bibliography <bibliography>`.
+
 Indices and tables
 ==================
 

docs/source/part1.rst

@@ -1,6 +1,22 @@
 Introduction
 ============
 
+`LiteBIRD <https://www.litebird.jp/eng/>`_ is a JAXA-led satellite mission that
+will map the polarization of the Cosmic Microwave Background over the full sky,
+searching for the primordial *B*-mode signature of cosmic inflation. Meeting its
+requirements relies on detailed end-to-end simulations of the instruments and of
+the observation process.
+
+The **LiteBIRD Simulation Framework** (LBS) is the collaboration's toolkit for
+those simulations: it turns *sky models* and *instrument descriptions* into
+*time-ordered detector data* (TOD), injects instrumental noise and systematic
+effects, and reduces the TOD back into sky maps. The rest of this manual
+documents these steps part by part.
+
+If you are new to the framework, start with the installation instructions below
+and then work through the :doc:`tutorial <tutorial>`. If you use LBS in your
+work, please cite it (:cite:`2025:litebird:tomasi`).
+
 .. toctree::
    :maxdepth: 1
 

docs/source/part2.rst

@@ -1,6 +1,11 @@
 Structure of the framework
 ==========================
 
+This part describes the core building blocks of the framework: the central
+``Simulation`` class, the objects it ties together (observations, detectors, the
+instrument model), and the cross-cutting facilities shared by every module —
+units, reports, reproducible random numbers, and MPI parallelism.
+
 .. toctree::
    :maxdepth: 1
 

docs/source/part3.rst

@@ -1,6 +1,10 @@
 Basic simulation modules
 ========================
 
+This part covers the modules that synthesize time-ordered data: how the
+spacecraft scans the sky (pointing and quaternions), how an input sky is sampled
+into a TOD, beam convolution, the CMB dipole, and instrumental noise.
+
 .. toctree::
    :maxdepth: 1
 

docs/source/part4.rst

@@ -1,6 +1,11 @@
 Systematic effects
 ==================
 
+This part documents the instrumental systematic effects that can be injected
+into the TOD — gain drifts, detector non-linearity, half-wave-plate systematics,
+and pointing disturbances — together with the spin-based :math:`h`-maps that
+encode the orientation information of the scanning strategy.
+
 .. toctree::
    :maxdepth: 1
 

docs/source/part5.rst

@@ -1,6 +1,9 @@
 Data-reduction modules
 ======================
 
+This part describes how time-ordered data is reduced back into sky maps,
+covering the binned, destriped, and pair-differenced map-making algorithms.
+
 .. toctree::
    :maxdepth: 1
 

docs/source/part6.rst

@@ -1,6 +1,9 @@
 Tutorials
 =========
 
+This part collects worked, end-to-end examples that walk through a complete
+simulation, from setting up the instrument to producing maps.
+
 .. toctree::
    :maxdepth: 1
 

docs/source/refs.bib

@@ -159,4 +159,21 @@ @ARTICLE{2019:duivenvoorden:beamconv
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{2025:litebird:tomasi,
+       author = {{Tomasi}, M. and {LiteBIRD Collaboration}},
+        title = "{A simulation framework for the LiteBIRD instruments}",
+      journal = {JCAP},
+     keywords = {Astrophysics - Instrumentation and Methods for Astrophysics},
+         year = 2025,
+        month = nov,
+       volume = {2025},
+       number = {11},
+          eid = {040},
+        pages = {040},
+          doi = {10.1088/1475-7516/2025/11/040},
+archivePrefix = {arXiv},
+       eprint = {2507.04918},
+ primaryClass = {astro-ph.IM}
+}
+
 

docs/source/scanning.rst

@@ -848,7 +848,7 @@ API reference
     :undoc-members:
     :show-inheritance:
 
-.. automodule:: litebird_sim.pointings_in_obs
+.. automodule:: litebird_sim.observation_utilities
     :members:
     :undoc-members:
     :show-inheritance:

litebird_sim/__init__.py

@@ -146,6 +146,14 @@
     apply_quadratic_nonlin_to_observations,
 )
 from .observations import Observation, TodDescription
+from .observation_utilities import (
+    normalize_observations,
+    for_each_observation,
+    for_each_observation_with_pointings,
+    prepare_pointings,
+    precompute_pointings,
+    apply_hwp_to_obs,
+)
 from .pointing_sys import (
     get_detector_orientation,
     FocalplaneCoord,
@@ -155,11 +163,6 @@
 from .pointings import (
     PointingProvider,
 )
-from .pointings_in_obs import (
-    prepare_pointings,
-    precompute_pointings,
-    apply_hwp_to_obs,
-)
 from .profiler import TimeProfiler, profile_list_to_speedscope
 from .quaternions import (
     quat_rotation_x,
@@ -363,7 +366,10 @@
     "RNGHierarchy",
     # pointings.py
     "PointingProvider",
-    # pointings_in_obs.py
+    # observation_utilities.py
+    "normalize_observations",
+    "for_each_observation",
+    "for_each_observation_with_pointings",
     "prepare_pointings",
     "precompute_pointings",
     "apply_hwp_to_obs",

litebird_sim/beam_convolution.py

@@ -9,7 +9,7 @@
 from .hwp import HWP
 from .mueller_convolver import MuellerConvolver
 from .observations import Observation
-from .pointings_in_obs import (
+from .observation_utilities import (
     _get_hwp_angle,
     _normalize_observations_and_pointings,
     _get_pointings_array,

litebird_sim/dipole.py

@@ -8,6 +8,7 @@
 from .maps_and_harmonics import SphericalHarmonics
 from .observations import Observation
 from .spacecraft import SpacecraftPositionAndVelocity
+from .observation_utilities import for_each_observation_with_pointings
 
 # Updated imports to match the new constants.py structure
 from .constants import C_LIGHT_KM_OVER_S, H_OVER_K_B, T_CMB_K
@@ -882,43 +883,14 @@ def add_dipole_to_observations(
     # For convolved types we keep the full (θ, φ, ψ) columns; otherwise strip to (θ, φ).
     ptg_cols = slice(None) if apply_convolution else slice(0, 2)
 
-    if pointings is None:
-        if isinstance(observations, Observation):
-            obs_list = [observations]
-            if hasattr(observations, "pointing_matrix"):
-                ptg_list = [observations.pointing_matrix[:, :, ptg_cols]]
-            else:
-                ptg_list = [observations.get_pointings]
-        else:
-            obs_list = observations
-            ptg_list = []
-            for ob in observations:
-                if hasattr(ob, "pointing_matrix"):
-                    ptg_list.append(ob.pointing_matrix[:, :, ptg_cols])
-                else:
-                    ptg_list.append(ob.get_pointings)
-    else:
-        if isinstance(observations, Observation):
-            assert isinstance(pointings, np.ndarray), (
-                "You must pass a list of observations *and* a list "
-                + "of pointing matrices to add_dipole_to_observations"
-            )
-            obs_list = [observations]
-            ptg_list = [pointings[:, :, ptg_cols]]
-        else:
-            assert isinstance(pointings, list), (
-                "When you pass a list of observations to add_dipole_to_observations"
-                + ", you must do the same for `pointings`"
-            )
-            assert len(observations) == len(pointings), (
-                f"The list of observations has {len(observations)} elements, but "
-                + f"the list of pointings has {len(pointings)} elements"
-            )
-            obs_list = observations
-            ptg_list = [point[:, :, ptg_cols] for point in pointings]
-
-    for cur_obs, cur_ptg in zip(obs_list, ptg_list):
-        tod = getattr(cur_obs, component)
+    for cur_obs, tod, cur_ptg in for_each_observation_with_pointings(
+        observations, pointings, component
+    ):
+        # Callables (lazy pointings) are forwarded untouched, matching the
+        # behaviour of the underlying normalizer. Convolved calculations need
+        # the psi column too, hence ptg_cols rather than a hardcoded 0:2.
+        if isinstance(cur_ptg, np.ndarray):
+            cur_ptg = cur_ptg[:, :, ptg_cols]
 
         # Resolve the beam S-parameters when convolution is requested.
         # 1. Use the explicitly provided beam_alms if given.

litebird_sim/gaindrifts.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 from .observations import Observation
-from .seeding import regenerate_or_check_detector_generators
+from .observation_utilities import for_each_observation
 
 
 class GainDriftType(IntEnum):
@@ -548,23 +548,13 @@ def apply_gaindrift_to_observations(
     if drift_params is None:
         drift_params = GainDriftParams()
 
-    if isinstance(observations, Observation):
-        obs_list = [observations]
-    elif isinstance(observations, list):
-        obs_list = observations
-    else:
-        raise TypeError(
-            "The parameter `observations` must be an `Observation` or a list of `Observation`."
-        )
-    dets_random = regenerate_or_check_detector_generators(
-        observations=obs_list,
-        comm=None,
+    for cur_obs, tod, dets_random in for_each_observation(
+        observations,
+        component,
         user_seed=user_seed,
         dets_random=dets_random,
-    )
-
-    for cur_obs in obs_list:
-        tod = getattr(cur_obs, component)
+        requires_rng=True,
+    ):
         sampling_freq_hz = cur_obs.sampling_rate_hz
         focalplane_attr = getattr(cur_obs, drift_params.focalplane_group)
 

litebird_sim/hwp_diff_emiss.py

@@ -5,7 +5,8 @@
 
 from .hwp import HWP
 from .observations import Observation
-from .pointings_in_obs import _get_hwp_angle
+from .observation_utilities import _get_hwp_angle
+from .observation_utilities import for_each_observation
 
 
 # We calculate the additive signal coming from hwp harmonics.
@@ -83,13 +84,8 @@ def add_2f_to_observations(
             # in `observations.hwp_2f_tod`
             add_2f_to_observations(sim.observations, component="hwp_2f_tod")
     """
-    if isinstance(observations, Observation):
-        obs_list = [observations]
-    else:
-        obs_list = observations
-
     # iterate through each observation
-    for cur_obs in obs_list:
+    for cur_obs, tod, _ in for_each_observation(observations, component):
         if amplitude_2f_k is None:
             amplitude_2f_k = getattr(cur_obs, "amplitude_2f_k", None)
         assert amplitude_2f_k is not None, (
@@ -104,7 +100,7 @@ def add_2f_to_observations(
         hwp_angle = _get_hwp_angle(obs=cur_obs, hwp=hwp, pointing_dtype=pointings_dtype)
 
         add_2f(
-            tod=getattr(cur_obs, component),
+            tod=tod,
             hwp_angle=hwp_angle,
             pol_angle_rad=cur_obs.pol_angle_rad,
             amplitude_2f_k=amplitude_2f_k,

litebird_sim/hwp_harmonics/hwp_harmonics.py

@@ -17,7 +17,7 @@
 from ..input_sky import SkyInput
 from ..maps_and_harmonics import HealpixMap, SphericalHarmonics, interpolate_alm
 from ..observations import Observation
-from ..pointings_in_obs import (
+from ..observation_utilities import (
     _get_pointings_array,
 )
 from .jones_methods import (

litebird_sim/madam.py

@@ -15,7 +15,7 @@
 from .hwp import HWP
 from .mapmaking import ExternalDestriperParameters
 from .observations import Observation
-from .pointings_in_obs import _get_pointings_and_pol_angles_det
+from .observation_utilities import _get_pointings_and_pol_angles_det
 from .simulations import MpiDistributionDescr, Simulation
 
 

litebird_sim/mapmaking/binner.py

@@ -21,7 +21,7 @@
 from litebird_sim.coordinates import CoordinateSystem
 from litebird_sim.hwp import HWP
 from litebird_sim.observations import Observation
-from litebird_sim.pointings_in_obs import (
+from litebird_sim.observation_utilities import (
     _get_hwp_angle,
     _normalize_observations_and_pointings,
 )

litebird_sim/mapmaking/common.py

@@ -9,7 +9,7 @@
 
 from litebird_sim.coordinates import CoordinateSystem
 from litebird_sim.observations import Observation
-from litebird_sim.pointings_in_obs import _get_pointings_array, _get_pol_angle
+from litebird_sim.observation_utilities import _get_pointings_array, _get_pol_angle
 
 # The threshold on the conditioning number used to determine if a pixel
 # was really “seen” or not

litebird_sim/mapmaking/destriper.py

@@ -16,7 +16,7 @@
 from litebird_sim.hwp import HWP
 from litebird_sim.mpi import MPI_ENABLED, MPI_COMM_WORLD
 from litebird_sim.observations import Observation
-from litebird_sim.pointings_in_obs import (
+from litebird_sim.observation_utilities import (
     _get_hwp_angle,
     _normalize_observations_and_pointings,
 )