> > 
META TOPICPARENT 
name="PacsDocumentation" 
======
Notes:

HISTORY:
30Apr2009: Dave Frayer Initial writeup
14Oct2009: Babar Ali Updated API change
12Jul2010: Cate Liu Update API documentation
16Jul2010: Cate Liu Added new parameter to makeTodArray
======================================
PACS MADmap Users Manual Documentation

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

MADmap uses a maximumlikelihood 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 lowfrequency 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
skypixels, p=1...m, s is the sky map s(p), for pixels p=1...m, and n
is the noise n(t).
The timetime 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,
"piecewise 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 timetime 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(sd) ~ 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 skymap (s) given input TOD (d). MADcap
solves the the inverse timetime 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 bandlimited
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
(5) runMadMap
==============
PhotOffsetCorr

I. Summary
Removes detectortodetector 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

See MADmap User's Guide.
=========================
photGlobalDriftCorrection

See MADmap User's Guide.
============
makeTodArray

I. Summary
Builds timeordered data (TOD) stream for input into MADmap and
derives meta header information of the output skymap. Input data is
assumed to be calibrated and flatfielded. Also prepares the "to's"
and "from's" header information for the InvNtt (inverse timetime
noise covariance matrix) calibration file. Assumes BAND and
OnTargetflag 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
metadata: (1) RA,Dec cubes associated with the frames
including the effects of distortion. Assumes this step has
been previously done by PhotAssignRaDec, (2) input mask cube
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, ontarget flags are used to chunk scan legs.
i.e. offtarget 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
header information and TOD data
(Reference: http://crd.lbl.gov/~cmc/MADmap/doc/man/MADmap.html).
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 onetoone 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, 8byte float) ==v
For each sky pixel observed:
weight (4byte float) == w
skypixel index (4byte 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
onetoone 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
RADESYS ICRS
EQUINOX 2000.
CTYPE1 RATAN
CTYPE2 DECTAN
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 Naxis (=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 unrotated with CROTA2=0.
(3) BADPIXEL  For "BAD" dead/bad detector which are always bad,
then do not include the detector in TOD calculations Dead/bad
detectors will not have good InvNtt data and should be discarded
for MADmap.
(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 nonzero 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
masks (BADPIXEL, SATURATION, GLITCH, UNCLEANCHOP). For (BLIND)
Default BADwtval =0.0, but may need to use a small nonzero value
(e.g., 0.001) in practice (to avoid MADmap precondition that
requires nonzero 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 turnaround 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 turnoff 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.
(2) fill in the small bad GAPS (<=maxGap), weights=BADwtval
(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].
=========
runMadMap

I. Summary
Wrapper that runs JAVA MADmap module. MADmap is a maximum likelihood
map making routine for converting time ordered data sets into
pixelized sky maps. Input TOD data is assumed to be calibrated and
flatfielded and input InvNtt is from calibration tree.
(Reference: http://crd.lbl.gov/~cmc/MADmap/doc/man/MADmap.html)
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 = 1e5, 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 perpixel 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:
HeaderLONG:
min_sample starting sample index
max_sample last sample index
n_correlation correlation width of matrix
DataDOUBLE:
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.
(2) Run MADmap module.
(3) If runNaiveMapper = true run NaiveMapper module.
(4) Put astrometric header information into output products.

EOF
