======
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


======================================
PACS MADmap Users Manual Documentation
--------------------------------------

=================
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 
(5) runMadMap


==============
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
-------------------------

See MADmap User's Guide.

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

See MADmap User's Guide.


============
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
         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, 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
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 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
    RADESYS   ICRS         
    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.

(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 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
   masks (BADPIXEL, SATURATION, GLITCH, UNCLEANCHOP).  For (BLIND)
   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.
   (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
flat-fielded 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 = 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:

Header-LONG:
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.


(2) Run MADmap module.

(3) If runNaiveMapper = true run NaiveMapper module.

(4) Put astrometric header information into output products.


---
EOF