Add possibility to convert Non Linearity factor from 1/K to 1/SR

CHANGELOG.md

@@ -1,5 +1,7 @@
 # HEAD
 
+-   Add conversion from Temperature to Spectral Radiance units for the non linearity factor [#524](https://github.com/litebird/litebird_sim/pull/524)
+
 -   Add complete HWP Jones formalism, including band integration [#499](https://github.com/litebird/litebird_sim/pull/499)
 
 -   Fix missing parallelizations in multiple places [#521](https://github.com/litebird/litebird_sim/pull/521)

litebird_sim/hwp_harmonics/hwp_harmonics.py

@@ -1,21 +1,17 @@
+import logging
+import os
+
 import numpy as np
 import numpy.typing as npt
 from astropy import constants as const
 from astropy.cosmology import Planck18 as cosmo
 from ducc0.healpix import Healpix_Base
 from numba import njit
-import os
-import logging
 
 from litebird_sim.hwp_jones_parameters import HWPJonesParams
-from .jones_methods import (
-    compute_signal_for_one_detector as compute_signal_for_one_detector_jones,
-    integrate_inband_signal_for_one_detector as integrate_inband_signal_for_one_detector_jones,
-)
-from .mueller_methods import (
-    compute_signal_for_one_detector as compute_signal_for_one_detector_mueller,
-)
+
 from ..bandpass_template_module import bandpass_profile
+from ..constants import NUM_THREADS_ENVVAR
 from ..coordinates import CoordinateSystem
 from ..hwp_non_ideal import HWPFormalism, NonIdealHWP
 from ..input_sky import SkyGenerationParams
@@ -24,7 +20,15 @@
 from ..pointings_in_obs import (
     _get_pointings_array,
 )
-from ..constants import NUM_THREADS_ENVVAR
+from .jones_methods import (
+    compute_signal_for_one_detector as compute_signal_for_one_detector_jones,
+)
+from .jones_methods import (
+    integrate_inband_signal_for_one_detector as integrate_inband_signal_for_one_detector_jones,
+)
+from .mueller_methods import (
+    compute_signal_for_one_detector as compute_signal_for_one_detector_mueller,
+)
 
 COND_THRESHOLD = 1e10
 
@@ -34,11 +38,11 @@
 Tcmb0 = getattr(cosmo, "Tcmb0").value  # CMB temperature today in K
 
 
-def _dBodTrj(nu: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+def _dBodTrj(nu: npt.NDArray[np.float64] | np.float64):
     return 2 * k_B * nu * nu * 1e18 / c / c
 
 
-def _dBodTth(nu: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+def _dBodTth(nu: npt.NDArray[np.float64] | np.float64):
     x = h * nu * 1e9 / k_B / Tcmb0
     ex = np.exp(x)
     exm1 = ex - 1.0e0

litebird_sim/non_linearity.py

@@ -1,6 +1,8 @@
-import numpy as np
 from dataclasses import dataclass
 
+import numpy as np
+
+from .hwp_harmonics.hwp_harmonics import _dBodTth
 from .observations import Observation
 from .seeding import regenerate_or_check_detector_generators
 
@@ -57,8 +59,11 @@ def apply_quadratic_nonlin_for_one_sample(
 
 def apply_quadratic_nonlin_for_one_detector(
     tod_det,
+    det_bandcenter_ghz: np.float64 | None = None,
+    det_bandwidth_ghz: np.float64 | None = None,
     nl_params: NonLinParams | None = None,
     random: np.random.Generator | None = None,
+    conv_K_to_SR: bool = False,
 ):
     """This function applies the quadratic non-linearity on the TOD corresponding to
     only one detector.
@@ -81,26 +86,47 @@ def apply_quadratic_nonlin_for_one_detector(
 
         random (np.random.Generator, optional): A random number generator.
           Defaults to None.
+
+        conv_K_to_SR (bool, optional): Flag for temperature to spectral radiance
+            units conversion. Defaults to False.
+
+        det_freq_ghz (np.float64, optional): Detector central frequency in GHz.
+            Used in the K to SR conversion. Defaults to None.
+
+        det_bandwidth_ghz (np.float64, optional): Detector bandwidth in GHz.
+            Used in the K to SR conversion. Defaults to None.
+
     """
     assert random is not None, (
         "You should pass a random number generator which implements the `normal` method."
     )
     if nl_params is None:
         nl_params = NonLinParams()
 
-    g_one_over_k = random.normal(
+    g_nonlin = random.normal(
         loc=nl_params.sampling_gaussian_loc,
         scale=nl_params.sampling_gaussian_scale,
     )
 
+    if conv_K_to_SR:
+        assert det_bandcenter_ghz is not None and det_bandwidth_ghz is not None, (
+            "You should pass det_bandcenter_ghz and det_bandwidth_ghz when conv_K_to_SR is set to True."
+        )
+        conv_factor = _dBodTth(det_bandcenter_ghz)
+        # convert sampled g (1/K units) to (1/spectral radiance) units
+        g_nonlin = 1 / (conv_factor * (1 / g_nonlin) * det_bandwidth_ghz * 1e9)
+
     for i in range(len(tod_det)):
-        tod_det[i] += g_one_over_k * tod_det[i] ** 2
+        tod_det[i] += g_nonlin * tod_det[i] ** 2
 
 
 def apply_quadratic_nonlin(
     tod: np.ndarray,
+    bandcenter_ghz: np.ndarray,
+    bandwidth_ghz: np.ndarray,
     nl_params: NonLinParams | None = None,
     dets_random: list[np.random.Generator] | None = None,
+    conv_K_to_SR: bool = False,
 ):
     """Apply a quadratic nonlinearity to some time-ordered data
 
@@ -115,6 +141,9 @@ def apply_quadratic_nonlin(
             tod_det=tod[detector_idx],
             nl_params=nl_params,
             random=dets_random[detector_idx],
+            conv_K_to_SR=conv_K_to_SR,
+            det_bandcenter_ghz=bandcenter_ghz[detector_idx],
+            det_bandwidth_ghz=bandwidth_ghz[detector_idx],
         )
 
 
@@ -124,6 +153,7 @@ def apply_quadratic_nonlin_to_observations(
     component: str = "tod",
     user_seed: int | None = None,
     dets_random: list[np.random.Generator] | None = None,
+    conv_K_to_SR: bool = False,
 ):
     """
     Apply a quadratic nonlinearity to some time-ordered data
@@ -162,6 +192,8 @@ def apply_quadratic_nonlin_to_observations(
         List of per-detector random number generators. If not provided, and
         `user_seed` is given, generators are created internally. One of
         `user_seed` or `dets_random` must be provided.
+    conv_K_to_SR (bool, optional): Flag for temperature to spectral radiance
+        units conversion. Defaults to False.
 
     Raises
     ------
@@ -192,9 +224,14 @@ def apply_quadratic_nonlin_to_observations(
     # iterate through each observation
     for cur_obs in obs_list:
         tod = getattr(cur_obs, component)
+        bandcenter_ghz = getattr(cur_obs, "bandcenter_ghz")
+        bandwidth_ghz = getattr(cur_obs, "bandwidth_ghz")
 
         apply_quadratic_nonlin(
             tod=tod,
             nl_params=nl_params,
             dets_random=dets_random,
+            conv_K_to_SR=conv_K_to_SR,
+            bandcenter_ghz=bandcenter_ghz,
+            bandwidth_ghz=bandwidth_ghz,
         )

litebird_sim/simulations.py

@@ -2094,13 +2094,17 @@ def apply_quadratic_nonlin(
         component: str = "tod",
         append_to_report: bool = False,
         rng_hierarchy: RNGHierarchy | None = None,
+        conv_K_to_SR: bool = False,
     ):
         """A method to apply non-linearity to the observation.
 
         This is a wrapper around
         :func:`.apply_quadratic_nonlin_to_observations()` that
         applies non-linearity to a list of :class:`.Observation` instance. Random number generators are obtained from the detector-level layer. As default it uses
         the `dets_random` field of a :class:`.Simulation` object for this.
+
+        conv_K_to_SR (bool, optional) is a flag for temperature to spectral radiance
+        units conversion. Defaults to False.
         """
         if rng_hierarchy is None:
             rng_hierarchy = self.rng_hierarchy
@@ -2116,6 +2120,7 @@ def apply_quadratic_nonlin(
             user_seed=user_seed,
             component=component,
             dets_random=dets_random,
+            conv_K_to_SR=conv_K_to_SR,
         )
 
         if append_to_report and MPI_COMM_WORLD.rank == 0:
@@ -2132,6 +2137,7 @@ def apply_quadratic_nonlin(
             self.append_to_report(
                 markdown_template,
                 g=g,
+                conv_K_to_SR=conv_K_to_SR,
             )
 
     @_profile

litebird_sim/templates/report_quad_nonlin.md

@@ -1,5 +1,14 @@
-## Addition of quadratic non-lineariy 
+# Detector Non Linearity
 
-Quadratic non-linearity has been added to the TOD
+{% if g -%}
+- Quadratic Non Linearity was applied to the detectors.
+- Detector Non-Linearity params: {{ g }}
+{% else %}
+- Quadratic Non Linearity was NOT applied to the detectors.
+{% endif %}
 
-- detector non-linearity factor : {{g}}
+{% if conv_K_to_SR -%}
+- Non-Linearity g factor was converted from temperature to spectral radiance units. 
+{% else %}
+- Non-Linearity g factor was used in temperature units.
+{% endif %}