Revision 12010-10-04 - BartVandenbussche

Line: 1 to 1
>
>
 META TOPICPARENT name="PacsDocumentation"
```

======
Notes:
------

HISTORY:
30-Apr-2009: Dave Frayer    Initial write-up
14-Oct-2009: Babar Ali      Updated API change
12-Jul-2010: Cate Liu       Update API documentation
16-Jul-2010: Cate Liu       Added new parameter to makeTodArray

======================================
--------------------------------------

=================
Summary of Method
-----------------

MADmap uses a maximum-likelihood technique to build a map from an
input Time Order Data (TOD) set by solving a system of linear
equations.  It is used to remove low-frequency drift ("1/f") noise
from bolometer data while preserving the sky signal on large spatial
scales.

The TOD d = d(t) for t=1...n. is given by

(1) d = As+n,

where A is the "pointing" matrix which relates each time sample to the
corresponding pixel in the sky; A(t,p) for time t=1...n and
sky-pixels, p=1...m, s is the sky map s(p), for pixels p=1...m, and n
is the noise n(t).

The time-time noise covariance matrix N is given by

N = ,

where "^T" represents transpose and "< >" represents the ensemble
average.  In practice, the noise is assumed to be Gaussian,
"piece-wise stationary", and "circulant" which enables the noise
properties to be described by its Fourier power spectrum.  Under these
assumptions, the noise depends only on the time difference between two
samples, and N is defined by its first row, N_i0, which is given by
the autocorrelation of the noise.

(2) N_i0 =  =

where "*" represents convolution and FFT is the Fourier transform.
Therefore, the time-time noise covariance matrix N is given by the
inverse Fourier transform of its Fourier power spectrum.

For Gaussian noise, the probability distribution [P(n)] is given by
(ignoring normalization factors):

(3) P(n) ~ exp[-0.5 n^T N^-1 n].

Using (1) and (3), the probability of sky (s) given data (d) is:

(4) P(s|d) ~ exp[-0.5 (d - As)^T N^-1 (d - As)].

Maximizing this probability yields A^T N^-1 As = A^T N^-1 d, which has
a solution:

(5) s = (A^T N^-1 A)^-1 A^T N^-1 d.

MADmap solves (5) for the sky-map (s) given input TOD (d).  MADcap
solves the the inverse time-time noise covariance matrix N^-1
(N_i0^-1) without requiring inversion under the assumptions for
equation (2) by

(6) N_i0^-1 = .

N^-1 effectively acts like a set of convolutions with band-limited
kernels.  N^-1 is the InvNtt calibration file used by MADmap and is
effectively the inverse FFT of 1/(power_spectrum).

References:

J.-C. Hamilton, 2003, C. R. Physique 4, 871

Ashdown, M. A. J., et al. 2007, A&A, 467, 761

D. Clements, et al. 2006, "SPIRE Mapmaking Algorithm Review Report"

==============
Implementation
--------------

The implementation is based on three modules:

(1) photOffsetCorr
(2) photModuleDriftCorrection
(3) photGlobalDriftCorrection
(4) makeTodArray

==============
PhotOffsetCorr
--------------

I. Summary

Removes detector-to-detector offsets before running MADmap.  Defaults
to a median of input frames cube, but supports removal of an offset
calibration file.

II. Usage

frames = PhotOffsetCorr(Frames frames, int copy, int zeroLevelCorr, Selection regionSelection)

frames -- Input frames, MANDATORY, NO default value
copy -- default 0: input frames are changed
1: new frames are copied and returned
zeroLevelCorr -- default 0: median removal
1: zero level calibration value subtracted
regionSelection -- default None: noOp; Otherwise, if zeroLevelCorr=0, do median
removal using the given region; if zeroLevelCorr=1, do zero
level calibration value subtraction using given region

Assumes BAND information is present for selection of proper
calibration file.  If BAND not present, checks dimensions for Red vs
Blue data and defaults to BS for for Blue.

III. Functionality

frames_OUT = frames_IN - OFFSET
where OFFSET is median per detector of frames_IN, or OFFSET is input calibration file.

=========================
photModuleDriftCorrection
-------------------------

=========================
photGlobalDriftCorrection
-------------------------

============
makeTodArray
------------

I. Summary

Builds time-ordered data (TOD) stream for input into MADmap and
derives meta header information of the output skymap. Input data is
assumed to be calibrated and flat-fielded.  Also prepares the "to's"
and "from's" header information for the InvNtt (inverse time-time
noise covariance matrix) calibration file.    Assumes BAND and
OnTarget-flag in the status.

II. Usage

PacsTodProduct todProd = makeTodArray(Frames frames, Double scale,
Double crota2, Boolean optimizedOrientation,
Double minRotation, Boolean chunkScanLegs, PacsCal calTree)

frames -- Data frames in units of mJy/pixel.  Required input
meta-data: (1) RA,Dec cubes associated with the frames
including the effects of distortion.  Assumes this step has
which identifies bad pixels, (3) information on band
(BS,BL,R), mode (scan/chopped raster), and locations between
scan legs for data "chunking".

scale -- Default = 1.0.  Multiplication pixel scale factor of output sky pixels
in relation to nominal PACS detector size (3.2" for blue and
6.4" red).  For scale = 1.0, the skymap has square pixels equal
to nominal PACS detector size; for scale = 0.5, then 1.6" for
blue and 3.2" for red.

crota2 -- Default = 0.0 degree.  CROTA2 of output skymap.

optimizeOrientation -- Default = false
when true: do auto map rotation if calculated auto angle > minRotation

minRotation -- Default = 15.0 degrees.  Minimial angle for auto rotation if
optimizeOrientation=true.

chunkScanLegs -- Default = true, on-target flags are used to chunk scan legs.
i.e. off-target data are not used.

calTree -- Default = null, Pacs calibration tree.

todProd -- Output product with TOD file name, final output map's Wcs, to/from index
for each good data stream segment and the correspondence between tod
indices and sky map pixels.

The intermediate output TOD file is saved in a directory specified by the
property var.hcss.workdir (please note that this property has default value
defined in HCSS).  It is removed after MADmap finishes and at the exit of HIPE.

The body of the TOD file is a byte stream binary data file consisting of

The binary header is four 8byte integers

(1) First sample index for TOD data, set to 0.

(2) Last sample index for TOD data chunk, set to
(n_good_detectors * n_samples) -1.

(3) nnObs = Number of detector values per sky pixel
during each time sample (for default one-to-one mapping of
detectors on to sky pixels, nnObs=1).

(4) total number of sky pixels with good data.

The binary header is followed by the data in the order of:

For each input GOOD detector pixel ("observation"):
value (double, 8-byte float) ==v
For each sky pixel observed:
weight (4-byte float)  == w
skypixel index (4-byte int) == p
(e.g.,.) for
good detectors ii=1,nd and time samples kk=1,nt, TOD order is given by:
for ii=1,nd
for kk=1,nt
v[ii,kk]
for jj=1,nnObs
w[ii,kk]
p[ii,kk]

Initially for the SPG, we will set nnOBS=1, i.e., use the default
one-to-one mapping of input detectors onto sky pixels.

The output tod product includes the astrometry of output map using Wcs, in
particular meta data keywords such as:

CRVAL1    RA Reference position of skymap
CRVAL2    Dec Reference position of skymap
EQUINOX   2000.
CTYPE1    RA---TAN
CTYPE2    DEC--TAN
CRPIX1    Pixel x value corresponding to CRVAL1
CRPIX2    Pixel y value corresponding to CRVAL2
CDELT1    pixel scale of sky map (=input as default, user parameter)
CDELT2    pixel scale of sky map (=input as default, user parameter)
CROTA2    PA of image N-axis (=0 as default, user parameter)

III. Functionality

(1) Build TOD binary data file with format given above.

(2) Define the astrometry of output map and save as keywords give
above.  Default CROTA2=0.0, but if optimizedOrientation=True, then
call Maptools to compute the optimial rotation (i.e., for
elongated maps).  If rotation less than minRotation, then leave
map un-rotated with CROTA2=0.

detectors will not have good InvNtt data and should be discarded

(4) SKYPIX INDICES --  Compute the skypix indices which are derived from
the projection of each input pixel onto the output sky grid.  The
skypix indices are increasing integers representing the location
in the sky map with good data.  The skypixel indices of the output
map must have some data with non-zero weights, must be continuous,
must start with 0, and must be sorted with 0 first and the largest
index last.

(5) GLITCHES -- Set weights=BADwtval for bad data as flagged in the
Default BADwtval =0.0, but may need to use a small non-zero value
(e.g., 0.001) in practice (to avoid MADmap precondition that
requires non-zero wts and data for each observed skypixel).  Good
data should have weights set to 1.0 for initial versions with
nnOBS=1.

(6) CHUNK per SCAN LEG -- Use the OnTarget status flag to define the
boundaries for data chunking per scan leg.  The sample start and
end boundaries of the TOD indices of each data chunk per detector
are needed for the InvNtt headers (the "tos" and "froms").  TOD
rules:
(1) Chunk for large GAPS (>maxGap) defined by OnTarget (it is
assumed that the scan turn-around time will be larger than
maxGap, but this is not a requirement).
(2) for small GAPS (<=maxGap) defined by OnTarget use the data
values for the TOD, but set the TOD weights=BADwtval (i.e.,
data not effectively used for map, but needed for continuous
noise estimate for MADmap).   This condition is not expected
for properly handled data products upstream, but could exists
if there are issues with the pointing products.
(3) Include an option to turn-off chunking by scan leg.

(7) CHUNK as function of TIME per DETECTOR based on the mask
information -- Check TOD stream per detector.  For "gaps" of bad
data in samples larger than maxGap, then chunked the data.  Ignore
data streams that have a number of samples less than minGOOD,
i.e., each TOD chuck should be larger or equal to minGOOD samples.
For gaps, smaller or equal to maxGap, linearly interpolate the TOD
values across the gap and set TOD weights=BADwtval (i.e., data not
used for map, but needed for continuous noise estimate for
MADmap).  TOD rules (in this order):
(1) throw out initial and end samples that are bad.
(3) chuck for large bad GAPS (>maxGap)
(4) throw out small chucks (maxGap.  For GAPS<=maxGap, linearly interpolate
data values across the gaps for (7) [i.e., glitches], but use the data
values for (6) [i.e., nothing wrong with the data values].

=========
---------

I. Summary

map making routine for converting time ordered data sets into
pixelized sky maps.  Input TOD data is assumed to be calibrated and
flat-fielded and input InvNtt is from calibration tree.

II. Usage

SimpleImage map = runMadMap(PacsTodProduct tod, PacsCal calTree,
Double maxerror, Integer maxiterations,
Boolean runNaiveMapper, Boolean useAvgInvntt)

Inputs:

tod -- PacsTodProduct from makeTodArray, MANDATORY.

calTree -- Pacs calibration tree containing calibration InvNtt
information stored as an array of size max(n_correlation+1) x
n_all_detectors.  Each row represents the InvNtt
information for each detector.  MANDATORY.

maxerror -- Default = 1e-5, maximum relative error allowed in PCG routine.

maxiterations -- Default = 200, maximum number of iterations in PCG routine.

runNaiveMapper -- Default = false, run MadMapper; when true, run NaiveMapper

useAvgInvntt -- Default = false, use per-pixel invntt data; when true,
use averaged invntt data

Output:
map -- Output product consisting of following:

(1) image -- Sky map image (either a Naive map or MADmap) with Wcs information

(2) coverage -- Coverage map corresponding to sky map, with Wcs

(3) error -- Uncertainty image associated with the map, with Wcs

The error map is currently only made properly for NaiveMapper.
Two calls are needed to produce both the MadMap and
NaiveMap simple image products (runNaiveMapper=true yields Naivemap
product and runNaivMapper=false yields MadMap product).

III.  Functionality

(1) Build InvNtt from input calibration tree.

Format of InvNtt chunk:

min_sample        starting sample index
max_sample        last sample index
n_correlation  correlation width of matrix
Data-DOUBLE:
invntt(n_correlation+1)

The min/max samples are the "tos" and "froms" calculated from a method
within makeTodArray.  The sample indices need to be consistent between
the TOD chunked files and the InvNtt file.