Source code for sofia_redux.instruments.hawc.steps.stepnodpolsub
# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Nod subtraction pipeline step."""
from astropy import log
import numpy as np
from sofia_redux.instruments.hawc.datafits import DataFits
from sofia_redux.instruments.hawc.stepparent import StepParent
__all__ = ['StepNodPolSub']
[docs]
class StepNodPolSub(StepParent):
r"""
Subtract low nods from high nods.
This step identifies low nods in the input data, and subtracts
them from the high nods. The data must have been taken in either
an ABBA pattern, or else must consist entirely of A nods (no
nodding occurred). Error images for all nods are propagated
appropriately.
For diagnostic purposes, the R and T data are also added across
all nod positions, and stored as a 'CMSIG' extension.
The input for this step is assumed to be generated by
`sofia_redux.instruments.hawc.steps.StepCombine`. The output is a
DataFits with several image and table extensions. Seven image extensions
for each HWP angle *M* are created: DATA R-T HWP M, ERROR R-T HWP M,
DATA R+T HWP M, ERROR R+T HWP M, VAR R HWP M, VAR T HWP M,
and CMSIG R+T HWP M. Also, a table, TABLE HWP M is created. Lastly,
BAD PIXEL MASK is copied from the input into the output file.
Notes
-----
For ABBA data, the output images (:math:`d'`) and associated errors
(:math:`\sigma'`) are calculated as:
.. math:: d' = \frac{1}{2}(d_A - d_B)
.. math:: \sigma' = \frac{1}{2}\sqrt{\sigma_A^2 + \sigma_B^2}
The variances in R and T are also propagated as:
.. math:: \sigma_R'^2 = \sigma_{R,A}'^2 + \sigma_{R,B}'^2
.. math:: \sigma_T'^2 = \sigma_{T,A}'^2 + \sigma_{T,B}'^2
"""
[docs]
def setup(self):
"""
Set parameters and metadata for the pipeline step.
Output files have PRODTYPE = 'nodpolsub', and are named with
the step abbreviation 'NPS'.
No parameters are currently defined for this step.
"""
# Name of the pipeline reduction step
self.name = 'nodpolsub'
self.description = 'Subtract Nods'
# Shortcut for pipeline reduction step and identifier for
# saved file names.
self.procname = 'nps'
# Clear Parameter list
self.paramlist = []
# Append parameters
# (none at this time)
[docs]
def run(self):
"""
Run the data reduction algorithm.
Because this step is single-in, single-out (SISO),
self.datain must be a DataFits object. The output
is also a DataFits object, stored in self.dataout.
The process is:
1. For each HWP angle, subtract the B nod from the A nod,
and divide by two (assuming a nod-match-chop pattern);
propagate errors accordingly.
2. Store the data in a DataFits object.
"""
self.dataout = DataFits(config=self.datain.config)
self.dataout.filename = self.datain.filename
self.dataout.setheader(self.datain.header)
nhwp = self.datain.getheadval('nhwp')
nodpatt = self.datain.getheadval("nodpatt")
nnod = len(set(nodpatt))
if nnod == 1:
log.info('Only 1 nod position in NODPATT={}'.format(nodpatt))
log.info('No nod subtraction performed; '
'only NOD 0 is propagated.')
# lists to hold the Left, Right data and sigma for a single HWP angle
for hwp in range(nhwp):
if nhwp > 1:
# R-T
rmt_data = [0.] * nnod
rmt_sigma = [0.] * nnod
else:
rmt_data = []
rmt_sigma = []
# R+T
rpt_data = [0.] * nnod
rpt_cmsig = [0.] * nnod
rpt_sigma = [0.] * nnod
r_var = [0.] * nnod
t_var = [0.] * nnod
table = [0.] * nnod
for nod in range(nnod):
root = 'HWP%d NOD%d' % (hwp, nod)
log.debug("Starting %s" % root)
if nhwp > 1:
rmt_data[nod] = self.datain.imageget(
'DATA R-T %s' % root)
rmt_sigma[nod] = self.datain.imageget(
'ERROR R-T %s' % root)
rpt_data[nod] = self.datain.imageget('DATA R+T %s' % root)
rpt_sigma[nod] = self.datain.imageget('ERROR R+T %s' % root)
r_var[nod] = self.datain.imageget('VAR R %s' % root)
t_var[nod] = self.datain.imageget('VAR T %s' % root)
table[nod] = self.datain.tableget('TABLE %s' % root)
rpt_cmsig[nod] = rpt_data[nod].copy()
if nodpatt == 'ABBA':
if nhwp > 1:
rmt_data[0] = 0.5 * (rmt_data[0] - rmt_data[1])
rmt_sigma[0] = 0.5 * np.sqrt(rmt_sigma[0]**2
+ rmt_sigma[1]**2)
rpt_data[0] = 0.5 * (rpt_data[0] - rpt_data[1])
rpt_sigma[0] = 0.5 * np.sqrt(rpt_sigma[0]**2
+ rpt_sigma[1]**2)
r_var[0] = 0.25 * (r_var[0] + r_var[1])
t_var[0] = 0.25 * (t_var[0] + t_var[1])
rpt_cmsig[0] = 0.5 * (rpt_data[0] + rpt_data[1])
elif nnod == 1:
# only 1 nod position, so nothing to subtract
pass
else:
msg = 'Can only process data with ABBA nod pattern'
log.error(msg)
raise ValueError(msg)
if nhwp > 1:
self.dataout.imageset(rmt_data[0], 'DATA R-T HWP%d' % hwp)
self.dataout.imageset(rmt_sigma[0], 'ERROR R-T HWP%d' % hwp)
self.dataout.imageset(rpt_data[0], 'DATA R+T HWP%d' % hwp)
self.dataout.imageset(rpt_sigma[0], 'ERROR R+T HWP%d' % hwp)
self.dataout.imageset(r_var[0], 'VAR R HWP%d' % hwp)
self.dataout.imageset(t_var[0], 'VAR T HWP%d' % hwp)
# R+T_Nod0 + R+T_Nod1, for diagnostic use only
self.dataout.imageset(rpt_cmsig[0], 'CMSIG R+T HWP%d' % hwp)
tbhdu = self.dataout.tablemergetables(table)
self.dataout.tableset(tbhdu.data, 'TABLE HWP%d' % hwp,
tbhdu.header)
self.dataout.copydata(self.datain, 'BAD PIXEL MASK')