#-----------------------------------------------------------------------------
# IRACproc: IRAC Post-BCD Processing v4.3 PSF MASKING
#-----------------------------------------------------------------------------
(c) 2005 Mike Schuster, Massimo Marengo (mmarengo@cfa.harvard.edu)
(c) 2005 SAO - Smithsonian Astrophysical Observatory

Harvard-Smithsonian Center for Astrophysics
Smithsonian Astrophysical Observatory

Mike Schuster (software integration, PDL software)

Massimo Marengo (PDL software, HDR PSFs)
email: mmarengo@cfa.harvard.edu

16 JAN 2007

The IRAC PSFs are affected by electronic artifacts which are not completely 
characterized, and are not proportional to the source brightness. Among them: 
banding, muxbleed, and column pull-up/down which affects the rows and columns 
with bright sources. Given their non-linearity, these artifacts will not be 
completely removed by a PSF subtraction. In addition, there appears to be a 
charge diffusion effect on channel 1, responsible for a broadening of source 
cores that is related to their fluence (flux times integration time, or 
photons per area; brighter sources are more broadened).

We have developed this software to mask these PSF residuals. In addition to 
having the light from a bright point-source removed, a masked image will be 
cleared of many artifacts that could lead to false source identification by 
source extraction software (such as MOPEX/APEX and DAOPHOT). This software 
may also be useful when preparing images for publication.

#-----------------------------------------------------------------------------
# PSF Masking
#-----------------------------------------------------------------------------
This step is for masking the residual artifacts from a PSF from an IRAC 
image. This software will work with any FITS image, whether PSF-subtracted or 
not, or IRAC or not.

Change dir to the target for which you wish to mask PSF residuals. If you 
have used IRACproc to subtract PSF(s) from your image(s), then go to:

user_name][ cd proc_dir/analysis/

Start PDL.

user_name][ perldl

perldl> psfmask_bcd('I<CH>_name_psfsub.fits','I<CH>_name_psfmask.fits',
{CROSS_MASK_WIDTH=>[<XWIDTH>,<YWIDTH>],CORE_MASK_RADIUS=><RADIUS>})

Type psfmask_bcd() at the PDL prompt to see available options.

The input and output images can be specified with both relative and absolute 
paths. The output image preserves the WCS grid of the input image. The output 
FITS file 'I<CH>_name_psfmask.fits' will be written to analysis/.

This is an iterative process where the optimal parameters are reached by 
visually inspecting the masked image. It may help to display the masked image 
with a linear (or square root) stretch, minimum=median_sky-0.5, and 
maximum=median_sky+0.5.

Repeat all of the above steps for channels 1-4.

Quit PDL - you will need to restart for the next target.

perldl> quit

#-----------------------------------------------------------------------------
# Available Options
#-----------------------------------------------------------------------------
 COORDS           image coordinates of PSF to be masked - [<X>,<Y>]
                   (first pix is no. 1, can be outside image boundary v0.5+)
                   default is input image center
 CENTER           PDL coordinates of PSF to be masked - [<X>,<Y>]
                   (first pix is no. 0, can be outside image boundary v0.5+)
                   default is input image center
 CROP_SIZE        size of cropped output image - centered on PSF coordinates 
                   - [<XSIZE>,<YSIZE>]
 CROP_SEC         section to copy from input image - [x1,x2,y1,y2]
                   (in image coordinates, i.e. first pix is pix no. 1)
                   this section does not need to be centered on the PSF
 CROSS_CROP_WIDTH X and Y widths to crop for cross - [<XWIDTH>,<YWIDTH>]
 CORE_MASK_RADIUS radius to mask for core
 CROSS_MASK_WIDTH X and Y widths to mask for cross - [<XWIDTH>,<YWIDTH>]
 MASK             value for MASK (default is NaN)
 BACKGRND         background pedestal level to subtract from image
 BKG_METHOD       method to use for background pedestal calculation
                  'mode' or 'median', default is 'median'
 BKG_RES          resolution parameter that increases (for BKG_RES>1), 
                  and decreases (for BKG_RES<1), the number of bins 
                  to use for mode calculation, default = 1
 BKG_CLIP         set this to use clipped histogram in mode calculation
                  default = 0
 BITPIX           The bit depth of the output pixels
 NAN              value for NaN

 If BACKGRND is set to "?" then the background pedestal level is calculated 
 from the image.

 The option "CROSS_CROP_WIDTH" will literally take a cross, with the 
 specified widths and centered on the PSF coordinates, out of the image. The 
 remaining masked core (if any) will be replaced with the median value of the 
 unmasked pixels in a square that circumscribes the core.

 When cropping the PSF cross the values for CRPIXi will be the same as though 
 the PSF cross were not cropped out. This option is meant to be used with 
 psfuncrop_bcd (v0.3+) which replaces the cropped cross with a mask. The 
 purpose of this option is so that FFT's can be run on an image with as few 
 sharp edges as possible.

#-----------------------------------------------------------------------------
# Masking Multiple PSFs in an Image
#-----------------------------------------------------------------------------
Use multimask_bcd along with the psfsub_I<CH>.par parameter file from PSF 
subtraction. This parameter file is used to set the input and output image 
filenames, and PSF coordinates. The routine will mask each PSF listed in 
psfsub_I<CH>.par. The wrapper script multimask_bcd can automatically 
determine the most appropriate PSF coordinates for masking.

Example psfsub_I<CH>.par file:
------------------------------
/<path>/file.fits  Xshift Yshift sky amplitude
proc_dir/analysis/I<CH>_name_mosaic.fits  0.0 0.0 0 1.0
/<path>/IRACproc/psf/I<CH>_psf_xsq2_best.fits  0.0 0.0 0 1.0

----------------

For each mask to have the same parameters, use the following syntax:

perldl> multimask_bcd('psfsub_I<CH>.par',
{CROSS_MASK_WIDTH=>[<XWIDTH>,<YWIDTH>],CORE_MASK_RADIUS=><RADIUS>})

Type multimask_bcd() at the PDL prompt to see available options.

----------------

For each mask to have unique parameters, use the following syntax:

perldl> multimask_bcd('psfsub_I<CH>.par',{MASK_PAR=>'multimask_I<CH>.par'})

To specify separate masking parameters for each PSF, create a PSF mask 
parameter file that can be specified with the option "MASK_PAR" (e.g. 
multimask_I<CH>.par):

user_name][ pico multimask_I<CH>.par

Example multimask_I<CH>.par file:
---------------------------------
# PSF_X PSF_Y RADIUS XWIDTH YWIDTH
 ? ? 10 5 6
 ? ? 7.2 3 2

Note: enter 0.0 for zero, to avoid problems distinguishing between 0 and null.

----------------

The wrapper multimask_bcd has the capability to set the positions of the PSFs 
in the MASK_PAR file directly. If PSF_X or PSF_Y are set to '?', then the 
positions are calculated automatically from the PSF subtraction parameter 
file psfsub_I<CH>.par.

The parameter files can be specified with either a relative or an absolute 
path. In addition, the input and output file names may be specified with the 
options "IN_IMAGE" and "OUT_IMAGE" (e.g. if different from *psfsub.fits and 
*psfmask.fits, respectively). The "BACKGRND" options are also available.

The center for the HDR PSF can be set in the wrapper multimask_bcd with the 
options "HDR_PSF_CENTER" (PDL coordinates) or "HDR_PSF_COORDS" (image 
coordinates) - the default is located at image coordinates [183,183] for the 
sqrt(2) sampled HDR PSFs.

#-----------------------------------------------------------------------------
# Software Status
#-----------------------------------------------------------------------------
Rotation of the mask is not available at this time. It is assumed that the 
PSF artifacts are registered with the rows and columns of the image.

#-----------------------------------------------------------------------------
# Bugs
#-----------------------------------------------------------------------------
During the iterative masking process DS9 may die.