Differences

This shows you the differences between two versions of the page.

--- calibration_simulations [2018/04/17 23:57] – created mauriziotomasi
+++ calibration_simulations [2019/02/20 04:56] (current) – Put an "outdated" warning message at the top of the page. mauriziotomasi
@@ Line 1: / Line 1: @@
+**This page contains outdated information**. Refer to the report «Simulating dipole calibration for PICO», by Maurizio Tomasi, for the most up-to-date information {{::tomasi-calibration-note.pdf|}}.
 ====== Calibration simulations ======
-In this page we collect information and results about the simulation of the photometric calibration procedure.
+In this page we collect information and results about the simulation
+of the photometric calibration procedure.
 ===== Hypothesis =====
-We assume to calibrate the instrument using the solar dipole of the CMB, whose peak-to-peak amplitude is 7 mK. We use [[https://github.com/hpc4cmb/toast|TOAST]] to produce the timelines, and [[https://github.com/ziotom78/dacapo_calibration|DaCapo]] to compute the calibration factors.
+I assume to calibrate the instrument using the solar dipole of the
+CMB, whose peak-to-peak amplitude is ~7 mK. I use
+[[https://github.com/hpc4cmb/toast|TOAST]] to produce the timelines,
+and [[https://github.com/ziotom78/dacapo_calibration|DaCapo]] to
+compute the calibration factors. The scripts used to run the
+simulation have been prepared by Andrea Zonca and are available at
+[[https://github.com/zonca/pico-simulations|pico-simulations]].
+As a reference, the dipole in Ecliptic coordinates has the following
+shape:
+{{ ::dipole-ecliptic.png?600 |}}
 ===== Observation of the dipole =====
-I used the scanning strategy parameters listed at the page [[optimizingscanstrategy:start|Optimizing Scan Strategy]], with **no** kinematic dipole and the Galactic signal simulated by PySM for a W-band detector, which should be the best case in terms of Galactic contamination.
+I used the scanning strategy parameters listed at the page
+[[optimizingscanstrategy:start|Optimizing Scan Strategy]], with **no**
+kinematic dipole and the Galactic signal simulated by PySM for a
+W-band detector, which should be the best case in terms of Galactic
+contamination.
-The following plots show a comparison between the proposed scanning strategy for PICO and the strategy used for Planck. This is the sky coverage as a function of time:
+The sky coverage as a function of time shows a sharp rise in the first
+few hours, reaching ~50% of the sky in one precession period. Then,
+the value increases slowly till it reaches 100% after nearly 6
+months. (For comparison, fsky varies more linearly in the case of
+Planck.) The following figure shows fsky as a function of time:
-{{ ::pico-fsky.png |}}
+{{ ::pico-fsky-6months.png?600 |}}
-The precession allows PICO to observe wider portions of the sky in the same time. This is good for dipole calibration as well, as it allows the detectors to potentially sample larger variations of the dipole, thus increasing the S/N of the dipole measurement.
+The precession allows PICO to observe wider portions of the sky in the
+same time. This is good for dipole calibration as well, as it allows
+the detectors to potentially sample larger variations of the dipole,
+thus increasing the S/N of the dipole measurement.
-{{ ::pico-dipole-visibility.png |}}
+A simple parameter that encapsulates much of the details about
+photometric calibration is the peak-to-peak variation in the dipole
+signal during some fixed time span. As PICO detectors are not meant to
+perform absolute measurements, the //modulation// of the dipole signal
+consititutes our calibration source. The peak-to-peak amplitude of the
+dipole is therefore a good tracer of the S/N for the calibration
+signal.
-The plot above shows the peak-to-peak variation of the dipole signal over periods of 6 hours. In general, high values (more than 5 mK) are good, lower values (less than 3) are bad, because the signal is not strong enough to allow for a good calibration. Planck followed a very regular scanning strategy, which caused the dipole to be poorly visible for several weeks (this effect was dramatic during the 2nd and 4th sky surveys). The PICO boresight direction often samples the equator of the dipolar signal (which is bad), but such instances are usually short (a few days, and this is good).
+Consider this image, which shows the simulated timeline produced by
+TOAST for one of the W-band detectors (temperature in K as a function
+of time in seconds, no Galaxy is present):
-{{ ::pico-dipole-histogram.png |}}
+{{ :pico-nofg-totaltod.png?600 |}}
+The low-frequency modulation is due to the dipole signal, while the
+high-frequency noise is detector noise. As the plot shows ~100 seconds
+of data, the peak-to-peak variation of the dipole signal within this
+time frame is less than 3 mK. The following image shows how the
+peak-to-peak dipole amplitude changes as a function of the time window
+used to estimate the peak-to-peak value:
+{{ ::dipole-amplitude.gif?600 |}}
+(Go to https://zzz.physics.umn.edu/ipsig/_media/dipole-amplitude.gif to see the animation.)
+There are several interesting facts that become apparent from the
+image above. First of all, there are high-frequency fluctuations which
+disappear for longer time windows. Specifically, as the length of the
+time window approaches the precession period, high-frequency
+fluctuations are present no longer. A second interesting fact is that
+the plot shows low-frequency fluctuations as well. These are
+unavoidable and are due to the yearly variation in the angle between
+the spacecraft spin axis and the CMB dipole axis. As a matter of fact,
+these fluctuations are present even for much longer time windows:
+{{ ::dipole-amplitude-overview.png?600 |}}
-This plot is a histogram of the dipole amplitudes shown in the previous plot.
 ===== Running DaCapo =====
+I assume that the calibration code to be used in a real mission would
+include a component-separation step able to reduce the bias due to the
+contamination of the dipole signal by Galactic large-scale
+structures. Therefore, in our simulations I assume that the only
+signal in the sky apart from the dipole is the CMB itself. From the
+point of view of calibration, CMB acts as a noise signal which is
+correlated among detectors.
+The DaCapo algorithm is explained in
+[[https://www.aanda.org/articles/aa/abs/2016/10/aa26632-15/aa26632-15.html|Planck
+results. V. LFI calibration]] and in the source code of the
+Python version of
+[[https://github.com/ziotom78/dacapo_calibration/raw/master/dacapo_calibration.pdf|DaCapo]],
+which has been written with the aim to be readable as an article. The
+algorithm works as follows:
+  - Apply an extended version of the destriping equation, which reconstructs the baselines of 1/f noise as well as the gain factors. This uses the following assumptions:
+    * 1/f baselines have zero mean
+    * Gains are constrained by a model of the dipole, which must be provided as input
+  - The equation is solved iteratively using the Conjugate Gradient algorithm. At the end, the following estimates are available:
+    * 1/f baselines
+    * Timeline of gain factors
+    * Sky map (temperature only) without the dipole
+  - Iterate the algorithm again, subtracting the sky map obtained in the previous step from the timelines, in order to remove any systematic due to spurious signals (CMB, foregrounds).
+  - When the 1/f baselines, the gain factors, and the sky map no longer change, save all the results and quit.
+  - DaCapo must be run once per each detector, so that the sky map produced as output is typically noisier than the maps expected from a typical configuration of detectors.
+I employed the DaCapo algorithm in the following pipeline:
+  - Use TOAST to produce timelines using two detectors (0A and 0B) along the boresight direction, the nominal scanning strategy and realistic noise, and save them to FITS files. The gain used in this step is 1.0 (constant) for both detectors. TOAST produces a sky map (including the dipole) using the Madam map-maker (a destriper) and data from both detectors 0A and 0B.
+  - Run DaCapo on the FITS files produced during the previous stage. As the algorithm works on single detectors, the code must be run twice (once for detector 0A and once for 0B).
+  - If DaCapo converges, re-run TOAST as in the previous step, but this time use the gains calculated by DaCapo instead of assuming a constant (1.0) gain.
+  - At the end of the simulation, take the difference between the maps produced by TOAST in steps 1. and 3. and compute the power spectrum. These spectra represent the error caused by an imperfect estimate of gain drifts during the mission.
+Our simulations assume that each detector has a calibration factor
+equal to 1.0. The calibration code is therefore expected to retrieve
+the same value from the simulated data, plus some statistical
+noise. Of course, the code does not know that the gain is constant, so
+it always tries to estimate the dependency of gain on time.
+I have run a few simulations using DaCapo on the simulated timelines
+produced by Andrea Zonca. Here are the details:
+  * 1 year of observation
+  * 2 boresight W-band detectors
+  * No Galaxy
+  * Planck2015 solar dipole, including the quadrupolar correction
+  * No orbital dipole
+  * 1/f + WN, with fknee = 10 mHz (corresponding to a period of 100 s).
+  * Constant gain (1.0)
+As DaCapo implements a destriper as well, let's check that the 1/f is
+correctly estimated. The baseline I used was 10 s:
+{{ :dacapo-cleaned-tod-dsp.png?600 |}}
+(“Original” is the TOD produced by TOAST and fed to DaCapo, while
+“Cleaned” is the same TOD minus the 1/f baselines calculated by
+DaCapo.) This plot is complicate to understand due to the presence of
+the sky signal (the dipole peak and its harmonics). Let's compute the
+PSD of the noise component alone: this is easy, since TOAST saves the
+TOD of the foreground signal in a separate column in the FITS files.
+{{ ::dacapo-cleaned-tod-dsp-nofg.png?600 |}}
+The peak at 167 mHz corresponds to the spin period (60 s): note that
+the height of this residual peak matches the level of the
+high-frequency WN tail. Apart from that, 1/f is correctly suppressed
+by the code.
+The estimated gains for a time window of 10 hours (the same as the
+precession period) are shown here:
+{{ ::dacapo-gains-nofg-10h.png?600 |}}
+You can see that the RMS of the gain factors varies with time. This is
+anticorrelated with the peak-to-peak amplitude of the dipole (see the
+plot above). The average value of the gains shows some residual bias,
+as it is not exactly 1.0: this bias is due to the fact that TOAST
+convolves the sky and the dipole signals with a beam of some finite
+size, but DaCapo considers a pencil beam (i.e., a Dirac delta). The
+usual way to deal with this is to use a 4π convolver, like the one
+used in Appendix A of
+[[https://www.aanda.org/articles/aa/abs/2016/10/aa26632-15/aa26632-15.html|Planck
+results. V. LFI calibration]], but for our purpose it's easier to
+subtract an offset so that the average value of the gain is 1.0. The
+RMS of the timeline is 0.04%, which can be considered as an estimate
+of the relative calibration error.
+===== Output maps =====
+As sketched in the previous section, I produced full-sky maps with
+TOAST, calibrate them with DaCapo and re-run TOAST with the gains
+estimated by DaCapo to get a new full-sky map which contains the
+effect of imperfect gain drift reconstruction. The difference between
+the two maps show the projected effect of gain drifts on the sky.
+I have run two sets of simulations, using the following assumptions:
+  - Two noiseless receivers (0A and 0B);
+  - Two W-band receivers with nominal white noise and 1/f noise with knee frequency 10 mHz.
+The differenced maps produced in the two cases are quite different:
+{{ :systematicswg:cal_error_maps.png?600 |}}
-==== No Galactic mask ====
+Of course, the scale of the effect considers the presence of two
+detectors on the whole focal plane, which is not the case for PICO.
-The following plot shows the results of running DaCapo on the 2-year mission, assuming to recompute the gain constant every 6 hours. (This is the same timescale used to produce the plots above, in the section "Observation of the dipole").
+Computing the BB spectrum of the map in the two cases and confronting
+it with the input map shows the scale of the effect (as well as two
+power spectra assuming only tensor modes, provided as a reference):
-{{ :pico-gain-timestream.png |}}
+{{ :systematicswg:calibration_spectrum_bb.png?600 |}}
-The true gain should be 1.000, but their average 0.994, with a standard deviation of each sample equal to 0.010. This means that the accuracy of the calibration is 0.6% and the error in the gain estimate is 1.0%. We could probably do better if we increase the time span used to compute the gain constant (e.g., from 6 to 24 or 48 hours). In order to do so, we need to understand how much stable we expect the PICO detectors to be.
+In order to produce this plot, I rescaled the spectrum by 2/N, where
+N is the number of W-band bolometers in the focal plane of PICO, and
+then I performed a NET-weighted average over all the cosmology bands
+from 60 GHz (band 7) to 220 GHz (band 22): this
+assumes that the gain fluctuations are negligibly correlated among
+detectors. This is the case for the case with realistic noise, but not
+for the noiseless simulation; this should be expected, as the only
+noise here is caused by the CMB, which is the same for detectors 0A
+and 0B:
-==== Using a Galactic mask ====
+{{ :systematicswg:gain_timeline.png?600 |}}
-The use of a Galactic mask should help to reduce the bias shown above (i.e., gain average is 0.994 instead of 1.0), at the price of increasing the noise (because less data are used in the calculation). The code is still running on NERSC, once results will be available I'll post them here. [This was written on Wednesday, April 4th 2018].

User Tools

Differences

Page Tools