Gain drift injection#

Gain drift is the systematic that is multiplicative to time-ordered data. The LiteBIRD Simulation Framework provides a gain drift simulation module based on the same module in toast3. Though the exact nature of the gain drift depends on the specifics of the electronics, the gain drift module provides the functions to simulate two kinds of gain drifts:

Linear gain drifts;
Thermal gain drifts.

For any gain drift, one can use either the method of Simulation class Simulation.apply_gaindrift(), or any of the low-level functions: apply_gaindrift_to_observations(), apply_gaindrift_to_tod(), apply_gaindrift_for_one_detector(). The following example shows the typical usage of the method and low-level functions:

import numpy as np
import litebird_sim as lbs
from astropy.time import Time

start_time = Time("2034-05-02")
duration_s = 2*24*3600
sampling_freq_Hz = 1

# Creating a list of detectors. The detector name is used as one of
# two seeds to introduce unique and reproducible randomness to the
# gain drift for each detector.
dets = [
    lbs.DetectorInfo(
        name="det_A_wafer_1", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_1"
    ),
    lbs.DetectorInfo(
        name="det_B_wafer_1", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_1"
    ),
    lbs.DetectorInfo(
        name="det_C_wafer_2", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_2"
    ),
]

# Defining the gain drift simulation parameters
drift_params = lbs.GainDriftParams(
    drift_type=lbs.GainDriftType.LINEAR_GAIN,
    sampling_uniform_low=0.2,
    sampling_uniform_high=0.7,
)

sim1 = lbs.Simulation(
    base_path="linear_gd_example",
    start_time=start_time,
    duration_s=duration_s,
    random_seed=12345,
)

sim1.create_observations(
    detectors=dets,
    split_list_over_processes=False,
    num_of_obs_per_detector=1,
)

# Creating fiducial TODs
sim1.observations[0].gain_2_self = np.ones_like(sim1.observations[0].tod)
sim1.observations[0].gain_2_obs = np.ones_like(sim1.observations[0].tod)
sim1.observations[0].gain_2_tod = np.ones_like(sim1.observations[0].tod)
sim1.observations[0].gain_2_det = np.ones_like(sim1.observations[0].tod)

# Applying gain drift using the `Simulation` class method
sim1.apply_gaindrift(
    drift_params=drift_params,
    component="gain_2_self",
)

# Applying gain drift on the given TOD component of an `Observation` object
lbs.apply_gaindrift_to_observations(
    observations=sim1.observations,
    drift_params=drift_params,
    component="gain_2_obs",
)

# Applying gain drift on the TOD array. One must pass the name of the
# associated detectors (or just an array of string objects) to seed
# the reproducible randomness
lbs.apply_gaindrift_to_tod(
    tod=sim1.observations[0].gain_2_tod,
    sampling_freq_hz=sampling_freq_Hz,
    det_name=sim1.observations[0].name,
    drift_params=drift_params,
)

# Applying gain drift on the TOD arrays of the individual detectors.
# One must pass the name of the detector (or a string object) to seed
# the reproducible randomness.
for idx, tod in enumerate(sim1.observations[0].gain_2_det):
    lbs.apply_gaindrift_for_one_detector(
        det_tod=tod,
        sampling_freq_hz=sampling_freq_Hz,
        det_name=sim1.observations[0].name[idx],
        drift_params=drift_params,
    )

# The four TODs we obtain this way are equal to each other.

One has to specify the gain drift simulation parameters as an instance of the GainDriftParams class. The type of the gain drift can be selected using the enum class GainDriftType. The GainDriftParams class also offers the facility to specify the distribution of the slope for the linear gain and the distribution of the detector mismatch for the thermal gain, which can be specified with the help of the enum class SamplingDist.

Following is an example of linear gain drift simulation parameters where the slope of gain for different detectors follows Gaussian distribution with a mean of 0.8 and a standard deviation of 0.2:

import litebird_sim as lbs

drift_params = lbs.GainDriftParams(
    drift_type = lbs.GainDriftType.LINEAR_GAIN,
    sampling_dist = lbs.SamplingDist.GAUSSIAN,
    sampling_gaussian_loc = 0.8,
    sampling_gaussian_scale = 0.2,
)

The following example shows the thermal gain drift simulation parameters where the detector mismatch within a detector group has a uniform distribution varying between the factors 0.2 to 0.8:

import litebird_sim as lbs

drift_params = lbs.GainDriftParams(
    drift_type = lbs.GainDriftType.THERMAL_GAIN,
    sampling_dist = lbs.SamplingDist.UNIFORM,
    sampling_uniform_low = 0.2,
    sampling_uniform_high = 0.8,
)

Refer to the API reference for the full list of gain drift simulation parameters.

Linear gain drift#

Linear gain drift is the linearly increasing factor for TODs. The gaindrifts module provides methods and functions to simulate the linear gain drift with the possibility of periodic calibration. The calibration event resets the gain factor to one periodically after every calibration period interval. You can specify the calibration period with the attribute GainDriftParams.calibration_period_sec. The following example shows the time evolution of the linear gain drift factor over four days with a calibration period of 24 hours:

import numpy as np
import litebird_sim as lbs
from astropy.time import Time
import matplotlib.pyplot as plt

start_time = Time("2034-05-02")
duration_s = 4 * 24 * 3600
sampling_freq_Hz = 1

# Creating a list of detectors.
dets = [
    lbs.DetectorInfo(
        name="det_A_wafer_1", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_1"
    ),
    lbs.DetectorInfo(
        name="det_B_wafer_1", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_1"
    ),
]

# Defining the gain drift simulation parameters
drift_params = lbs.GainDriftParams(
    drift_type=lbs.GainDriftType.LINEAR_GAIN,
    calibration_period_sec=24 * 3600,
)

sim1 = lbs.Simulation(
    start_time=start_time,
    duration_s=duration_s,
    random_seed=12345,
)

sim1.create_observations(
    detectors=dets,
    split_list_over_processes=False,
    num_of_obs_per_detector=1,
)

sim1.observations[0].lingain_tod = np.ones_like(sim1.observations[0].tod)

# Applying gain drift using the `Simulation` class method
sim1.apply_gaindrift(
    drift_params=drift_params,
    component="lingain_tod",
)

plt.figure(figsize=(8, 5))

time_domain = (
    np.arange(sim1.observations[0].tod.shape[1]) / sampling_freq_Hz / 24 / 3600
)

for idx in range(sim1.observations[0].tod.shape[0]):
    plt.plot(
        time_domain,
        sim1.observations[0].lingain_tod[idx],
        label=sim1.observations[0].name[idx],
    )

plt.xlabel("Time (in days)")
plt.ylabel("Linear gain factor amplitude")
plt.legend()

(Source code)

Note that the figure above shows only the nature of the linear gain drift factor to be multiplied by the sky TOD.

The module generates different gain slopes for different detectors. The factor \(\sigma_{drift}\times\delta\) determines the slope (or the peak amplitude) of the linear gain, where \(\sigma_{drift}\) is a dimensionless parameter specified by GainDriftParams.sigma_drift and \(\delta\) is the random factor generated uniquely for each detector. The distribution of \(\delta\) or, conversely, the distribution of the gain slopes over all the detectors can be specified with attributes of SamplingDist enum class and the associated parameters listed in GainDriftParams.

Thermal gain drift#

The thermal gain drift is modelled as the gain drift due to \(1/f\) fluctuation in focalplane temperature. In the first step, the \(1/f\) noise timestream is generated from oversampled power spectral density given by

\[S(f) = \sigma_{drift}^2\left(\frac{f_{knee}}{f}\right)^{\alpha_{drift}}\]

The noise timestream is considered to be the same for all the detectors belonging to a given detector group. One may specify which detector parameter to use to make a detector group, using the attribute GainDriftParams.focalplane_group. Valid values are “wafer”, “pixtype”, or “channel”. For example, if GainDriftParams.focalplane_group = "wafer", all the detectors with the same wafer name will be considered in one group and have the same noise timestream.

Once the noise timestreams are obtained for all the groups in the focal plane, the code inserts a mismatch for the detectors within a group by a random factor and the detector mismatch factor. We can express the noise timestream with detector mismatch using the following expression:

\[t^{(mis)}_{stream} = (1 + \delta\times\alpha_{mis})t_{stream}\]

where \(\alpha_{mis}\) is the detector mismatch factor specified using the attribute GainDriftParams.detector_mismatch and \(\delta\) is the random factor generated uniquely for each detector. The distribution of \(\delta\) or, conversely, the distribution of noise timestream mismatch can be specified with attributes of SamplingDist enum class and the associated parameters listed in GainDriftParams.

The mismatched timestream is then scaled and passed through a responsivity function to obtain the thermal gain factor (\(\sigma\)):

\[\sigma = \text{responsivity_function}\left(1+\frac{ t^{(mis)}_{stream} \times \delta_T }{T_{bath}}\right)\]

where \(\delta_T\) is the amplitude of the thermal gain fluctuation in Kelvin unit, specified with attribute GainDriftParams.thermal_fluctuation_amplitude_K, and \(T_{bath}\) is the temperature of the focalplane in Kelvin unit specified with the attribute GainDriftParams.focalplane_Tbath_K.

The following example shows the comparison of thermal gain drift factor with or without detector mismatch over 100 seconds.

import numpy as np
import litebird_sim as lbs
from astropy.time import Time
import matplotlib.pyplot as plt

start_time = Time("2034-05-02")
duration_s = 100
sampling_freq_Hz = 1

# Creating a list of detectors. The three detectors belong to two
# different wafer groups.
dets = [
    lbs.DetectorInfo(
        name="det_A_wafer_1", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_1"
    ),
    lbs.DetectorInfo(
        name="det_B_wafer_1", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_1"
    ),
    lbs.DetectorInfo(
        name="det_C_wafer_2", sampling_rate_hz=sampling_freq_Hz, wafer="wafer_2"
    ),
]

# Defining the gain drift simulation parameters with no detector mismatch
drift_params_no_mismatch = lbs.GainDriftParams(
    drift_type=lbs.GainDriftType.THERMAL_GAIN,
    focalplane_group="wafer",
    detector_mismatch=0.0,
)

# Defining the gain drift simulation parameters with detector mismatch
drift_params_with_mismatch = lbs.GainDriftParams(
    drift_type=lbs.GainDriftType.THERMAL_GAIN,
    focalplane_group="wafer",
    detector_mismatch=1.0,
)

sim1 = lbs.Simulation(
    start_time=start_time,
    duration_s=duration_s,
    random_seed=12345,
)

sim1.create_observations(
    detectors=dets,
    split_list_over_processes=False,
    num_of_obs_per_detector=1,
)

sim1.observations[0].thermalgain_tod_no_mismatch = np.ones_like(
    sim1.observations[0].tod
)
sim1.observations[0].thermalgain_tod_with_mismatch = np.ones_like(
    sim1.observations[0].tod
)

# Generating the gain drift factors with no detector mismatch
sim1.apply_gaindrift(
    drift_params=drift_params_no_mismatch,
    component="thermalgain_tod_no_mismatch",
    user_seed=12345,
)

# Generating the gain drift factors with detector mismatch
sim1.apply_gaindrift(
    drift_params=drift_params_with_mismatch,
    component="thermalgain_tod_with_mismatch",
    user_seed=12345,
)

plt.figure(figsize=(8, 10))

plt.subplot(211)
for idx in range(sim1.observations[0].tod.shape[0]):
    plt.plot(
        sim1.observations[0].thermalgain_tod_no_mismatch[idx],
        label=sim1.observations[0].name[idx],
    )

plt.xlabel("Time (in seconds)")
plt.ylabel("Linear gain factor amplitude")
plt.title("Thermal gain drift factor with no detector mismatch")
plt.legend()

plt.subplot(212)
for idx in range(sim1.observations[0].tod.shape[0]):
    plt.plot(
        sim1.observations[0].thermalgain_tod_with_mismatch[idx],
        label=sim1.observations[0].name[idx],
    )

plt.xlabel("Time (in seconds)")
plt.ylabel("Linear gain factor amplitude")
plt.title("Thermal gain drift factor with detector mismatch")
plt.legend()
plt.tight_layout()

(Source code)

In the plots above, when there is no detector mismatch, det_A_wafer_1 and det_B_wafer_1 have the same gain drift factor as they belong to the same focal plane group (grouped by wafer name). However, when the detector mismatch is enabled, the two gain drift factors have the same shape due to the same underlying noise timestream but differ slightly in amplitude due to an additional random mismatch factor.

API reference#

class litebird_sim.gaindrifts.GainDriftParams(sampling_uniform_low: float = 0.0, sampling_uniform_high: float = 1.0, sampling_gaussian_loc: float = 0.7, sampling_gaussian_scale: float = 0.5, drift_type: GainDriftType = GainDriftType.LINEAR_GAIN, sigma_drift: float = 0.01, sampling_dist: SamplingDist = SamplingDist.UNIFORM, calibration_period_sec: int = 86400, focalplane_group: str = 'wafer', oversample: int = 2, fknee_drift_mHz: float = 20.0, alpha_drift: float = 1.0, detector_mismatch: float = 1.0, thermal_fluctuation_amplitude_K: float = 1.0, focalplane_Tbath_K: float = 0.1)#

Bases: object

A class to store the gain drift injection parameters.

The gain drift type can be one of the following:

Linear: It simulates gain drift that increases linearly in time. The gain factor resets to one periodically with time interval specified by the attribute GainDriftParams.calibration_period_sec.
Thermal: It simulates the gain drift as the fluctuation in the focalplane temperature. It offers the possibility to inject common mode drift to the TODs of detectors belonging to the same group of detectors identified by the attribute GainDriftParams.focalplane_group. This is enabled by setting the attribute GainDriftParams.detector_mismatch to 0.

The complete list of parameters is provided here:

Parameters common for the simulation of all types of gain drifts:
- drift_type (GainDriftType): Enumeration to determine the type of gain drift to be simulated. See GainDriftType.
- sigma_drift (float): A dimensionless parameter that determines the slope of gain drift in case of linear gain drift, and amplitude of thermal fluctuation in case of thermal gain drift.
- sampling_dist (SamplingDist): Enumeration to specify the distribution of the random scaling/mismatch factor applied on the gain drift. See SamplingDist.
Parameters that are specific to the simulation of linear gain drift:
- calibration_period_sec (int): This is the time period in seconds after which the linear gain drift resets periodically.
Parameters that are specific to the simulation of thermal gain drift:
- focalplane_group (str): Detector attribute to group the detectors. It is used to simulate same noise timestream for all the detectors belonging to a given group. It can be any of the detector attributes like “wafer”, “pixtype” or “channel”.
- oversample (int): The factor by which to oversample thermal noise FFT beyond the TOD size.
- fknee_drift_mHz (float): \(f_{knee}\) of the thermal drift power spectral density given in mHz.
- alpha_drift (float): The spectral index of thermal drift power spectral density.
- detector_mismatch (float): The factor that determines the degree of mismatch in thermal fluctuation of detectors belonging to same focalplane group. A value other than 0 implies no common gain. Whereas a value 0 sets the thermal gain to be same for all detectors in a focalplane group.
- thermal_fluctuation_amplitude_K (float): Amplitude of thermal gain fluctuation in Kelvin.
- focalplane_Tbath_K (float): Temperature of the focalplane in Kelvin.
Parameters for the sampling distributions:
- sampling_uniform_low (float): Lower boundary of the output for uniform distribution.
- sampling_uniform_high (float): Upper boundary of the output for uniform distribution.
- sampling_gaussian_loc (float): Mean of the Gaussian distribution.
- sampling_gaussian_scale (float): Standard deviation of the Gaussian distribution.

alpha_drift: float = 1.0#

calibration_period_sec: int = 86400#

detector_mismatch: float = 1.0#

drift_type: GainDriftType = 0#

fknee_drift_mHz: float = 20.0#

focalplane_Tbath_K: float = 0.1#

focalplane_group: str = 'wafer'#

oversample: int = 2#

sampling_dist: SamplingDist = 0#

sampling_gaussian_loc: float = 0.7#

sampling_gaussian_scale: float = 0.5#

sampling_uniform_high: float = 1.0#

sampling_uniform_low: float = 0.0#

sigma_drift: float = 0.01#

thermal_fluctuation_amplitude_K: float = 1.0#

class litebird_sim.gaindrifts.GainDriftType(value)#

Bases: IntEnum

An enumeration class to specify the type of gain drift injection.

The gain drift type can be:

LINEAR_GAIN: To inject linear gain drift in time with the possibility to calibrate the detectors at periodic interval
THERMAL_GAIN: To inject a gain drift with \(1/f\) psd mimicking the fluctuations in the focalplane temperature

LINEAR_GAIN = 0#

SLOW_GAIN = 2#

THERMAL_GAIN = 1#

class litebird_sim.gaindrifts.SamplingDist(value)#

Bases: IntEnum

An enumeration class to specify the distribution for the random scaling factor applied on the gain drift. For linear gain drift, it specifies the distribution of the slope of the gain drift. In case of thermal gain drift, it specifies the distribution of the detector mismatch.

The implemented distributions are:

UNIFORM: Uniform distribution. The lower and upper bound of the uniform distribution can be specified by the attributes GainDriftParams.sampling_uniform_low and GainDriftParams.sampling_uniform_high.
GAUSSIAN: Normal (Gaussian) distribution. The mean and standard deviation of the Gaussian distribution can be specified by the attributes GainDriftParams.sampling_gaussian_loc and GainDriftParams.sampling_gaussian_scale.

GAUSSIAN = 1#

UNIFORM = 0#

litebird_sim.gaindrifts.apply_gaindrift_for_one_detector(det_tod: ndarray, sampling_freq_hz: float, drift_params: GainDriftParams | None = None, noise_timestream: ndarray | None = None, random: Generator | None = None)#

This function applies the gain drift on the TOD corresponding to only one detector.

The linear drift is applied on the TODs in a periodic way with the period size specified in seconds by drift_params.callibration_period_sec. This is by assuming that the detectors are calibrated for linear gain drift periodically. The slope of the linear gain is determined randomly based on the detector name and the user-provided seed.

The thermal gain drift, on the other hand, is based on the fluctuation of the focalplane temperature modeled after \(1/f\) power spectral density (PSD). This \(1/f\) PSD is common to all the detectors belonging to the focalplane group identified by drift_params.focalplane_group. The function provides an option to introduce a mismatch between the individual detectors within a focalplane group with the parameter drift_params.detector_mismatch. This mismatch parameter along with a random number determines the extent of the mismatch of the thermal fluctuation within the focalplane group. Finally the thermal fluctuation is applied to the TODs according to the responsivity function of the detectors.

Parameters:

det_tod (np.ndarray) – The TOD array corresponding to only one detector.
sampling_freq_hz (float) – The sampling frequency of the detector in Hz.
det_name (str) – The name of the detector to which the TOD belongs. This name is used with user_seed to generate hash. This hash is used to set random slope in case of linear drift, and randomized detector mismatch in case of thermal gain drift.
drift_params (GainDriftParams, optional) – The gain drift injection parameters object. Defaults to None.
focalplane_attr (str, optional) – This is the parameter corresponding to the drift_params.focalplane_group attribute. For example, if drift_params.focalplane_group = 'wafer', the focalplane_attr will be the name of the detector wafer. Defaults to None.
noise_timestream (np.ndarray, optional) – The thermal noise time stream. Defaults to None.
user_seed (int, optional) – A seed provided by the user. Defaults to None.
random (np.random.Generator, optional) – A random number generator. Defaults to None.

litebird_sim.gaindrifts.apply_gaindrift_to_observations(observations: Observation | list[Observation], drift_params: GainDriftParams | None = None, user_seed: int | None = None, component: str = 'tod', dets_random: list[Generator] | None = None)#

Apply gain drift to one or more observations.

This function injects gain drift into the time-ordered data (TOD) of one or more Observation instances. It wraps apply_gaindrift_to_tod(), and ensures proper setup of per-detector random number generators using either a user-provided seed or a list of pre-initialized RNGs.

Parameters:

observations (Observation or list of Observation) – A single Observation instance or a list of them.
drift_params (GainDriftParams, optional) – Parameters defining the gain drift injection (e.g., linear or thermal). If not provided, a default configuration is used.
user_seed (int, optional) – Base seed to build the RNG hierarchy and generate detector-level RNGs that overwrite any eventual dets_random. Required if dets_random is not provided.
component (str, optional) – Name of the TOD attribute to modify. Defaults to “tod”.
dets_random (list of np.random.Generator, optional) – 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.

Raises:

TypeError – If observations is neither an Observation nor a list of them.
ValueError – If neither user_seed nor dets_random is provided.
AssertionError – If the number of random generators does not match the number of detectors.

litebird_sim.gaindrifts.apply_gaindrift_to_tod(tod: ndarray, sampling_freq_hz: float, drift_params: GainDriftParams | None = None, focalplane_attr: list | ndarray | None = None, dets_random: list[Generator] | None = None)#

The function to apply the gain drift to all the detectors of a given TOD object.

This function is a wrapper around apply_gaindrift_for_one_detector() that applies the gain drift on each detector TODs of the TOD object. In case of thermal gain drift injection, this function computes the thermal noise fluctuations at once for all the detectors belonging to the focalplane group specified by drift_params.focalplane_group and passes them to apply_gaindrift_for_one_detector() with individual TOD arrays to inject thermal gain drift.

Parameters:

tod (np.ndarray) – The TOD object consisting TOD arrays for multiple detectors.
sampling_freq_hz (float) – The sampling frequency of the detector in Hz.
det_name (list | np.ndarray) – The list of the name of the detectors to which the TOD arrays correspond. The detector names are used to generate unique and reproducible random numbers for each detector.
drift_params (GainDriftParams, optional) – The gain drift injection parameters object. Defaults to None.
focalplane_attr (list | np.ndarray | None) – This is the parameter corresponding to the drift_params.focalplane_group attribute. For example, if drift_params.focalplane_group = 'wafer', the focalplane_attr will be the list of the names of detector wafer. Defaults to None.
user_seed (int, optional) – A seed provided by the user. Defaults to None.
dets_random – list of np.random.Generator List of per-detector random number generators. Must match the number of detectors. Typically obtained from an RNGHierarchy.