Difference between revisions of "TOI processing LFI"

From Planck PLA 2015 Wiki
Jump to: navigation, search
(4 K total-power and Fix)
(Photometric Calibration)
 
(72 intermediate revisions by 4 users not shown)
Line 2: Line 2:
 
==Overview==
 
==Overview==
  
The LFI Level2 Pipeline analyzes each horn of the instrument separately, one pointing period at time, and store results in object the length of an OD. Each diode of the horn is corrected from systematic, differentiated and then combined with its complementary diode in the same radiometer. The horn is then calibrated and the photometric calibration is applied.
+
The LFI Level2 Pipeline analyzes data from each horn of the instrument separately, one pointing period at time, and stores the results in an object the length of an OD. Each diode of the horn is corrected for systematic effects.  Next, measurements of the sky and the 4K load are differenced, then the signals from one diode are combined with signals from the complementary diode in the same radiometer. Finally, photometric calibration is applied for each horn.
  
 
===Pre-processing===
 
===Pre-processing===
  
Before the run of the Level2 pipeline and to improve the analysis the Mission information and data sampling divisions are stored in the database.
+
Before the Level 2 pipeline is run, the Mission information and data sampling divisions are stored in the database, in order to improve the analysis.
  
The Mission information is a set of objects, one for each Operational Day (OD, as defined in the [[glossary]]), in which are stored Pointing Period data: DPC pointing ID (where 1 is the first pointing of the nominal mission), PSO pointing ID, start OBT of the pointing maneuver, start OBT of the stable pointing, end OBT of the pointing, spin axis ecliptic longitude and latitude.
+
The Mission information is a set of objects, one for each Operational Day (OD, as defined in the [[glossary]]), in which are stored pointing period data: the DPC pointing ID (where 1 is the first pointing of the nominal mission), PSO pointing ID, start OBT of the pointing maneuver, start OBT of the stable pointing, OBT of the end of the stable pointing, and spin axis ecliptic longitude and latitude.
  
The sampling information is a set of objects, one for each LFI frequency, in which are stored for each pointing ID: start OBT of the pointing maneuver, start OBT of the stable pointing, end OBT of the pointing, number of samples of the pointing, number of stable samples of the pointing, start sample of the stable pointing and sample number from the start of the nominal mission. Valid samples and OBTs are defined where any of the radiometers from that frequency cohort contain valid data.
+
The sampling information is a set of objects, one for each LFI frequency, in which are stored for each pointing ID: start OBT of the pointing maneuver, start OBT of the stable pointing, end OBT of the pointing, number of samples of the pointing, number of stable samples of the pointing, start sample of the stable pointing and sample number from the start of the nominal mission. Valid samples and OBTs are defined to be those where any of the radiometers in a frequency band contain valid data.
  
 
==ADC Correction==
 
==ADC Correction==
  
During analysis it appeared that white noise and calibration seemed affected by something in common. It turn out to be a non linearity in the Analogic/Digital Converter on board. More on {{PlanckPapers|planck2013-p02}}{{PlanckPapers|planck2013-p02a}}.
+
During the analysis of LFI date, we discovered a systematic common to both white noise and calibration. It turn out to be a non linearity in the Analogue to Digital Converter (ADC) on board. More detail is supplied in {{PlanckPapers|planck2013-p02}}{{PlanckPapers|planck2013-p02a}}, {{PlanckPapers|planck2014-a03||Planck-2015-A03}} and {{PlanckPapers|planck2014-a04||Planck-2015-A04}}.
  
 
===Evaluation===
 
===Evaluation===
Line 24: Line 24:
 
where <math> V </math> is the voltage input, <math> \gamma </math> is the DAE gain, <math> \Delta </math> is the DAE offset and <math> x_0 </math> is the DAE <math> T_{zero} </math>.
 
where <math> V </math> is the voltage input, <math> \gamma </math> is the DAE gain, <math> \Delta </math> is the DAE offset and <math> x_0 </math> is the DAE <math> T_{zero} </math>.
  
We can model the non-linearity as a function of the input voltage <math> R(V) </math>. So we have the apparent inferred voltage <math>V^'</math> and we can link it to the actual input voltage with:
+
We can model the non-linearity as a function of the input voltage <math> R(V) </math>. So we have the apparent inferred voltage <math>V'</math> and we can link it to the actual input voltage with:
  
:<math> ((V - \Delta) \gamma + x_0) R(V) = X = ((V^' - \Delta) \gamma + x_0 </math>
+
:<math> ((V - \Delta) \gamma + x_0) R(V) = X = ((V' - \Delta) \gamma + x_0 </math>
  
 
so that:
 
so that:
  
:<math> R(V) = (V^' - \Delta) \gamma + x_0 \over (V - \Delta) \gamma + x_0 </math>
+
:<math> R(V) = (V' - \Delta) \gamma + x_0 \over (V - \Delta) \gamma + x_0 </math>
  
 
Since <math> V >> \Delta </math> and <math> V \gamma >> x_0 </math> we can use the much simpler relation:
 
Since <math> V >> \Delta </math> and <math> V \gamma >> x_0 </math> we can use the much simpler relation:
:<math> R(V) = {V^' \over V} </math>
+
:<math> R(V) = \frac{V'}{V} </math>
  
 
and we expect it to be very near to unity for all <math> V </math>.
 
and we expect it to be very near to unity for all <math> V </math>.
  
To find the response curve we have only the apparent voltage to work with, so we had to use the inverse response function <math> R^'(V^') </math> and replace the real input voltage with <math> T_{sys} </math> times the time varying gain factor <math> G(t) </math>.
+
To find the response curve we have only the apparent voltage to work with, so we had to use the inverse response function <math> R'(V') </math> and replace the real input voltage with <math> T_{sys} </math> times the time varying gain factor <math> G(t) </math>.
  
:<math> V = V^'R^'(V^') = G(t)T_{sys} </math>
+
:<math> V = V'R'(V') = G(t)T_{sys} </math>
  
 
If we introduce a small signal on top of <math> T_{sys} </math> which leads to increased detected voltage and corresponding apparent voltage increment:
 
If we introduce a small signal on top of <math> T_{sys} </math> which leads to increased detected voltage and corresponding apparent voltage increment:
  
:<math> V + \delta V = (V^' + \delta V^') = (V^' + \delta V^') R^' (V^' + \delta V^') = G(t) (T_{sys} + \delta T) </math>
+
:<math> V + \delta V = (V' + \delta V') R' (V' + \delta V') = G(t) (T_{sys} + \delta T) </math>
  
so carrying out the differentiation respect to <math> V^' </math> to the relation between true and apparent signal voltage leads to:
+
so carrying out the differentiation with respect to <math> V' </math> to the relation between true and apparent signal voltage leads to:
  
:<math> \delta V = \left( V^' {dR^'(V^') \over dV^' + R^'(V^')} \right) \delta V^' = G(t) \delta T </math>
+
:<math> \delta V = \left( V' {dR'(V') \over dV' + R'(V')} \right) \delta V' = G(t) \delta T </math>
  
 
We now assume <math> T_{sys} </math> and <math> \delta T </math> are fixed and that the variations are due to slow drifts in the gain. So we can isolate the terms:
 
We now assume <math> T_{sys} </math> and <math> \delta T </math> are fixed and that the variations are due to slow drifts in the gain. So we can isolate the terms:
  
:<math> V^' = {G(t) T_{sys} \over R^'(V^')} </math>
+
:<math> V' = {G(t) T_{sys} \over R'(V')} </math>
:<math> \delta V^' = {G(t) \delta T \over V^' {dR^'(V^') \over dV} + R^'(V^')} </math>
+
:<math> \delta V' = {G(t) \delta T \over V' {dR'(V') \over dV} + R'(V')} </math>
  
 
Combining the equations through the gain factor to remove it:
 
Combining the equations through the gain factor to remove it:
  
:<math> {V^' R^'(V^') \over T_{sys}} = {\delta V^' {dR^'(V^') \over dV^'} + R^'(V^') \over \delta T } </math>
+
:<math> {V' R'(V') \over T_{sys}} = {\delta V' {dR'(V') \over dV'} + R'(V') \over \delta T } </math>
  
 
Rearranging and putting <math> a = {\delta T \over T_{sys}} </math>
 
Rearranging and putting <math> a = {\delta T \over T_{sys}} </math>
  
:<math> {dR^'(V^') \over dV^' } = \left( {a \over dV^'} - {1 \over V^'} \right) R^'(V^') </math>
+
:<math> {dR'(V') \over dV' } = \left( {a \over dV'} - {1 \over V'} \right) R'(V') </math>
  
So there is the expected direct proportionality of <math> \delta V^' </math> to <math> V^' </math> due to the assumption that the variations in voltage are due to overall gain drift, so the amplitude of voltage and signal will vary together. Then there is the additional differential term which will pull signal amplitude away from the linear relationship. So if we plot measured white noise or dipole gain factor against recovered voltage we should see this linear curve with variations due to local slope changes at particular voltages. The linear part can be taken out and the differential part fitted. This was numerically integrated up to get the inverse response curve, what we need to convert the measured voltages to corrected voltages.
+
So there is the expected direct proportionality of <math> \delta V' </math> to <math> V' </math> based on the assumption that the variations in voltage are due to overall gain drift, so the amplitude of voltage and signal will vary together. Then there is the additional differential term which will pull the signal amplitude away from the linear relationship. So if we plot measured white noise or dipole gain factor against recovered voltage, we should see a linear curve with variations due to local slope changes at particular voltages. The linear part can be taken out and the differential part fitted. This was numerically integrated up to get the inverse response curve, which we then used to convert the measured voltages to corrected voltages.
  
 
===Application===
 
===Application===
  
For each of the 44 LFI diodes there is the corresponding object in the Database. Each object contains 4 columns: the input voltages coming from the sky channel and the corresponding linearized output, the input voltages coming from the reference channel and the corresponding linearized output.  
+
For each of the 44 LFI diodes there is the corresponding corrected object in the Database. Each object contains 4 columns: the input voltages coming from the sky channel and the corresponding linearized output, the input voltages coming from the 4K reference channel and the corresponding linearized output.  
  
 
Data loaded by the module are used to initialize two different interpolators using CSPLINE and the functions from gsl ([http://www.gnu.org/software/gsl/ GNU Scientific Libraries]) libraries. The interpolators are then used to correct each sample.
 
Data loaded by the module are used to initialize two different interpolators using CSPLINE and the functions from gsl ([http://www.gnu.org/software/gsl/ GNU Scientific Libraries]) libraries. The interpolators are then used to correct each sample.
Line 72: Line 72:
 
==Spikes Removal==
 
==Spikes Removal==
  
Some of the LFI receivers exhibit a small artifact with exactly 1 second repetition, which visible in the power spectra. The effect is a set of spikes at 1 Hz and harmonics. The spurious signal is very well modeled and is removed from the timelines. More information can be found in {{PlanckPapers|planck2013-p02}} {{PlanckPapers|planck2013-p02a}}.
+
Some of the LFI receivers exhibit a small artifact with exactly 1 second period, that produces effects visible in the power spectra. The effect is a set of spikes at 1 Hz and harmonics. The spurious signal is very well modeled and is removed from the timelines. More information can be found in {{PlanckPapers|planck2013-p02}} {{PlanckPapers|planck2013-p02a}} {{PlanckPapers|planck2014-a03||Planck-2015-A03}} {{PlanckPapers|planck2014-a04||Planck-2015-A04}}.
  
 
===Modeling===
 
===Modeling===
  
The cause of the spikes at 1 Hz and harmonics is a tiny 1 second square wave embedded in affected channels. The method to estimate the 1 Hz signal is to build a template in time domain synchronized with the spurious signal. The first step is dividing each second of data into time bins using OBT. The number of bins is computed using:
+
The cause of the spikes at 1 Hz and harmonics is a tiny 1 second square wave embedded in signals from the affected channels. The method to estimate the 1 Hz signal is to build a template in the time domain synchronized with the spurious signal. The first step is dividing each second of data into time bins using OBT. The number of bins is computed using:
  
 
:<math> nbins = fsamp * template\_resolution</math>
 
:<math> nbins = fsamp * template\_resolution</math>
  
where fsamp is the sampling frequency and is 136 at 70 GHz, 80 at 44 GHz and 56 at 30 GHz. Then the bins vector is initialized with time intervals. To avoid aliasing effects template resolution is <math> \sqrt {3} </math>
+
where fsamp is the sampling frequency and is 136 Hz at 70 GHz, 80 at 44 GHz and 56 at 30 GHz. Then the bins vector is initialized with time intervals. To avoid aliasing effects the template resolution is <math> \sqrt {3} </math>
 
.
 
.
We can write the process adding an index to the time sample: lower index denotes the particular time sample, while the upper index labels the bin into which the sample falls. The linear filter can be written as:
+
We can write the process adding an indices to the time sample: the lower index denotes the particular time sample, while the upper index labels the bin into which the sample falls. The linear filter can be written as:
  
 
:<math> s(t_{i}^{j}) = a_j  \left(1- \Delta x (t_{i}^{j}) \right) + a_{j+1} \Delta x (t_{i}^{j})</math>
 
:<math> s(t_{i}^{j}) = a_j  \left(1- \Delta x (t_{i}^{j}) \right) + a_{j+1} \Delta x (t_{i}^{j})</math>
  
Here <math> \Delta x (t_{i}^{j})</math> is the filter weight which is determined by where within the bin sample lies. If we use <math> t^j </math> with only an upper index to denote the start of each bin, then we can write the filter weight as follows:
+
Here <math> \Delta x (t_{i}^{j})</math> is the filter weight which is determined by where within the bin the sample lies. If we use <math> t^j </math> with only an upper index to denote the start of each bin, then we can write the filter weight as follows:
  
 
:<math>  \Delta x (t_{i}^{j}) = {{{t_i^j - t^j} \over {t^{j+1} - t^j}}} </math>
 
:<math>  \Delta x (t_{i}^{j}) = {{{t_i^j - t^j} \over {t^{j+1} - t^j}}} </math>
Line 112: Line 112:
 
===Application===
 
===Application===
  
For each of the 44 LFI diodes there is the corresponding object in the Database. Because of the amplitude of the spikes we choose to apply the correction only on the 44 GHz radiometers. Each object contains 3 columns: the bins start time vector, the sky amplitudes and the reference amplitudes.  
+
For each of the 44 LFI diodes there is the corresponding object in the Database. After studying the amplitude of the spikes at other LFI frequencies, we chose to apply the correction only to the 44 GHz radiometers. Each object contains 3 columns: the bins start time vector, the sky amplitudes and the reference amplitudes.  
  
 
For each sample the value to be subtracted is computed using:
 
For each sample the value to be subtracted is computed using:
Line 120: Line 120:
 
where k is the index of the bins at a given time.
 
where k is the index of the bins at a given time.
  
==Gaps Filling==
+
==Filling Gaps in the Data==
  
During the mission some of the data packets were lost ({{PlanckPapers|planck2013-p02}}). Moreover in two different and very peculiar situations LFI was shutdown and restarted, giving inconsistencies in data sampling. All of those data aren't used for scientific purpose but to avoid discrepancies in data analysis all of the radiometers at the same frequency must have the same samples.  
+
During the mission, a small number of data packets were lost ({{PlanckPapers|planck2013-p02}}, {{PlanckPapers|planck2014-a03||Planck-2015-A03}}). Moreover in two different and very peculiar situations, LFI was shut down and restarted, giving inconsistencies in data sampling. None of these data are used for scientific purpose but to avoid discrepancies in data analysis all of the radiometers at the same frequency must have the same samples.  
  
To accomplish this the length of the data stream to be reduced in a specific pointing period is compared with the data stored in the sample information object. If the length is not the same the OBT vector is filled with missing sample times, the data vector is filled with zeros and in the flag column the bit for gap is raised.
+
To ensure this, we compare the length of the data stream to be reduced in a specific pointing period to the data stored in the sample information object. If the length is not the same, the OBT vector is filled with missing sample times, the data vector is filled with zeros, and in the flag column the bit for gap is entered.
  
 
==Gain Modulation Factor==
 
==Gain Modulation Factor==
Line 130: Line 130:
 
The pseudo-correlation design of the LFI radiometers allows a dramatic reduction of <math> 1/f </math> noise when the <math> V_{sky} </math> and <math> V_{load} </math> outputs are differenced (see [[LFI design, qualification, and performance#Overview|LFI instrument]] description). The two streams are slightly unbalanced, as one looks at the 2.7 K sky and the other looks at the ~4.5 K reference load. To force the mean of the difference to zero, the load signal is multiplied by the Gain Modulation Factor (R). For each pointing period this factor is computed using (see eq. (3) in [[LFI design, qualification, and performance#Radiometer Chain Assembly (RCA)|LFI]] description):
 
The pseudo-correlation design of the LFI radiometers allows a dramatic reduction of <math> 1/f </math> noise when the <math> V_{sky} </math> and <math> V_{load} </math> outputs are differenced (see [[LFI design, qualification, and performance#Overview|LFI instrument]] description). The two streams are slightly unbalanced, as one looks at the 2.7 K sky and the other looks at the ~4.5 K reference load. To force the mean of the difference to zero, the load signal is multiplied by the Gain Modulation Factor (R). For each pointing period this factor is computed using (see eq. (3) in [[LFI design, qualification, and performance#Radiometer Chain Assembly (RCA)|LFI]] description):
  
:<math> R = {<V_{sky}> \over <V_{load}>} </math>
+
:<math> R = { < V_{sky}> \over < V_{load}>} </math>
  
 
Then the data are differenced using:
 
Then the data are differenced using:
Line 136: Line 136:
 
:<math> TOI_{diff}^i = V_{sky}^i – R  V_{load}^i </math>
 
:<math> TOI_{diff}^i = V_{sky}^i – R  V_{load}^i </math>
  
This value for R minimizes the 1/f and the white noise in the difference timestream. The i index represents the diode and can be 0 or 1.
+
This value for R minimizes both the 1/f noise and the white noise in the difference timestream. The i index represents the diode and can be 0 or 1.
  
At this point the maneuver flag bit is set to identify which samples have missing data, using the information stored in the sampling information object. This identifies which data to ignore in the next step of the Pipeline.  
+
At this point, we also set the maneuver flag bit to identify which samples have missing data, using the information stored in the sampling information object. This identifies which data to ignore in the next step of the Pipeline.  
  
The R values are stored in the database. At the same time the mean values of <math> V_{sky} </math> and <math> V_{load} </math> are stored in order to be used in other steps of the analysis.
+
The R values are stored in the database. At the same time the mean values of <math> V_{sky} </math> and <math> V_{load} </math> are stored so they can be used in other steps of the analysis.
  
 
==Diode Combination==
 
==Diode Combination==
Line 158: Line 158:
 
The weights are fixed to a single value per diode for the entire dataset. Small variations in the relative noise of the diodes would in principle suggest recalculating the weights on shorter timescales, however, we decided a time varying weight could possibly induce more significant subtle systematics, so chose a single best estimate for the weights for each diode pair.
 
The weights are fixed to a single value per diode for the entire dataset. Small variations in the relative noise of the diodes would in principle suggest recalculating the weights on shorter timescales, however, we decided a time varying weight could possibly induce more significant subtle systematics, so chose a single best estimate for the weights for each diode pair.
  
 +
[[File:LFI_weights.jpg|400px]]
 +
:'''Weights used in combining diodes'''
 +
 +
<!--
 
{| border="1" cellpadding="5" cellspacing="0" align="center"
 
{| border="1" cellpadding="5" cellspacing="0" align="center"
 
|-
 
|-
Line 189: Line 193:
 
|-
 
|-
 
|}
 
|}
 +
-->
  
 
===Application===
 
===Application===
Line 204: Line 209:
 
===Extraction Method===
 
===Extraction Method===
  
The planets Temperature have been estimated from chunk of samples affected, plus a surrounding region, projected onto a grid (microstripes), by assuming an elliptical Gaussian beam using parameters from instrument database.
+
Measurements of planets have been formed from samples containing flux from the planet, plus a surrounding region, projected onto a grid (microstripes), by assuming an elliptical Gaussian beam using parameters from instrument database.
  
Microstripes are a way to extract and store relevant samples for planets detection. Relevant samples are samples affected by the planet plus samples in the neighbor. The search radius to select samples as relevant is 5 deg around the planet position, computed at the pointing period mid time. For each sample we store SCET (Spacecraft Event Time), pointing directions and calibrated temperature. Destriping is applied during application.
+
Microstripes are a way to extract and store relevant samples for planet detection. Relevant samples are samples affected by the planet plus samples in the neighborhood (to establish a background level). The search radius to select samples as relevant is 5 deg around the planet's position, computed at the pointing period mid time. For each sample we store SCET (Spacecraft Event Time), pointing directions and calibrated temperature. Destriping is applied.
  
Random errors are estimated by taking the variance of samples entering each micromap pixel. This is fast and the major problems (near a bright source the noise gives a larger value and it is difficult to extract the correlation matrix) causes the noise to be overestimated by a factor of two that in this situation is not a major drawback.
+
Random errors are estimated by taking the variance of samples entering each micromap pixel. This is fast but not exact, since the variance is larger near a bright source.  This can cause the noise to be overestimated by a factor of.  Given the large S/N for planetary observations, that is not a major drawback.
  
The apparent position of Planets as seen from Planck at a given time is derived from [http://ssd.jpl.nasa.gov/?horizons JPL Horizon]. Position are sampled in tables at steps of 15 minutes and then linearly interpolated at the sampling frequency of each detector. JPL Horizons tables allow also to derive other quantities such as the Planet-Planck distance and the Planet-Sun distance nad the planet angular diameter affecting the apparent brightness of the planet.
+
The apparent position of a planet as seen from Planck at a given time is derived from [http://ssd.jpl.nasa.gov/?horizons JPL Horizon]. Positions are tabulated in steps of 15 minutes and then linearly interpolated at the sampling frequency of each detector. JPL Horizons tables allow also to derive other quantities such as the planet-Planck distance and the planet-Sun distance and the planet angular diameter, which affects the apparent brightness of the planet.
  
 
The antenna temperature is a function of the dilution factor, according to:
 
The antenna temperature is a function of the dilution factor, according to:
Line 216: Line 221:
 
:<math> T_{ant,obs} = 4 log 2T_{ant,1} \left( {\theta \over b_{fwhm} } \right) ^2 </math>
 
:<math> T_{ant,obs} = 4 log 2T_{ant,1} \left( {\theta \over b_{fwhm} } \right) ^2 </math>
  
where <math> T_{ant,obs} </math> and <math> T_{ant,1} </math> are the observed and reduced <math> T_{ant} </math>, <math> \theta </math> the instantaneous planets angular diameter and <math> b_{fwhm} </math> the beam full width half maximum.
+
where <math> T_{ant,obs} </math> and <math> T_{ant,1} </math> are the observed and reduced <math> T_{ant} </math>, <math> \theta </math> the instantaneous angular diameter of the planet and <math> b_{fwhm} </math> the beam full width at half maximum.
  
With the above definition <math> T_{ant,1} </math> could be considered as the <math> T_{ant} </math> for a planet with <math> b_{fwhm} = \theta </math>, but a more convenient view is to take a Reference Dilution factor <math> D_0 </math>, as the dilution factor for a standardized planet angular diameter and beam fwhm <math> b_{fwhm} </math>, <math> \theta_0 </math>, to have:
+
With the above definition <math> T_{ant,1} </math> could be considered as the <math> T_{ant} </math> for a planet with <math> b_{fwhm} = \theta </math>, but a more convenient view is to take a reference dilution factor <math> D_0 </math>, as the dilution factor for a standardized angular diameter for the planet and fiducial beam fwhm <math> b_{fwhm} </math>, <math> \theta_0 </math>, to have:
  
 
:<math> D_0 = \left( {\theta_0 \over b_0 } \right) ^2 </math>
 
:<math> D_0 = \left( {\theta_0 \over b_0 } \right) ^2 </math>
Line 230: Line 235:
 
===Application===
 
===Application===
  
The OBT vector found by the search are saved in a set of object, one for each horn. In Level2 Pipeline those OBTs are compared with the OBT vector of the data to raise planet bit flag where needed.
+
The OBT vectors found by the search are saved in a set of objects, one for each horn. In the Level 2 pipeline those OBTs are compared with the OBT vector of the data to set the planet bit flag where needed.
  
 
==Photometric Calibration==
 
==Photometric Calibration==
  
Photometric calibration is the procedure used to convert data from volts to kelvin. The source of the calibration is the well known CMB dipole, caused by the motion of the Solar System with respect to the CMB reference frame. To this signal we add the modulation induced by the orbital motion of Planck around the Sun. The resulting signal is then convoluted with the horn beam to get the observed Dipole.
+
Photometric calibration is the procedure used to convert data from volts to kelvin. The source of the calibration is the well known CMB dipole, caused by the motion of the Solar System with respect to the CMB reference frame. To this signal we add the modulation induced by the orbital motion of Planck around the Sun. The resulting signal is then convoluted with the horn beam to get the observed dipole. For details refer to {{PlanckPapers|planck2014-a06||Planck-2015-A06}}.  
  
 
===Beam Convolved Dipole===
 
===Beam Convolved Dipole===
Line 244: Line 249:
 
where <math> \mathbf{P}_E(t) </math> is the pointing direction, in the observer reference frame and <math> \mathbf{D}_E </math> is the dipole axis scaled by the dipole amplitude again in the same reference frame.
 
where <math> \mathbf{P}_E(t) </math> is the pointing direction, in the observer reference frame and <math> \mathbf{D}_E </math> is the dipole axis scaled by the dipole amplitude again in the same reference frame.
  
In general the true signal would have to be convolved with the beam pattern of the given radiometer, usually described as a fixed map in the beam reference frame or as a time dependent map in the observer reference frame. In this case it is easiest to describe the convolution in the beam reference frame, since the function to be convolved is described by a single vector.
+
In general the true signal would have to be convolved with the beam pattern of the given radiometer, usually described as a fixed map in the beam reference frame or as a time dependent map in the observer reference frame. In this case it is easiest to describe the convolution in the beam reference frame, since the function to be convolved then can be described by a single vector.
  
 
Denoting with <math> \mathcal{U}(t) </math> the matrix converting from the observer to the beam reference frame, so that:
 
Denoting with <math> \mathcal{U}(t) </math> the matrix converting from the observer to the beam reference frame, so that:
Line 258: Line 263:
 
:<math> \Delta T_D(t) = N \int_{4\pi} B(\mathbf{P})\mathbf{P} \cdot \mathbf{D}(t) d^3\mathbf{P} </math>
 
:<math> \Delta T_D(t) = N \int_{4\pi} B(\mathbf{P})\mathbf{P} \cdot \mathbf{D}(t) d^3\mathbf{P} </math>
  
where <math> \mathbf{N} </math> is a normalization costant.
+
where <math> \mathbf{N} </math> is a normalization constant.
  
 
:<math> N^{-1} = \int_{4\pi} B(\mathbf{P}) d^3\mathbf{P} </math>
 
:<math> N^{-1} = \int_{4\pi} B(\mathbf{P}) d^3\mathbf{P} </math>
Line 277: Line 282:
 
! <math> S_z </math>
 
! <math> S_z </math>
 
|-
 
|-
| LFI18S || 1.4105692317321994e-03 || -3.7689062388084022e-04 || 9.9999893412338192e-01
+
| LFI18S || 4.6465302801434851e-03 || -1.3271241134658798e-03 || 9.9577083251002374e-01  
 
|-
 
|-
| LFI18M || 1.1200251268914613e-03 || -3.2838598619563524e-04 || 9.9999931885294768e-01
+
| LFI18M || 2.8444949872572308e-03 || -8.9511701920330633e-04 || 9.9735308499993403e-01  
 
|-
 
|-
| LFI19S || 1.7861136968831050e-03 || -4.4036975450455066e-04 || 9.9999830793473898e-01
+
| LFI19S || 4.2948823053610341e-03 || -1.4037771886936600e-03 || 9.9612844031547154e-01  
 
|-
 
|-
| LFI19M || 1.4292780457919835e-03 || -4.7454175238335579e-04 || 9.9999886598655352e-01
+
| LFI19M || 4.5097043389558588e-03 || -1.6255119872940569e-03 || 9.9573258887316529e-01  
 
|-
 
|-
| LFI20S || 1.7008692096818349e-03 || -6.1036624911600191e-04 || 9.9999836724715374e-01
+
| LFI20S || 5.2797843121430788e-03 || -1.9479348764780897e-03 || 9.9514557074549859e-01  
 
|-
 
|-
| LFI20M || 1.5548897911626446e-03 || -5.9289001736737262e-04 || 9.9999861539862389e-01
+
| LFI20M || 4.7420554683343680e-03 || -1.8462064956921880e-03 || 9.9548302435765323e-01  
 
|-
 
|-
| LFI21S || 1.6975720932854463e-03 || 6.0961185087824777e-04 || 9.9999837330986663e-01
+
| LFI21S || 5.2377151034493320e-03 || 1.9337533884005056e-03 || 9.9517552967578515e-01  
 
|-
 
|-
| LFI21M || 1.5486274949897787e-03 || 5.9228926426513112e-04 || 9.9999862547220986e-01
+
| LFI21M || 4.3713123245201109e-03 || 1.7075534076511198e-03 || 9.9585142295482576e-01  
 
|-
 
|-
| LFI22S || 1.7861136968831245e-03 || 4.4036975450366470e-04 || 9.9999830793473898e-01
+
| LFI22S || 3.6272964935149194e-03 || 1.2030860597901591e-03 || 9.9679007501368400e-01  
 
|-
 
|-
| LFI22M || 1.4292780457920242e-03 || 4.7454175238250377e-04 || 9.9999886598655352e-01
+
| LFI22M || 3.1817723013115090e-03 || 1.1383448292256770e-03 || 9.9705527860728405e-01  
 
|-
 
|-
| LFI23S || 1.4105692317321714e-03 || 3.7689062387997129e-04 || 9.9999893412338203e-01
+
| LFI23S || 3.2582522957523585e-03 || 9.1080217792813213e-04 || 9.9706066300815721e-01  
 
|-
 
|-
| LFI23M || 1.1200251268914476e-03 || 3.2838598619481239e-04 || 9.9999931885294757e-01
+
| LFI23M || 2.6384956780749246e-03 || 8.0545189089709709e-04 || 9.9755224385088148e-01  
 
|-
 
|-
| LFI24S || 3.4636411743209074e-04 || -2.8530917087092225e-07 || 9.9999994001590664e-01
+
| LFI24S || 1.2007669925034473e-03 || -2.9877396130232452e-07 || 9.9897924004634964e-01  
 
|-
 
|-
| LFI24M || 4.3939553230170735e-04 || -2.9414231975370517e-07 || 9.9999990346573508e-01
+
| LFI24M || 1.1966952022302408e-03 || -1.3590311231894260e-06 || 9.9894545519286426e-01  
 
|-
 
|-
| LFI25S || -1.0428719495964051e-04 || 1.9328051933678115e-04 || 9.9999997588341061e-01
+
| LFI25S || -2.1006042848582498e-04 || 4.4287041964311393e-04 || 9.9965683952381934e-01  
 
|-
 
|-
| LFI25M || -1.1004766833423990e-04 || 2.7656488668259429e-04 || 9.9999995570068612e-01
+
| LFI25M || -2.5133433631240997e-04 || 5.5803053337715499e-04 || 9.9954614333864866e-01  
 
|-
 
|-
| LFI26S || -1.0428719495970346e-04 || -1.9328051933760877e-04 || 9.9999997588341061e-01
+
| LFI26S || -1.9766749938836072e-04 || -4.2202164851420096e-04 || 9.9966955259454382e-01  
 
|-
 
|-
| LFI26M || -1.1004766833430009e-04 || -2.7656488668343130e-04 || 9.9999995570068612e-01
+
| LFI26M || -2.5426738260127387e-04 || -5.7426678736881309e-04 || 9.9952104084618332e-01  
 
|-
 
|-
| LFI27S || 1.6613273546973915e-03 || 6.6518363019636186e-04 || 9.9999839875979735e-01
+
| LFI27S || 5.8555501689062294e-03 || 1.8175261751324087e-03 || 9.9448302844644199e-01  
 
|-
 
|-
| LFI27M || 1.5583345016298123e-03 || 6.4183510236962536e-04 || 9.9999857981963269e-01
+
| LFI27M || 4.9962842224387377e-03 || 1.5766304210447924e-03 || 9.9529240281463061e-01  
 
|-
 
|-
| LFI28S || 1.6633788116048607e-03 || -6.6629002345089925e-04 || 9.9999839461297824e-01
+
| LFI28S || 6.3750745068891614e-03 || -1.9617207193203534e-03 || 9.9396633709784610e-01  
 
|-
 
|-
| LFI28M || 1.5571200481047094e-03 || -6.4144198187461837e-04 || 9.9999858196366442e-01
+
| LFI28M || 4.7877703093970429e-03 || -1.5442874049905993e-03 || 9.9551694709996730e-01  
 
|-
 
|-
 +
 
|}
 
|}
  
  
  
By using this characteristic vector the calculation of the convolved dipole is simply defined by a dot product of the vector <math> \mathbf{S} </math> by the dipole axis rotated in the beam reference frame.
+
By using this characteristic vector the calculation of the convolved dipole is simply defined by a dot product of the vector <math> \mathbf{S} </math> and the dipole axis rotated in the beam reference frame.
  
 
:<math> \Delta T (t) = \mathbf{S}^T \mathcal{U}(t) \mathbf{D}_E </math>
 
:<math> \Delta T (t) = \mathbf{S}^T \mathcal{U}(t) \mathbf{D}_E </math>
Line 331: Line 337:
 
===Binning===
 
===Binning===
  
In order to simplify the computation and to reduce the amount of data used in the calibration procedure the data are phase binned in map with Nside 256. During phase binning all the data with flagged for maneuvers, planets, gaps and the ones flagged in Level1 analysis as not recoverable are discharged.
+
In order to simplify the computation and to reduce the amount of data used in the calibration procedure the data are phase binned in map with <math>N_{side}</math> 256. During phase binning all the data flagged for maneuvers, planets and gaps, as well as the ones flagged in Level 1 analysis as not recoverable, are ignored.
  
 
===Fit===
 
===Fit===
  
The first order calibration values are given by a Least Square Fit between the signal and the dipole. For each pointing a gain (<math> g_k </math>) and an offset (<math> b_k </math>) values are computed minimizing:
+
The first order calibration values are given by a least squares fit between the signal and the dipole. For each pointing gain (<math> g_k </math>) and offset (<math> b_k </math>) values are computed by minimizing:
  
 
:<math> \chi^2 = \sum_{i \in k} {[ \Delta V (t_i) - \Delta V_m (t_i|g_k, b_k)]^2 \over rms_i^2 } </math>
 
:<math> \chi^2 = \sum_{i \in k} {[ \Delta V (t_i) - \Delta V_m (t_i|g_k, b_k)]^2 \over rms_i^2 } </math>
  
The sum includes samples outside a Galactic mask.
+
The sum includes only samples outside a Galactic mask.
 +
 
 +
===DaCapo===
 +
 
 +
The largest source of error in the fit arises from unmodeled sky signal <math> \Delta T_a </math> from CMB anisotropy.
 +
 
 +
The procedure adopted to correct this effect is described below.
 +
 
 +
The uncalibrated toi <math>y_{i}</math> belonging to a pointing period  <math>k</math> is modeled as 
 +
:<math>y_{i}=g_{k}*(s_{i}+d_{i})+b_{k}+n</math>
 +
 
 +
where the unknowns <math>s_{i}</math>,<math>g_{k}</math> and <math>b_{k}</math> are respectively the signal per OBT, the gain per PID and the offset per PID, and <math>d_{i}</math> is the total dipole.
 +
 
 +
This is a non linear <math>\chi^{2}</math>  minimization problem. It can be linearized with an iterative procedure: for each iteration step a linearized model toi is built as
 +
:<math>y_{i,model}=g_{k}(-1)*(s_{i}-s_{i}(-1))+g_{k}*(d_{i}+a_{i}(-1))+b_{k}</math>
  
===Mademoiselle===
+
where <math>g_{k}(-1)</math> and <math>s_{i}(-1)</math>  are the gain and sky signal computed at the previous step of iteration.
  
The largest source of error in the fit arises from unmodeled sky signal <math> \Delta T_a </math> from CMB anisotropy. To correct this we iteratively project the calibrated data (without the dipole) onto a map, scan this map to produce a new TOD with astrophysical signal removed, and finally run a simple destriping algorithm to find the corrections to the gain and offset factors.
+
A linear fit is performed between <math>y_{i}</math> and <math>y_{i,model}</math> to get the new <math>g_{k}</math>  and <math>s_{i}</math> . The procedure is iterated until convergence.  
  
 
To reduce the impact of the noise during the iterative procedure the sky estimation is built using data from both radiometers of the same horn.
 
To reduce the impact of the noise during the iterative procedure the sky estimation is built using data from both radiometers of the same horn.
Line 349: Line 369:
 
===Smoothing===
 
===Smoothing===
  
To improve accuracy given by the iterative algorithm and remove noise from the solution a smoothing algorithm must be performed. We used two different algorithms: OSG for the 44 and 70 GHz radiometers, and DV/V Fix for the 30 GHz. The reasons behind this choice can be found in {{PlanckPapers|planck2013-p02b}}.
+
To improve accuracy given by the iterative algorithm and remove noise from the solution a smoothing algorithm must be performed.  
 +
 
 +
====OSGTV====
 +
 
 +
OSGTV is a 3 step smoothing algorithm, implemented with a C++ code.  
  
====OSG====
+
The gain reconstructed with DaCapo can be expressed as
 +
:<math>G_{raw}(t,T(t))=G_{true}(t,T(t))+G_{noise}(t)</math>
  
OSG is a python code that performs smoothing with a 3 step algorithm.  
+
where the <math>T(t)</math> function represents temperature fluctuations of the Focal Plane.
 +
The time dependence of the “real” gain <math>G_{true}</math>  is modeled as the superposition of a “slow” component, with a time scale of ~3 months, and a “fast” component with a time scale of few PIDs:
 +
:<math>G_{true}=G_{slow}+G_{fast}</math>.
  
The first step is a Moving Average Window: the gain and offset factors are streams containing one value for each pointing period, that we call dipole fit raw streams. The optimized window has a length of 600 pointing periods.
+
The slow component  takes into account the seasonal variations of the thermal structure of the spacecraft due to the orbital motion, while the “fast” component describes the thermal effects of the electronics and compressors, as well as single events like the sorption cooler switchover.
  
The second step is a wavelet algorithm, using pywt ([http://www.pybytes.com/pywavelets/ Discrete Wavelet Transform in Python]) libraries. Both dipole fit raw streams and averaged streams are denoised using wavelets of the Daubechies family extending the signals using  symmetric-padding.
+
To disentangle the components of  <math>G_{true}</math>, we need a parametric “hybrid” approach in three steps.
  
The third step is the combination of dipole fit raw and averaged denoised signal using knowledge about the instrument performance during the mission.
+
<b>Step 1. </b> For each PID <math>i</math>, we used the 4K total-power <math>V_{load}(i)</math> and the signal from temperature detectors in the focal plane <math>T_{i}</math>, subsampled at 1 sample/PID, to track gain changes. This is implemented through a linear fit between <math>V_{load}(i) * G_{raw}</math> and <math>T_{i}</math>:
  
====4 K total-power and Fix====
+
:<math>V_{load}(i)*G_{raw}=K_{0}+K_{1}*(T_{i}-\langle T_{MW}\rangle)</math>
  
For the 30 GHz channels we used 4K total-power to track gain changes. The theory and explanation of the choice can be found in {{PlanckPapers|planck2013-p02b}}.
+
where the window length <math>MW</math> of the moving average is proportional to the variance of the dipole in the considered PID.
 +
The resulting gain is:
 +
:<math> G_{i} = \frac{ K_{0} + K_{1} * (T_{i} - \langle T_{MW} \rangle ) }{ V_{load}(i) } </math>
  
The algorithm uses <math> V_{load} </math> mean values computed during differentiation and raw gains as they are after iterative calibration, performing a linear weighted fit between the two streams using as weight the dipole variance in single pointing periods. The fit is a single parameter fit, so the offsets are put to zero in this smoothing method. It uses the gsl libraries.
+
<b>Step 2. </b> The "fast" component <math>G_{fast}</math> is recovered as follows: we define a maximum moving average window length <math>WL_{max}</math>,and for each PID we compute the variance <math>\sigma_{i}^2</math> of <math>G_{i}</math> in this window.
 +
We define a percentile on the ordered variance array and we compute the corresponding value of the variance <math>\sigma_{p}^2</math>. The window length for each PID is then computed as
 +
:<math>WL_{i}=\frac{\sigma_{p}^2*WL_{max}}{\sigma_{i}^2}</math>.
  
In addition to the smoothing, to better follow sudden gain changes due to instrument configuration changes, a fix algorithm is implemented. The first step is the application of the 4k total-power smoothed gains to the data and the production of single radiometer maps in the periods between events. The resulting maps are then fit with dipole maps covering the same period of time producing two factor for each radiometer: <math> corrM </math> is the result of the fit using the main  radiometer and <math> corrS </math> the one coming from the side radiometer. The correction to be applied to the gain values is then computed as:
+
If <math>WL_{i}\gt WL_{max}</math> we impose <math>WL_{i}=WL_{max}</math>. With these window lengths
 +
a moving window average is performed on <math>G_{i}</math>.
 +
The averaged gain vector is subtracted from the raw hybrid gain to get <math>G_{fast}</math>.
  
:<math> corr = {1 \over {1 + {{corrM+corrS} \over 2 }}} </math>
+
<b>Step 3. </b> We perform a moving window average on <math>G_{i}</math> and we compute the variances of the smoothed array. The window length is computed with a linear interpolation between a minimum length <math>WL_{min}</math> defined in the dipole minima and a maximum <math>WL_{max}</math> defined in the dipole maxima. The array of gain variances is weighted with the variance of the dipole.
 +
We set a percentile on the variance array and we find the corresponding variance value <math>\sigma_{p}^2</math>. Around this value we search for local maxima of the variance array, and we split the domain of the gain in subsets between consecutive maxima.
 +
For each subset we perform a moving average with the corresponding window length.
 +
The "slow" component <math>G_{slow}</math> is given by the union of these subsets.
  
 
===Gain Application===
 
===Gain Application===
Line 375: Line 411:
 
The last step in TOI processing is the creation of the calibrated stream. For each sample we have:
 
The last step in TOI processing is the creation of the calibrated stream. For each sample we have:
  
:<math> TOI_{cal}(t) = (TOI_{diff}(t) – offset(k)) g(k) convDip(t) </math>
+
:<math> TOI_{cal}(t) = (TOI_{diff}(t) – offset(k)) \cdot g(k) \, - \, convDip(t)\,  - \, T_{stray}(t)</math>
 +
 
 +
where t is the time and k is the pointing period, <math> convDip </math> is the CMB Dipole convolved with the beam, and <math>T_{stray}</math> is the straylight.
  
where t is the time and k is the pointing period. <math> convDip </math> is the CMB Dipole convolved with the beam.
 
 
== Noise ==
 
== Noise ==
This pipeline step aims at the reconstruction of the noise parameters from calibrated flight TOI. The goal is two-folds: one the one side we need to know the actual noise properties of
+
This pipeline step aims at the reconstruction of the noise parameters from calibrated flight TOI. The goal is two-fold.  On the one hand we need to know the actual noise properties of
our instrument in order to properly take them into account especially during the following processing and analysis steps like map-making and power spectrum estimation. On the other
+
our instrument in order to properly take them into account, especially during later processing and analysis steps like map-making and power spectrum estimation. On the other
side evaluation of noise properties along the instrument life-time is a way to track down possible variations, anomalies and general deviations from the expected behaviour.
+
hand evaluation of noise properties during the instrument life-time is a way to track down possible variations, anomalies and general deviations from the expected behaviour.
  
 
=== Operations ===
 
=== Operations ===
Noise estimation is performed on calibrated data and since we would like to track possible noise variations along mission life-time, we
+
Noise estimation is performed on calibrated data; since we would like to track possible noise variations during the mission life-time, we
 
select data in chunks of 5 ODs (Operational Days). These data are processed by the ROMA Iterative Generalized Least Square (IGLS)
 
select data in chunks of 5 ODs (Operational Days). These data are processed by the ROMA Iterative Generalized Least Square (IGLS)
map-making algorithm which includes a noise estimation tool. In general an IGLS map-making is a quite consuming in terms of time and  
+
map-making algorithm which includes a noise estimation tool. In general IGLS map-making is a quite costly in terms of time and  
resources required. However the length of the data is such that running on the DPC cluster in very short time (~1-2 minutes).
+
resources required. However the length of the data is such that it can run on the DPC cluster in very short time (~1-2 minutes).
  
 
The method implemented can be summerized as follows. We model the calibrated TOI as
 
The method implemented can be summerized as follows. We model the calibrated TOI as
:$\mathbf{\Delta T} = \mathbf{P} \mathbf{m} + \mathbf{n}$
+
:<math> \mathbf{\Delta T} = \mathbf{P} \mathbf{m} + \mathbf{n} </math>
where $\mathbf{n}$ is the noise vector and $\mathbf{P}$ is the pointing matrix that links a pixel in the map $\mathbf{m}$
+
where <math> \mathbf{n} </math> is the noise vector and <math> \mathbf{P} </math> is the pointing matrix that links a pixel in the map <math> \mathbf{m} </math>
with a sample in the TOI $\mathbf{d}$. The zero-th order estimation of the signal is obtained simply rebinning TOI into a map.
+
with a sample in the TOI <math> \mathbf{d} </math>. The zero-th order estimation of the signal is obtained by simply rebinning TOI into a map.
 
Then an iterative approach follows in which both signal and noise are estimated according to
 
Then an iterative approach follows in which both signal and noise are estimated according to
:$\mathbf{\hat{n}_i} = \mathbf{\Delta T} - \mathbf{P\hat{m}_i}$
+
:<math> \mathbf{\hat{n}_i} = \mathbf{\Delta T} - \mathbf{P\hat{m}_i} </math>
: $\mathbf{\hat{m}_{i+1}} = \mathbf{(P^T\hat{N}_i^{-1}P)^{-1}P^T\hat{N}_i^{-1}\Delta T}$
+
: <math> \mathbf{\hat{m}_{i+1}} = \mathbf{(P^T\hat{N}_i^{-1}P)^{-1}P^T\hat{N}_i^{-1}\Delta T} </math>
where $\mathbf{\hat{N}_i}$ is the noise covariance matrix in time domain out from iteration $i$. After three iterations convergence is  
+
where <math> \mathbf{\hat{N}_i} </math> is the noise covariance matrix in the time domain resulting from iteration <math>i </math>. After three iterations, convergence is  
 
achieved.
 
achieved.
  
We then perform an FFT (Fast Fourier Transform) on the noise time stream out from the iterative approach and then fit the resulting spectrum.
+
We then perform an FFT (Fast Fourier Transform) on the noise time stream from the iterative approach and fit the resulting spectrum.
  
 
=== Fitting Pipeline ===
 
=== Fitting Pipeline ===
In the very first release of Planck data, once noise spectra were extracted a simply log-periodogram fitting approach was applied to derive the
+
As already done in the 2013 release, we estimate the parameters of the noise properties  (i.e. white noise level, knee-frequency and slope of the
most important noise parameters (white noise level, knee-frequency and slope of the low-frequency noise component). However during mission life-time
+
low-frequency noise component) by means of a MCMC approach. Therefore on the spectra just described we first compute the white noise level taking the mean of the last 10% of data (at 30 GHz due to the higher value of the
there were some specific events (e.g. the switch over of the sorption coolers) that we expect were able to cause variation in instrument behaviour and hence in
+
knee-frequency this quantity is reduced to 5%). Once the white noise level has been determined we proceed with the actual fitting
its noise properties. In this respect we have improved our fitting pipeline adding a Monte Carlo Markov Chain approach to estimate noise parameters.
+
of knee-frequency and slope. The resulting values reported are the medians of the fitted values for our 5 ODs chunk along the
 +
whole mission lifetime.
  
==== MCMC approach ====
+
{| border="1" cellpadding="5" cellspacing="0" align="center"
This new approach allows us to improve our noise model. Indeed this can be parametrized by the usual combination of white plus <math> 1/f </math> noise
+
|-
:<math> P(f)=\sigma^2\left[1+\left(\frac{f}{f_k}\right)^\beta\right] </math>
+
! Horn
with three basic noise parameters.
+
! White Noise M [<math>\mu</math>K s<math>^{1/2}</math>]
However it is also possible to work with  a
+
! White Noise S [<math>\mu</math>K s<math>^{1/2}</math>]
functional form with two more parameters as
+
! Knee-frequency M [mHz]
:<math> P(f) = \sigma^2\left[1 + \left(\frac{f}{f_{k1}}\right)^{\beta_{1}} +\left(\frac{f}{f_{k2}}\right)^{\beta_{2}}\right] </math>
+
! Knee-frequency S [mHz]
 
+
! Slope M
This latter could be useful when there are clearly two different
+
! Slope S
behaviour in the low-frequency part of the spectrum where, beside
+
|-
usual radiometric <math> 1/f </math> noise, appears signature of thermal
+
| 18 || 513.0<math>\pm</math>2.1 || 467.2<math>\pm</math>2.3 || 14.8<math>\pm</math>2.5 || 17.8<math>\pm</math>1.5 || -1.06<math>\pm</math>0.10 || -1.18<math>\pm</math>0.13
fluctuations induced noise.
+
|-
 
+
| 19 || 579.6<math>\pm</math>2.2 || 555.0<math>\pm</math>2.2 || 11.7<math>\pm</math>1.2 || 13.7<math>\pm</math>1.3 || -1.21<math>\pm</math>0.26 || -1.11<math>\pm</math>0.14
As for the white noise part, this is, as
+
|-
before, computed making a simple average of noise spectrum on the last 10% of
+
| 20 || 587.3<math>\pm</math>2.1 || 620.5<math>\pm</math>2.7 || 8.0<math>\pm</math>1.9 || 5.7<math>\pm</math>1.5 || -1.20<math>\pm</math>0.36 || -1.30<math>\pm</math>0.41
frequency bins. This percentage works well for almost all radiometers
+
|-
at 44 and 70 GHz but it is indeed quite delicate for the
+
| 21 || 451.0<math>\pm</math>1.7 || 560.1<math>\pm</math>2.0 || 37.9<math>\pm</math>5.2 ||13.3<math>\pm</math>1.5 || -1.25<math>\pm</math>0.09 ||-1.21<math>\pm</math>0.09
30 GHz radiometers which show typical values of knee-frequency around
+
|-
100 mHz and, therefore, require a smaller number to get an un-biased
+
| 22 || 490.8<math>\pm</math>1.5 || 531.3<math>\pm</math>2.3 || 9.7<math>\pm</math>2.3 ||14.8<math>\pm</math>6.7 || -1.42<math>\pm</math>0.23 || -1.24<math>\pm</math>0.30
white noise estimation.
+
|-
Once white noise is computed, the code creates Markov Chains for the
+
| 23 || 504.3<math>\pm</math>1.8 || 539.7<math>\pm</math>1.8 || 29.7<math>\pm</math>1.1 ||59.0<math>\pm</math>1.4 || -1.07<math>\pm</math>0.03 || -1.21<math>\pm</math>0.02
other parameters. Discarding the burn-in period of the chains we can
+
|-
directly get from the chain samples distribution, the expected value
+
| 24 ||463.0<math>\pm</math>1.4 || 400.7<math>\pm</math>1.3 || 26.8<math>\pm</math>1.3 ||88.3<math>\pm</math>8.9 || -0.94<math>\pm</math>0.01 || -0.91<math>\pm</math>0.01
and variance of each noise parameters sampled.
+
|-
 +
| 25 || 415.3<math>\pm</math>1.5 || 395.4<math>\pm</math>2.9 || 20.1<math>\pm</math>0.7 ||46.4<math>\pm</math>1.8 || -0.85<math>\pm</math>0.01 || -0.90<math>\pm</math>0.01
 +
|-
 +
| 26 || 483.0<math>\pm</math>1.9 || 423.2<math>\pm</math>2.5 || 64.4<math>\pm</math>1.9 ||68.2<math>\pm</math>9.5 || -0.92<math>\pm</math>0.01 || -0.76<math>\pm</math>0.01
 +
|-
 +
| 27 || 281.5<math>\pm</math>2.1 || 303.2<math>\pm</math>1.8 || 174.5<math>\pm</math>2.9 ||108.8<math>\pm</math>2.5 ||-0.93<math>\pm</math>0.01|| -0.91<math>\pm</math>0.01
 +
|-
 +
| 28 || 317.1<math>\pm</math>2.4|| 286.5<math>\pm</math>2.3 || 130.1<math>\pm</math>4.4 ||43.1<math>\pm</math>2.4 ||-0.93<math>\pm</math>0.01|| -0.90<math>\pm</math>0.02
 +
|-
 +
|}
  
The left panel of the following Figure shows a typical spectrum at 70 GHz with
+
Time variations of noise parameters are a good tracer of possible modifications in the instrument behaviour and we know that
superimposed the simple log-periodogram fit (purple line) and the new
+
some events capable of affecting instrument behaviour had happened during the mission. Both variations of the physical temperature
MCMC derived spectrum (blue line). The right panel instead shows
+
of the instrument due to the transition in the operations from the first sorption cooler to the second one, as well as the
distribution for knee-frequency and slope derived from the example spectrum.
+
observed degradation in the performances of the first cooler are events that clearly show their fingerprints in the
 +
variation of the noise spectra.
  
 +
In the following figure we report a representative sample of noise spectra, one for each frequency channel (from left to right LFi 18M, LFI 25S and LFI 27M), covering
 +
the whole mission lifetime. The white noise is extremely stable at the level of 0.3%. As already noted in the
 +
2013 release, knee-frequencies and slopes are stable until OD 326 and show clear variations and deviations from the
 +
simple one knee-frequency one slope model. This is a sign of the degradation of the first cooler inducing thermal variations
 +
with a characteristic knee-frequency (different from the radiometric one) and with a steeper slope. Once the
 +
second cooler became operational and performed as expected, the noise spectra gradually returned to the original
 +
shape at the beginning of the mission. This behaviour is visible, although not identical at each frequency for
 +
several reasons e.g. intrinsic thermal susceptibility and position on the focal plane that determines the
 +
actual thermal transfer function.
  
=== The final noise parameters ===
 
As already reported we know that during the nominal operations there was a quite dramatic change in LFI
 
induced by the switch over of the two sorption coolers and particularly we expect to see the effect of degradation of the
 
performance of the first sorption cooler and the onset of the redundant one.
 
  
In the following figure we report a set of noise frequency spectra for three LFI radiometers (LFI28M, LFI24S and LFI18M)
+
[[File:noise_spectra_LFI18M_DX11D_legend.png|300px|'"'LFI 18M.'"']]
from the beginning of the operation till the time of the current data release. Some comments are in order.
+
[[File:noise_spectra_LFI25S_DX11D_legend.png|300px |'"'LFI 25S.'"']]
First of all the white noise level is extremely stable in all the three cases (but this is also true for all the LFI  
+
[[File:noise_spectra_LFI27M_DX11D_legend.png|300px|'"'LFI 27M.'"']]
radiometer). Also knee-frequency and low-frequency slope are quite stable till OD 326. After that period
 
spectra show a
 
noise increase and two slopes for the low-frequency part which become
 
more evident for spectra around OD 366 and OD 466 where the first
 
cooler starts to be less effective and produces low-frequency thermal
 
noise. After the switch-over to the redundant cooler data still
 
present (the very last spectrum) thermal noise at very low-frequency.
 
This behaviour is almost present in all radiometers with different
 
trends ranging from the small effect shown by LFI24S to more prominent
 
effect as shown by LFI28M and LFI18M.
 
  
 
== References ==
 
== References ==

Latest revision as of 15:34, 4 February 2015

Overview[edit]

The LFI Level2 Pipeline analyzes data from each horn of the instrument separately, one pointing period at time, and stores the results in an object the length of an OD. Each diode of the horn is corrected for systematic effects. Next, measurements of the sky and the 4K load are differenced, then the signals from one diode are combined with signals from the complementary diode in the same radiometer. Finally, photometric calibration is applied for each horn.

Pre-processing[edit]

Before the Level 2 pipeline is run, the Mission information and data sampling divisions are stored in the database, in order to improve the analysis.

The Mission information is a set of objects, one for each Operational Day (OD, as defined in the glossary), in which are stored pointing period data: the DPC pointing ID (where 1 is the first pointing of the nominal mission), PSO pointing ID, start OBT of the pointing maneuver, start OBT of the stable pointing, OBT of the end of the stable pointing, and spin axis ecliptic longitude and latitude.

The sampling information is a set of objects, one for each LFI frequency, in which are stored for each pointing ID: start OBT of the pointing maneuver, start OBT of the stable pointing, end OBT of the pointing, number of samples of the pointing, number of stable samples of the pointing, start sample of the stable pointing and sample number from the start of the nominal mission. Valid samples and OBTs are defined to be those where any of the radiometers in a frequency band contain valid data.

ADC Correction[edit]

During the analysis of LFI date, we discovered a systematic common to both white noise and calibration. It turn out to be a non linearity in the Analogue to Digital Converter (ADC) on board. More detail is supplied in Planck-2013-II[1]Planck-2013-III[2], Planck-2015-A03[3] and Planck-2015-A04[4].

Evaluation[edit]

The mathematical model represents the digital ADC output as:

[math] X = (V - \Delta) \gamma + x_0 [/math]

where [math] V [/math] is the voltage input, [math] \gamma [/math] is the DAE gain, [math] \Delta [/math] is the DAE offset and [math] x_0 [/math] is the DAE [math] T_{zero} [/math].

We can model the non-linearity as a function of the input voltage [math] R(V) [/math]. So we have the apparent inferred voltage [math]V'[/math] and we can link it to the actual input voltage with:

[math] ((V - \Delta) \gamma + x_0) R(V) = X = ((V' - \Delta) \gamma + x_0 [/math]

so that:

[math] R(V) = (V' - \Delta) \gamma + x_0 \over (V - \Delta) \gamma + x_0 [/math]

Since [math] V \gt \gt \Delta [/math] and [math] V \gamma \gt \gt x_0 [/math] we can use the much simpler relation:

[math] R(V) = \frac{V'}{V} [/math]

and we expect it to be very near to unity for all [math] V [/math].

To find the response curve we have only the apparent voltage to work with, so we had to use the inverse response function [math] R'(V') [/math] and replace the real input voltage with [math] T_{sys} [/math] times the time varying gain factor [math] G(t) [/math].

[math] V = V'R'(V') = G(t)T_{sys} [/math]

If we introduce a small signal on top of [math] T_{sys} [/math] which leads to increased detected voltage and corresponding apparent voltage increment:

[math] V + \delta V = (V' + \delta V') R' (V' + \delta V') = G(t) (T_{sys} + \delta T) [/math]

so carrying out the differentiation with respect to [math] V' [/math] to the relation between true and apparent signal voltage leads to:

[math] \delta V = \left( V' {dR'(V') \over dV' + R'(V')} \right) \delta V' = G(t) \delta T [/math]

We now assume [math] T_{sys} [/math] and [math] \delta T [/math] are fixed and that the variations are due to slow drifts in the gain. So we can isolate the terms:

[math] V' = {G(t) T_{sys} \over R'(V')} [/math]
[math] \delta V' = {G(t) \delta T \over V' {dR'(V') \over dV} + R'(V')} [/math]

Combining the equations through the gain factor to remove it:

[math] {V' R'(V') \over T_{sys}} = {\delta V' {dR'(V') \over dV'} + R'(V') \over \delta T } [/math]

Rearranging and putting [math] a = {\delta T \over T_{sys}} [/math]

[math] {dR'(V') \over dV' } = \left( {a \over dV'} - {1 \over V'} \right) R'(V') [/math]

So there is the expected direct proportionality of [math] \delta V' [/math] to [math] V' [/math] based on the assumption that the variations in voltage are due to overall gain drift, so the amplitude of voltage and signal will vary together. Then there is the additional differential term which will pull the signal amplitude away from the linear relationship. So if we plot measured white noise or dipole gain factor against recovered voltage, we should see a linear curve with variations due to local slope changes at particular voltages. The linear part can be taken out and the differential part fitted. This was numerically integrated up to get the inverse response curve, which we then used to convert the measured voltages to corrected voltages.

Application[edit]

For each of the 44 LFI diodes there is the corresponding corrected object in the Database. Each object contains 4 columns: the input voltages coming from the sky channel and the corresponding linearized output, the input voltages coming from the 4K reference channel and the corresponding linearized output.

Data loaded by the module are used to initialize two different interpolators using CSPLINE and the functions from gsl (GNU Scientific Libraries) libraries. The interpolators are then used to correct each sample.

Spikes Removal[edit]

Some of the LFI receivers exhibit a small artifact with exactly 1 second period, that produces effects visible in the power spectra. The effect is a set of spikes at 1 Hz and harmonics. The spurious signal is very well modeled and is removed from the timelines. More information can be found in Planck-2013-II[1]Planck-2013-III[2]Planck-2015-A03[3]Planck-2015-A04[4].

Modeling[edit]

The cause of the spikes at 1 Hz and harmonics is a tiny 1 second square wave embedded in signals from the affected channels. The method to estimate the 1 Hz signal is to build a template in the time domain synchronized with the spurious signal. The first step is dividing each second of data into time bins using OBT. The number of bins is computed using:

[math] nbins = fsamp * template\_resolution[/math]

where fsamp is the sampling frequency and is 136 Hz at 70 GHz, 80 at 44 GHz and 56 at 30 GHz. Then the bins vector is initialized with time intervals. To avoid aliasing effects the template resolution is [math] \sqrt {3} [/math] . We can write the process adding an indices to the time sample: the lower index denotes the particular time sample, while the upper index labels the bin into which the sample falls. The linear filter can be written as:

[math] s(t_{i}^{j}) = a_j \left(1- \Delta x (t_{i}^{j}) \right) + a_{j+1} \Delta x (t_{i}^{j})[/math]

Here [math] \Delta x (t_{i}^{j})[/math] is the filter weight which is determined by where within the bin the sample lies. If we use [math] t^j [/math] with only an upper index to denote the start of each bin, then we can write the filter weight as follows:

[math] \Delta x (t_{i}^{j}) = {{{t_i^j - t^j} \over {t^{j+1} - t^j}}} [/math]

In other words, the filter weight is the time sample value minus the start of the bin divided by the width of the bin.

We must estimate the parameters [math] a_j [/math] from the data. With the assumption that the instrument has stable noise properties, we can use a least square algorithm to estimate the bin values:

[math] {\partial \over \partial a_k} \sum_{i,j} \left( s(t_i^j) – d_i^j \right)^2 = 0 [/math]

This can be represented in matrix equation:

[math] M_{jk}a_k = b_j [/math]

with the following definitions:

[math] M_{k,k-1} = \sum_i (1 - \Delta x (t_i^{k-1})) \Delta x (t_i^{k-1}) [/math]
[math] M_{k,k} = \sum_i (1 - \Delta x (t_i^k))^2 \Delta x (t_i^{k-1})^2 [/math]
[math] M_{k,k+1} = \sum_i (1 - \Delta x (t_i^k)) \Delta x (t_i^k) [/math]
[math] M_{k,k+n} (|n| \gt 1) = 0 [/math]
[math] b_j = \sum_i d_i^k (1- \Delta x (t_k^i)) + d_i^{k-1}\Delta x (t_i^{k-1}) [/math]

With these definitions we have to make use of periodic boundary conditions to obtain the correct results, such that if [math] k = 0 [/math], [math] k-1 = n-1 [/math] and [math] k = n-1 [/math], [math] k+1 = 0 [/math]. Once this is done, we have a symmetric tridiagonal matrix with additional values at the upper right and lower left corners of the matrix. The matrix is solved with LU decomposition. In order to be certain of the numerical accuracy of the result, we can perform a simple iteration. The solving of the linear system and the iterative improvement of the solution are implemented as suggested in Numerical Recipes.

Application[edit]

For each of the 44 LFI diodes there is the corresponding object in the Database. After studying the amplitude of the spikes at other LFI frequencies, we chose to apply the correction only to the 44 GHz radiometers. Each object contains 3 columns: the bins start time vector, the sky amplitudes and the reference amplitudes.

For each sample the value to be subtracted is computed using:

[math] V = skyAmp_k (1 - \Delta x (t_k)) + skyAmp_{k+1} \Delta x (t_k) [/math]

where k is the index of the bins at a given time.

Filling Gaps in the Data[edit]

During the mission, a small number of data packets were lost (Planck-2013-II[1], Planck-2015-A03[3]). Moreover in two different and very peculiar situations, LFI was shut down and restarted, giving inconsistencies in data sampling. None of these data are used for scientific purpose but to avoid discrepancies in data analysis all of the radiometers at the same frequency must have the same samples.

To ensure this, we compare the length of the data stream to be reduced in a specific pointing period to the data stored in the sample information object. If the length is not the same, the OBT vector is filled with missing sample times, the data vector is filled with zeros, and in the flag column the bit for gap is entered.

Gain Modulation Factor[edit]

The pseudo-correlation design of the LFI radiometers allows a dramatic reduction of [math] 1/f [/math] noise when the [math] V_{sky} [/math] and [math] V_{load} [/math] outputs are differenced (see LFI instrument description). The two streams are slightly unbalanced, as one looks at the 2.7 K sky and the other looks at the ~4.5 K reference load. To force the mean of the difference to zero, the load signal is multiplied by the Gain Modulation Factor (R). For each pointing period this factor is computed using (see eq. (3) in LFI description):

[math] R = { \lt V_{sky}\gt \over \lt V_{load}\gt } [/math]

Then the data are differenced using:

[math] TOI_{diff}^i = V_{sky}^i – R V_{load}^i [/math]

This value for R minimizes both the 1/f noise and the white noise in the difference timestream. The i index represents the diode and can be 0 or 1.

At this point, we also set the maneuver flag bit to identify which samples have missing data, using the information stored in the sampling information object. This identifies which data to ignore in the next step of the Pipeline.

The R values are stored in the database. At the same time the mean values of [math] V_{sky} [/math] and [math] V_{load} [/math] are stored so they can be used in other steps of the analysis.

Diode Combination[edit]

The two complementary diodes of each radiometer are combined. The relative weights of the diodes in the combination are chosen for optimal noise. We assign relative weights to the uncalibrated diode streams based on their first order calibrated noise.

Evaluation[edit]

From first order calibration we compute an absolute gain [math] G_0 [/math] and [math] G_1 [/math], subtract an estimated sky signal and calculate the calibrated white noise [math] \sigma_0 [/math] and [math] \sigma_1 [/math], for the pair of diodes (see eq. (6) in LFI description):. The weights for the two diodes ([math] i [/math] = 0 or 1) are:

[math] W_i = {\sigma_i^2 \over G_{01}} {1 \over {\sigma_0^2 + \sigma_1^2}} [/math]

where the weighted calibration constant is given by:

[math] G_{01} = {1 \over {\sigma_0^2 + \sigma_1^2}} [G_0 \sigma_1^2 + G_1 \sigma_0^2] [/math]

The weights are fixed to a single value per diode for the entire dataset. Small variations in the relative noise of the diodes would in principle suggest recalculating the weights on shorter timescales, however, we decided a time varying weight could possibly induce more significant subtle systematics, so chose a single best estimate for the weights for each diode pair.

LFI weights.jpg

Weights used in combining diodes


Application[edit]

The weights in the table above are used in the formula:

[math] TOI_{diff} = w_0 TOI_{diff0} + w_1 TOI_{diff1} [/math]

Planet Flagging[edit]

Extraction Method[edit]

Measurements of planets have been formed from samples containing flux from the planet, plus a surrounding region, projected onto a grid (microstripes), by assuming an elliptical Gaussian beam using parameters from instrument database.

Microstripes are a way to extract and store relevant samples for planet detection. Relevant samples are samples affected by the planet plus samples in the neighborhood (to establish a background level). The search radius to select samples as relevant is 5 deg around the planet's position, computed at the pointing period mid time. For each sample we store SCET (Spacecraft Event Time), pointing directions and calibrated temperature. Destriping is applied.

Random errors are estimated by taking the variance of samples entering each micromap pixel. This is fast but not exact, since the variance is larger near a bright source. This can cause the noise to be overestimated by a factor of. Given the large S/N for planetary observations, that is not a major drawback.

The apparent position of a planet as seen from Planck at a given time is derived from JPL Horizon. Positions are tabulated in steps of 15 minutes and then linearly interpolated at the sampling frequency of each detector. JPL Horizons tables allow also to derive other quantities such as the planet-Planck distance and the planet-Sun distance and the planet angular diameter, which affects the apparent brightness of the planet.

The antenna temperature is a function of the dilution factor, according to:

[math] T_{ant,obs} = 4 log 2T_{ant,1} \left( {\theta \over b_{fwhm} } \right) ^2 [/math]

where [math] T_{ant,obs} [/math] and [math] T_{ant,1} [/math] are the observed and reduced [math] T_{ant} [/math], [math] \theta [/math] the instantaneous angular diameter of the planet and [math] b_{fwhm} [/math] the beam full width at half maximum.

With the above definition [math] T_{ant,1} [/math] could be considered as the [math] T_{ant} [/math] for a planet with [math] b_{fwhm} = \theta [/math], but a more convenient view is to take a reference dilution factor [math] D_0 [/math], as the dilution factor for a standardized angular diameter for the planet and fiducial beam fwhm [math] b_{fwhm} [/math], [math] \theta_0 [/math], to have:

[math] D_0 = \left( {\theta_0 \over b_0 } \right) ^2 [/math]

leading to the following definition of a standardized [math] T_{ant} [/math]:

[math] T_{ant,obs} = 4 log 2 T_{ant,0} \left( {b_{fwhm,0} \over b_{fwhm}} {\theta \over \theta_0} \right) ^2 [/math]

with the advantage of removing variations among different detectors and transits while keeping the value of [math] T_{ant} [/math] similar to that seen by the instrument and then allowing a prompt comparison of signals and sensitivities.

Application[edit]

The OBT vectors found by the search are saved in a set of objects, one for each horn. In the Level 2 pipeline those OBTs are compared with the OBT vector of the data to set the planet bit flag where needed.

Photometric Calibration[edit]

Photometric calibration is the procedure used to convert data from volts to kelvin. The source of the calibration is the well known CMB dipole, caused by the motion of the Solar System with respect to the CMB reference frame. To this signal we add the modulation induced by the orbital motion of Planck around the Sun. The resulting signal is then convoluted with the horn beam to get the observed dipole. For details refer to Planck-2015-A06[5].

Beam Convolved Dipole[edit]

In computing the beam convolved dipole we used an elegant algorithm to save time and computing power. In computing the cosmological dipole signal it is common to assume a pencil-like beam acting as a Dirac delta function. In this case a dipole timeline is defined as:

[math] \Delta T_{D,\delta}(t) = \mathbf{P}_E(t) \cdot \mathbf{D}_E [/math]

where [math] \mathbf{P}_E(t) [/math] is the pointing direction, in the observer reference frame and [math] \mathbf{D}_E [/math] is the dipole axis scaled by the dipole amplitude again in the same reference frame.

In general the true signal would have to be convolved with the beam pattern of the given radiometer, usually described as a fixed map in the beam reference frame or as a time dependent map in the observer reference frame. In this case it is easiest to describe the convolution in the beam reference frame, since the function to be convolved then can be described by a single vector.

Denoting with [math] \mathcal{U}(t) [/math] the matrix converting from the observer to the beam reference frame, so that:

[math] \mathcal{U}(t) \mathbf{P}_E(t) = \mathbf{e}_z [/math]

the instantaneous dipole direction in the beam reference frame is:

[math] \mathbf{D}(t) = \mathcal{U}(t) \mathbf{D}_E [/math]

By denoting with [math] \mathbf{P} [/math] a pointing direction in the beam reference frame then:

[math] \Delta T_D(t) = N \int_{4\pi} B(\mathbf{P})\mathbf{P} \cdot \mathbf{D}(t) d^3\mathbf{P} [/math]

where [math] \mathbf{N} [/math] is a normalization constant.

[math] N^{-1} = \int_{4\pi} B(\mathbf{P}) d^3\mathbf{P} [/math]

Denoting with [math] \mathbf{P}_x [/math], [math] \mathbf{P}_y [/math], [math] \mathbf{P}_z [/math] the three cartesian components of the [math] \mathbf{P} [/math] the integral of the dot product can be decomposed into three independent integrals:

[math] S_x = N \int_{4\pi} B(\mathbf{P}) P_x d^3\mathbf{P} [/math]
[math] S_y = N \int_{4\pi} B(\mathbf{P}) P_y d^3\mathbf{P} [/math]
[math] S_z = N \int_{4\pi} B(\mathbf{P}) P_z d^3\mathbf{P} [/math]

those integrals define a time independent vector characteristic of each radiometer and constant over the mission.

Detector ID [math] S_x [/math] [math] S_y [/math] [math] S_z [/math]
LFI18S 4.6465302801434851e-03 -1.3271241134658798e-03 9.9577083251002374e-01
LFI18M 2.8444949872572308e-03 -8.9511701920330633e-04 9.9735308499993403e-01
LFI19S 4.2948823053610341e-03 -1.4037771886936600e-03 9.9612844031547154e-01
LFI19M 4.5097043389558588e-03 -1.6255119872940569e-03 9.9573258887316529e-01
LFI20S 5.2797843121430788e-03 -1.9479348764780897e-03 9.9514557074549859e-01
LFI20M 4.7420554683343680e-03 -1.8462064956921880e-03 9.9548302435765323e-01
LFI21S 5.2377151034493320e-03 1.9337533884005056e-03 9.9517552967578515e-01
LFI21M 4.3713123245201109e-03 1.7075534076511198e-03 9.9585142295482576e-01
LFI22S 3.6272964935149194e-03 1.2030860597901591e-03 9.9679007501368400e-01
LFI22M 3.1817723013115090e-03 1.1383448292256770e-03 9.9705527860728405e-01
LFI23S 3.2582522957523585e-03 9.1080217792813213e-04 9.9706066300815721e-01
LFI23M 2.6384956780749246e-03 8.0545189089709709e-04 9.9755224385088148e-01
LFI24S 1.2007669925034473e-03 -2.9877396130232452e-07 9.9897924004634964e-01
LFI24M 1.1966952022302408e-03 -1.3590311231894260e-06 9.9894545519286426e-01
LFI25S -2.1006042848582498e-04 4.4287041964311393e-04 9.9965683952381934e-01
LFI25M -2.5133433631240997e-04 5.5803053337715499e-04 9.9954614333864866e-01
LFI26S -1.9766749938836072e-04 -4.2202164851420096e-04 9.9966955259454382e-01
LFI26M -2.5426738260127387e-04 -5.7426678736881309e-04 9.9952104084618332e-01
LFI27S 5.8555501689062294e-03 1.8175261751324087e-03 9.9448302844644199e-01
LFI27M 4.9962842224387377e-03 1.5766304210447924e-03 9.9529240281463061e-01
LFI28S 6.3750745068891614e-03 -1.9617207193203534e-03 9.9396633709784610e-01
LFI28M 4.7877703093970429e-03 -1.5442874049905993e-03 9.9551694709996730e-01


By using this characteristic vector the calculation of the convolved dipole is simply defined by a dot product of the vector [math] \mathbf{S} [/math] and the dipole axis rotated in the beam reference frame.

[math] \Delta T (t) = \mathbf{S}^T \mathcal{U}(t) \mathbf{D}_E [/math]

Binning[edit]

In order to simplify the computation and to reduce the amount of data used in the calibration procedure the data are phase binned in map with [math]N_{side}[/math] 256. During phase binning all the data flagged for maneuvers, planets and gaps, as well as the ones flagged in Level 1 analysis as not recoverable, are ignored.

Fit[edit]

The first order calibration values are given by a least squares fit between the signal and the dipole. For each pointing gain ([math] g_k [/math]) and offset ([math] b_k [/math]) values are computed by minimizing:

[math] \chi^2 = \sum_{i \in k} {[ \Delta V (t_i) - \Delta V_m (t_i|g_k, b_k)]^2 \over rms_i^2 } [/math]

The sum includes only samples outside a Galactic mask.

DaCapo[edit]

The largest source of error in the fit arises from unmodeled sky signal [math] \Delta T_a [/math] from CMB anisotropy.

The procedure adopted to correct this effect is described below.

The uncalibrated toi [math]y_{i}[/math] belonging to a pointing period [math]k[/math] is modeled as

[math]y_{i}=g_{k}*(s_{i}+d_{i})+b_{k}+n[/math]

where the unknowns [math]s_{i}[/math],[math]g_{k}[/math] and [math]b_{k}[/math] are respectively the signal per OBT, the gain per PID and the offset per PID, and [math]d_{i}[/math] is the total dipole.

This is a non linear [math]\chi^{2}[/math] minimization problem. It can be linearized with an iterative procedure: for each iteration step a linearized model toi is built as

[math]y_{i,model}=g_{k}(-1)*(s_{i}-s_{i}(-1))+g_{k}*(d_{i}+a_{i}(-1))+b_{k}[/math]

where [math]g_{k}(-1)[/math] and [math]s_{i}(-1)[/math] are the gain and sky signal computed at the previous step of iteration.

A linear fit is performed between [math]y_{i}[/math] and [math]y_{i,model}[/math] to get the new [math]g_{k}[/math] and [math]s_{i}[/math] . The procedure is iterated until convergence.

To reduce the impact of the noise during the iterative procedure the sky estimation is built using data from both radiometers of the same horn.

Smoothing[edit]

To improve accuracy given by the iterative algorithm and remove noise from the solution a smoothing algorithm must be performed.

OSGTV[edit]

OSGTV is a 3 step smoothing algorithm, implemented with a C++ code.

The gain reconstructed with DaCapo can be expressed as

[math]G_{raw}(t,T(t))=G_{true}(t,T(t))+G_{noise}(t)[/math]

where the [math]T(t)[/math] function represents temperature fluctuations of the Focal Plane. The time dependence of the “real” gain [math]G_{true}[/math] is modeled as the superposition of a “slow” component, with a time scale of ~3 months, and a “fast” component with a time scale of few PIDs:

[math]G_{true}=G_{slow}+G_{fast}[/math].

The slow component takes into account the seasonal variations of the thermal structure of the spacecraft due to the orbital motion, while the “fast” component describes the thermal effects of the electronics and compressors, as well as single events like the sorption cooler switchover.

To disentangle the components of [math]G_{true}[/math], we need a parametric “hybrid” approach in three steps.

Step 1. For each PID [math]i[/math], we used the 4K total-power [math]V_{load}(i)[/math] and the signal from temperature detectors in the focal plane [math]T_{i}[/math], subsampled at 1 sample/PID, to track gain changes. This is implemented through a linear fit between [math]V_{load}(i) * G_{raw}[/math] and [math]T_{i}[/math]:

[math]V_{load}(i)*G_{raw}=K_{0}+K_{1}*(T_{i}-\langle T_{MW}\rangle)[/math]

where the window length [math]MW[/math] of the moving average is proportional to the variance of the dipole in the considered PID. The resulting gain is:

[math] G_{i} = \frac{ K_{0} + K_{1} * (T_{i} - \langle T_{MW} \rangle ) }{ V_{load}(i) } [/math]

Step 2. The "fast" component [math]G_{fast}[/math] is recovered as follows: we define a maximum moving average window length [math]WL_{max}[/math],and for each PID we compute the variance [math]\sigma_{i}^2[/math] of [math]G_{i}[/math] in this window. We define a percentile on the ordered variance array and we compute the corresponding value of the variance [math]\sigma_{p}^2[/math]. The window length for each PID is then computed as

[math]WL_{i}=\frac{\sigma_{p}^2*WL_{max}}{\sigma_{i}^2}[/math].

If [math]WL_{i}\gt WL_{max}[/math] we impose [math]WL_{i}=WL_{max}[/math]. With these window lengths a moving window average is performed on [math]G_{i}[/math]. The averaged gain vector is subtracted from the raw hybrid gain to get [math]G_{fast}[/math].

Step 3. We perform a moving window average on [math]G_{i}[/math] and we compute the variances of the smoothed array. The window length is computed with a linear interpolation between a minimum length [math]WL_{min}[/math] defined in the dipole minima and a maximum [math]WL_{max}[/math] defined in the dipole maxima. The array of gain variances is weighted with the variance of the dipole. We set a percentile on the variance array and we find the corresponding variance value [math]\sigma_{p}^2[/math]. Around this value we search for local maxima of the variance array, and we split the domain of the gain in subsets between consecutive maxima. For each subset we perform a moving average with the corresponding window length. The "slow" component [math]G_{slow}[/math] is given by the union of these subsets.

Gain Application[edit]

The last step in TOI processing is the creation of the calibrated stream. For each sample we have:

[math] TOI_{cal}(t) = (TOI_{diff}(t) – offset(k)) \cdot g(k) \, - \, convDip(t)\, - \, T_{stray}(t)[/math]

where t is the time and k is the pointing period, [math] convDip [/math] is the CMB Dipole convolved with the beam, and [math]T_{stray}[/math] is the straylight.

Noise[edit]

This pipeline step aims at the reconstruction of the noise parameters from calibrated flight TOI. The goal is two-fold. On the one hand we need to know the actual noise properties of our instrument in order to properly take them into account, especially during later processing and analysis steps like map-making and power spectrum estimation. On the other hand evaluation of noise properties during the instrument life-time is a way to track down possible variations, anomalies and general deviations from the expected behaviour.

Operations[edit]

Noise estimation is performed on calibrated data; since we would like to track possible noise variations during the mission life-time, we select data in chunks of 5 ODs (Operational Days). These data are processed by the ROMA Iterative Generalized Least Square (IGLS) map-making algorithm which includes a noise estimation tool. In general IGLS map-making is a quite costly in terms of time and resources required. However the length of the data is such that it can run on the DPC cluster in very short time (~1-2 minutes).

The method implemented can be summerized as follows. We model the calibrated TOI as

[math] \mathbf{\Delta T} = \mathbf{P} \mathbf{m} + \mathbf{n} [/math]

where [math] \mathbf{n} [/math] is the noise vector and [math] \mathbf{P} [/math] is the pointing matrix that links a pixel in the map [math] \mathbf{m} [/math] with a sample in the TOI [math] \mathbf{d} [/math]. The zero-th order estimation of the signal is obtained by simply rebinning TOI into a map. Then an iterative approach follows in which both signal and noise are estimated according to

[math] \mathbf{\hat{n}_i} = \mathbf{\Delta T} - \mathbf{P\hat{m}_i} [/math]
[math] \mathbf{\hat{m}_{i+1}} = \mathbf{(P^T\hat{N}_i^{-1}P)^{-1}P^T\hat{N}_i^{-1}\Delta T} [/math]

where [math] \mathbf{\hat{N}_i} [/math] is the noise covariance matrix in the time domain resulting from iteration [math]i [/math]. After three iterations, convergence is achieved.

We then perform an FFT (Fast Fourier Transform) on the noise time stream from the iterative approach and fit the resulting spectrum.

Fitting Pipeline[edit]

As already done in the 2013 release, we estimate the parameters of the noise properties (i.e. white noise level, knee-frequency and slope of the low-frequency noise component) by means of a MCMC approach. Therefore on the spectra just described we first compute the white noise level taking the mean of the last 10% of data (at 30 GHz due to the higher value of the knee-frequency this quantity is reduced to 5%). Once the white noise level has been determined we proceed with the actual fitting of knee-frequency and slope. The resulting values reported are the medians of the fitted values for our 5 ODs chunk along the whole mission lifetime.

Horn White Noise M [[math]\mu[/math]K s[math]^{1/2}[/math]] White Noise S [[math]\mu[/math]K s[math]^{1/2}[/math]] Knee-frequency M [mHz] Knee-frequency S [mHz] Slope M Slope S
18 513.0[math]\pm[/math]2.1 467.2[math]\pm[/math]2.3 14.8[math]\pm[/math]2.5 17.8[math]\pm[/math]1.5 -1.06[math]\pm[/math]0.10 -1.18[math]\pm[/math]0.13
19 579.6[math]\pm[/math]2.2 555.0[math]\pm[/math]2.2 11.7[math]\pm[/math]1.2 13.7[math]\pm[/math]1.3 -1.21[math]\pm[/math]0.26 -1.11[math]\pm[/math]0.14
20 587.3[math]\pm[/math]2.1 620.5[math]\pm[/math]2.7 8.0[math]\pm[/math]1.9 5.7[math]\pm[/math]1.5 -1.20[math]\pm[/math]0.36 -1.30[math]\pm[/math]0.41
21 451.0[math]\pm[/math]1.7 560.1[math]\pm[/math]2.0 37.9[math]\pm[/math]5.2 13.3[math]\pm[/math]1.5 -1.25[math]\pm[/math]0.09 -1.21[math]\pm[/math]0.09
22 490.8[math]\pm[/math]1.5 531.3[math]\pm[/math]2.3 9.7[math]\pm[/math]2.3 14.8[math]\pm[/math]6.7 -1.42[math]\pm[/math]0.23 -1.24[math]\pm[/math]0.30
23 504.3[math]\pm[/math]1.8 539.7[math]\pm[/math]1.8 29.7[math]\pm[/math]1.1 59.0[math]\pm[/math]1.4 -1.07[math]\pm[/math]0.03 -1.21[math]\pm[/math]0.02
24 463.0[math]\pm[/math]1.4 400.7[math]\pm[/math]1.3 26.8[math]\pm[/math]1.3 88.3[math]\pm[/math]8.9 -0.94[math]\pm[/math]0.01 -0.91[math]\pm[/math]0.01
25 415.3[math]\pm[/math]1.5 395.4[math]\pm[/math]2.9 20.1[math]\pm[/math]0.7 46.4[math]\pm[/math]1.8 -0.85[math]\pm[/math]0.01 -0.90[math]\pm[/math]0.01
26 483.0[math]\pm[/math]1.9 423.2[math]\pm[/math]2.5 64.4[math]\pm[/math]1.9 68.2[math]\pm[/math]9.5 -0.92[math]\pm[/math]0.01 -0.76[math]\pm[/math]0.01
27 281.5[math]\pm[/math]2.1 303.2[math]\pm[/math]1.8 174.5[math]\pm[/math]2.9 108.8[math]\pm[/math]2.5 -0.93[math]\pm[/math]0.01 -0.91[math]\pm[/math]0.01
28 317.1[math]\pm[/math]2.4 286.5[math]\pm[/math]2.3 130.1[math]\pm[/math]4.4 43.1[math]\pm[/math]2.4 -0.93[math]\pm[/math]0.01 -0.90[math]\pm[/math]0.02

Time variations of noise parameters are a good tracer of possible modifications in the instrument behaviour and we know that some events capable of affecting instrument behaviour had happened during the mission. Both variations of the physical temperature of the instrument due to the transition in the operations from the first sorption cooler to the second one, as well as the observed degradation in the performances of the first cooler are events that clearly show their fingerprints in the variation of the noise spectra.

In the following figure we report a representative sample of noise spectra, one for each frequency channel (from left to right LFi 18M, LFI 25S and LFI 27M), covering the whole mission lifetime. The white noise is extremely stable at the level of 0.3%. As already noted in the 2013 release, knee-frequencies and slopes are stable until OD 326 and show clear variations and deviations from the simple one knee-frequency one slope model. This is a sign of the degradation of the first cooler inducing thermal variations with a characteristic knee-frequency (different from the radiometric one) and with a steeper slope. Once the second cooler became operational and performed as expected, the noise spectra gradually returned to the original shape at the beginning of the mission. This behaviour is visible, although not identical at each frequency for several reasons e.g. intrinsic thermal susceptibility and position on the focal plane that determines the actual thermal transfer function.


'"'LFI 18M.'"' '"'LFI 25S.'"' '"'LFI 27M.'"'

References[edit]

  1. 1.01.11.2 Planck 2013 results. II. Low Frequency Instrument data processing, Planck Collaboration, 2014, A&A, 571, A2
  2. 2.02.1 Planck 2013 results. III. Low Frequency Instrument systematic uncertainties, Planck Collaboration, 2014, A&A, 571, A3
  3. 3.03.13.2 Planck 2015 results. II. LFI processing, Planck Collaboration, 2016, A&A, 594, A2.
  4. 4.04.1 Planck 2015 results. III. LFI systematics, Planck Collaboration, 2016, A&A, 594, A3.
  5. Planck 2015 results. V. LFI calibration, Planck Collaboration, 2016, A&A, 594, A5.

(Planck) Low Frequency Instrument

Operation Day definition is geometric visibility driven as it runs from the start of a DTCP (satellite Acquisition Of Signal) to the start of the next DTCP. Given the different ground stations and spacecraft will takes which station for how long, the OD duration varies but it is basically once a day.

Data Processing Center

Planck Science Office

On-Board Time

analog to digital converter

[LFI meaning]: absolute calibration refers to the 0th order calibration for each channel, 1 single number, while the relative calibration refers to the component of the calibration that varies pointing period by pointing period.

LFI Data Acquisition Electronics

Cosmic Microwave background