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 that is based on the implementation of the same 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:

  1. Linear gain drift

  2. Thermal gain drift

For any kind of 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(
    obs=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 specified 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 another enum class SamplingDist.

Following is an example of linear gain drift simulation parameters where the slope of gain for different detectors follow Gaussian distribution with mean 0.8 and standard deviation 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 show the thermal gain drift simulation parameters where the detector mismatch within a detector group has 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 to the TODs. The gaindrifts module provides method 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. The calibration period can be specified with the attribute GainDriftParams.calibration_period_sec. The following example shows the time evolution of the linear gain drift factor over four days with 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 linear gain drift factor that is to be multiplied with the sky TOD, to obtain the sky TOD with linear gain.

The module is written in a way to generate different gain slopes for different detectors. The slope (or the peak amplitude) of the linear gain is determined by the factor \(\sigma_{drift}\times\delta\), 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 same for all the detectors belonging to a given detector group. One may specify which detector parameter to be used to make detector group, using the attribute GainDriftParams.focalplane_group. It can be set to “wafer”, “pixtype” or even “channel”. For example, if GainDriftParams.focalplane_group = "wafer", all the detectors with same wafer name will be considered in one group and will have same noise timestream.

Once the noise timestreams are obtained for all focalplane groups, a mismatch for the detectors within a group is introduced by a random factor and the detector mismatch factor. The noise timestream with detector mismatch can be given as following

\[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 finally 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 same gain drift factor as they belong to the same focalplane group (grouped by wafer name). But when the detector mismatch is enabled, the two gain drift factors have 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
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:

GAUSSIAN = 1
UNIFORM = 0
litebird_sim.gaindrifts.apply_gaindrift_for_one_detector(det_tod: ndarray, sampling_freq_hz: float, det_name: str, drift_params: GainDriftParams | None = None, focalplane_attr: str | None = None, noise_timestream: ndarray | None = None, user_seed: int = 12345)

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 12345.

litebird_sim.gaindrifts.apply_gaindrift_to_observations(obs: Observation | List[Observation], drift_params: GainDriftParams | None = None, user_seed: int = 12345, component: str = 'tod')

The function to apply gain drift to the TOD of a Observation instance or a list of observations.

This function is a wrapper around apply_gaindrift_to_tod() that injects gain drift to the TOD object.

Parameters:

  • obs (Union[Observation, List[Observation]]) – An instance or a list of instances of Observation.

  • drift_params (GainDriftParams, optional) – The gain drift injection parameters object. Defaults to None.

  • user_seed (int, optional) – A seed provided by the user. Defaults to 12345.

  • component (str, optional) – The name of the TOD on which the gain drift has to be injected. Defaults to “tod”.

litebird_sim.gaindrifts.apply_gaindrift_to_tod(tod: ndarray, sampling_freq_hz: float, det_name: List | ndarray, drift_params: GainDriftParams | None = None, focalplane_attr: List | ndarray | None = None, user_seed: int = 12345)

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 (Union[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 (Union[List, np.ndarray], 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 list of the names of detector wafer. Defaults to None.

  • user_seed (int, optional) – A seed provided by the user. Defaults to 12345.