%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@ Reducing the Lockman Hole Mosaic Images II 10.16.00 %@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@%@ The prior set of notes takes you through sky flattening. Only one additional comment is merited: Fixing the bias jump in chip 6 is actually a two step process. First, on the raw (untrimmed/unexploded) frames, run mscloverscan. This will just about fix things, but there are a few narrow residuals. Now process the data through sky flattening. (You don't have to go that far, but it's okay to do so.) Go back and look at the affected frames and note those that show the tell-tale skinny residuals. At this point, you can run fixchip6 on them (with the usual exploding of the extensions, etc), and the residuals should diminish. I found that some residuals did not extend all the way across the chip; these remain problematic, but they are few in number and extent. Now on with it! 1. Improve the wcs. a) This step is pretty easy, unless you have a catastrophic problem (which I did -- see below). The idea is to improve the astrometric solution from the raw RA and Dec in the header (which merely reflects where the telescope thinks it's pointing). Do this by comparing to a bunch of star with known positions from the USNO A2.0 catalogue. First, you set the parameters in mscgetcalalog like so: PACKAGE = mscred TASK = mscgetcatalog input = List of Mosaic files output = Output file of sources (magmin = 10.) Minimum magnitude (magmax = 21.) Maximum magnitude (catalog= NOAO:USNO-A2) Catalog (rmin = 21.) Minimum radius (arcmin) (mode = ql) This will actually fetch the catalogue for each iteration of msccmatch, which actually computes the astrometric solution. Set msccmatch like so: PACKAGE = mscred TASK = msccmatch input = obj151pfs.fits List of input mosaic exposures coords = !mscgetcat $I $C Coordinate file (ra/dec) (outcoor= ) List of updated coordinate files (usebpm = yes) Use bad pixel mask? (nsearch= 60) Maximum number of positions to use in search (search = 60.) Translation search radius (arcsec) (rsearch= 0.2) Rotation search radius (deg) (nfit = 50) Min for fit (>0) or max not found (<=0) (rms = 0.5) Maximum fit RMS to accept (arcsec) (maxshif= 3.5) Maximum centering shift (arcsec) (fitgeom= general) Fitting geometry (reject = 3.) Fitting rejection limit (sigma) (update = yes) Update coordinate systems? (interac= yes) Interactive? (fit = yes) Interactive fitting? (verbose= yes) Verbose? (listcoo= yes) List centered coordinates in verbose mode? (graphic= stdgraph) Graphics device (cursor = ) Graphics cursor accept = yes Accept solution? (mode = ql) The funny business at the "coords" parameter is the call to mscgetcat. Note that I ran this one file at a time. This turned out to be mere superstition... (It was the *second* frame that had the catastrophic error... Long story.) Go ahead and try it in batch mode. For every image, msccmatch will do a coarse fit and then a more exact fit to the USNO stars in the list generated by the call to mscgetcat. This involves an interactive window. The key buttons are: g - show you the bird's eye view of the stars being fitted f - make a new fit x - plots some residuals r - plots some residuals s - plots some residuals y - plots some residuals q - quit, and get prompted for update Basically, I cycled through x, r, s, and y and deleted outliers until my residuals sank below 0.4 arcseconds. Then I told msccmatch to go ahead and update the headers. b) If the discrepancy between where the telescope thinks it's pointed (i.e. the position written to TELRA, TELDEC, TELEQUIN) and where it's actually pointed is too large (e.g. greater then a minute or so), then msccmatch will fail with "automatic search failed". This means you need to update the wcs of the affected frames with a better estimate of where the telescope was pointed at the time that the image was taken. There are lots of ways to do this. I was lucky and had one good frame upon which msccmatch successfully ran. Hence, I knew the absolute positions of the stars in that frame and I could visually identify the same stars in the affected frames. So here's what I did: (i) Run mscexamine and hit "a" to get good Ra/Dec for one star. (ii) Run msczero on bad frames and hit "z" on the same star. Enter good Ra/Dec at prompt, then hit "q" to update. (This resets the zero point in the wcs of the affected frames with a shift in Ra and Dec.) (iii) Re-run msccmatch. It should now work. Hooray! My problem came from the fact that between obj145 and obj147 was obj146 -- which was badly trailed. The pointing of the telescope was screwed up thereafter, by about a minute of time in Ra and 2 minutes of arc in Dec. 2. Improve the Bad Pixel Masks in preparation for stacking. There are several steps that go in to getting the images ready for stacking. Improving the bad pixel masks is just one of them, and that itself involves several sub-steps. First you take care of saturated pixels (which, I think, are limited to bleed trails); then you handle cosmic rays. Once you've identified these pixels, you update the BPMs and apply the new mask to the image with ccdproc "fixpix". a) Use mkbleed4 to mask out saturated pixels. First, you need to determine the the saturation levels. Actually, there are two levels to set: the first (threshold 1) is the level above which pixels are in bleed trails; the second (threshold 2) is the level above which pixels are saturated. (The reason for this is that sometimes the bleed trails are themselves not saturated, but they still need to be removed.) Note that each threshold must be set independently for each chip! So use implot on each extension of a couple of frames to determine the levels of bleed trails and saturations. b) Now enter these levels into mkbleed4 and execute. (This step and the step that follows take a long time!) This'll make something called bleed_obj*.pl for every extension of every frame -- that's your mask of saturated pixels. Here's what I ended up using: PACKAGE = user TASK = mkbleed4 inlist = MEF image(s) for cosmic ray cleaning thrc1t1 = 39000. Value above which pixels are in bleeds for chip 1. thrc2t1 = 21000. Value above which pixels are in bleeds for chip 2. thrc3t1 = 33000. Value above which pixels are in bleeds for chip 3. thrc4t1 = 20000. Value above which pixels are in bleeds for chip 4. thrc5t1 = 12900. Value above which pixels are in bleeds for chip 5. thrc6t1 = 23000. Value above which pixels are in bleeds for chip 6. thrc7t1 = 27000. Value above which pixels are in bleeds for chip 7. thrc8t1 = 21500. Value above which pixels are in bleeds for chip 8. (thrc1t2= 39500.) Value above which pixels are saturated for chip 1. (thrc2t2= 22000.) Value above which pixels are saturated for chip 2. (thrc3t2= 36000.) Value above which pixels are saturated for chip 3. (thrc4t2= 24000.) Value above which pixels are saturated for chip 4. (thrc5t2= 27000.) Value above which pixels are saturated for chip 5. (thrc6t2= 27000.) Value above which pixels are saturated for chip 6. (thrc7t2= 30000.) Value above which pixels are saturated for chip 7. (thrc8t2= 23500.) Value above which pixels are saturated for chip 8. (radius = 1.) Grow radius (uses crgrow). (shift = 10.) Shift to be applied. (verbose= yes) Verbose. (inimgli= ) (mode = ql) c) Now use xzap to add cosmic rays to the mask. Xzap can't handle multi-extension files, so you'll have to split your frames with mscsplit. This'll make something called crmask_obj*.pl for every extension of every frame -- that's your mask of cosmic rays. PACKAGE = user TASK = xzap inlist = @list.xz_in Image(s) for cosmic ray cleaning outlist = @list.xz_out Output image(s) (zboxsz = 5) Box size for zapping (nsigma = 5.) Number of sky sigma for zapping threshold (nnegsig= 10.) Number of sky sigma for negative zapping (nrings = 1) Number of pixels to flag as buffer around CRs (nobjsig= 2.) Number of sky sigma for object identification (skyfilt= 15) Median filter size for local sky evaluation (skysubs= 1) Block averaging factor before median filtering (ngrowob= 0) Number of pixels to flag as buffer around objects (statsec= [550:1500,600:1500]) Image section to use for computing sky sigma (deletem= no) Delete CR mask after execution? (cleanpl= yes) Delete other working .pl masks after execution? (cleanim= yes) Delete working .imh images after execution? (verbose= yes) Verbose output? (checkli= yes) Check min and max pix values before filtering? (zmin = -32768.) Minimum data value for fmedian (zmax = 32767.) Minimum data value for fmedian (inimgli= ) (outimgl= ) (statlis= ) (mode = ql) with list.xz_in: obj145pfs_1 obj145pfs_2 obj145pfs_3 obj145pfs_4 obj145pfs_5 obj145pfs_6 obj145pfs_7 obj145pfs_8 obj146pfs_1 obj146pfs_2... etc. and list.xz_out: obj145pfsz_1 obj145pfsz_2 obj145pfsz_3 obj145pfsz_4 obj145pfsz_5 obj145pfsz_6 obj145pfsz_7 obj145pfsz_8 obj146pfsz_1 obj146pfsz_2... etc. d) Go ahead and rejoin the obj???obj145pfsz. e) Now you actually want to update the default BPMs with the bleed template and the CR template you just created. Do this with combinemasks: PACKAGE = mscred TASK = combinemasks inlist = @list.comb_masks Images for wich to generate new BPMs (Assumes CCDMosaThin1) (imglist= tmp$bld6335hs) (mode = ql) Note that the inlist wants the MEF filename, *not* the split extensions, with list.comb_masks: obj145pfs obj146pfs obj147pfs obj148pfs... etc. This'll make something called crmask3_obj*.pl for every extension of every frame. I think this is the default bad pixel mask updated with both saturated pixels and cosmic rays. f) Finally, use fixpix to apply the masks to the images. Fixpix does require the individual extensions, though there's no need to have them split up if you use the following syntax: PACKAGE = proto TASK = fixpix images = @list.fixpix List of images to be fixed masks = BPM List of bad pixel masks (linterp= INDEF) Mask values for line interpolation (cinterp= INDEF) Mask values for column interpolation (verbose= yes) Verbose output? (pixels = no) List pixels? (mode = ql) with list.fixpix: obj148pfsz[1] obj149pfsz[1] obj150pfsz[1] obj151pfsz[1] obj232pfsz[1] obj233pfsz[1] obj234pfsz[1] obj235pfsz[1] obj236pfsz[1] obj237pfsz[1] obj238pfsz[1] obj239pfsz[1] obj148pfsz[2] obj149pfsz[2] obj150pfsz[2] obj151pfsz[2]... etc. g) You should look at these individually when you are through! I found big problems. And the woods are DARK if you stray from the path. Here's the situation in which I found myself, and how I got out of it. I found that some bleed trails were left over, while some stars were zapped entirely. This meant that I needed to go back and drop threshold 2 in some chips and raise threshold 1 in others. However, just because I was running mkbleed4 again (thereby making a new bleed_obj*), I DID NOT want to run xzap again (thereby making the same old crmask_obj*). Xzap takes 40 minutes a pop on an ultra 2! The problem is, combinemasks looks for files with an extension z, which means it expects that you ran xzap. Rather than fool it, I just forgot about it. All combinemasks does is add the default BPM to the bleed_obj* to the crmask_obj* to make crmask3_obj*. You can do this yourself with imarith! Then you want to run fixpix again. I'm not quite sure whats going on here, but it seems that if I run fixpix without having made and extension z file, it tries to apply just the default pixel mask. However, you can be explicit and feed fixpix a list of the crmask3_obj* that you just made. Of course, you don't want to apply this to your obj???pfs file, because fixpix overwrites rather than generating a new file. So, you'll have to do a cp obj???pfs to something else before running the fixpix. I'm doing this now on obj145. We'll see.... IT WORKS. H E L L Y E A H. Ok -- here's a step-by-step account: a) Remove former bleed_obj*.pl. b) Remove former crmask3_obj*.pl. (Don't touch crmask_obj*.pl!) c) Copy obj???pfs.fits to obj???pfsz.fits. d) Run mkbleed4 on obj???pfs.fits. e) Run an imarith script like the one below. f) Run fixpix with explicit mask names: e.g. fixpix @list.fixpix @list.mask g) UPDATE HEADERS! with list.fixpix: obj148pfsz[1] obj148pfsz[2] obj148pfsz[3] obj148pfsz[4] obj148pfsz[5] obj148pfsz[6] obj148pfsz[7] obj148pfsz[8] with list.mask: crmask3_obj148pfs_1.pl crmask3_obj148pfs_2.pl crmask3_obj148pfs_3.pl crmask3_obj148pfs_4.pl crmask3_obj148pfs_5.pl crmask3_obj148pfs_6.pl crmask3_obj148pfs_7.pl crmask3_obj148pfs_8.pl And here's the imarith script (makemask.cl): imarith bleed_obj148pfs_1.pl + mscdb$noao/CCDMosaThin1/bpm_im1.pl temp_1.pl imarith bleed_obj148pfs_2.pl + mscdb$noao/CCDMosaThin1/bpm_im2.pl temp_2.pl imarith bleed_obj148pfs_3.pl + mscdb$noao/CCDMosaThin1/bpm_im3.pl temp_3.pl imarith bleed_obj148pfs_4.pl + mscdb$noao/CCDMosaThin1/bpm_im4.pl temp_4.pl imarith bleed_obj148pfs_5.pl + mscdb$noao/CCDMosaThin1/bpm_im5.pl temp_5.pl imarith bleed_obj148pfs_6.pl + mscdb$noao/CCDMosaThin1/bpm_im6.pl temp_6.pl imarith bleed_obj148pfs_7.pl + mscdb$noao/CCDMosaThin1/bpm_im7.pl temp_7.pl imarith bleed_obj148pfs_8.pl + mscdb$noao/CCDMosaThin1/bpm_im8.pl temp_8.pl imarith crmask_obj148pfs_1.pl + temp_1.pl crmask3_obj148pfs_1.pl imarith crmask_obj148pfs_2.pl + temp_2.pl crmask3_obj148pfs_2.pl imarith crmask_obj148pfs_3.pl + temp_3.pl crmask3_obj148pfs_3.pl imarith crmask_obj148pfs_4.pl + temp_4.pl crmask3_obj148pfs_4.pl imarith crmask_obj148pfs_5.pl + temp_5.pl crmask3_obj148pfs_5.pl imarith crmask_obj148pfs_6.pl + temp_6.pl crmask3_obj148pfs_6.pl imarith crmask_obj148pfs_7.pl + temp_7.pl crmask3_obj148pfs_7.pl imarith crmask_obj148pfs_8.pl + temp_8.pl crmask3_obj148pfs_8.pl imdel temp_1.pl verify=no imdel temp_2.pl verify=no imdel temp_3.pl verify=no imdel temp_4.pl verify=no imdel temp_5.pl verify=no imdel temp_6.pl verify=no imdel temp_7.pl verify=no imdel temp_8.pl verify=no fixpix @list.fixpix148 @list.masks148 Script to update headers: hedit obj145pfsz[1] BPM crmask3_obj145pfs_1.pl add- dele- verify- show+ update+ hedit obj145pfsz[2] BPM crmask3_obj145pfs_2.pl add- dele- verify- show+ update+ hedit obj145pfsz[3] BPM crmask3_obj145pfs_3.pl add- dele- verify- show+ update+ hedit obj145pfsz[4] BPM crmask3_obj145pfs_4.pl add- dele- verify- show+ update+ hedit obj145pfsz[5] BPM crmask3_obj145pfs_5.pl add- dele- verify- show+ update+ hedit obj145pfsz[6] BPM crmask3_obj145pfs_6.pl add- dele- verify- show+ update+ hedit obj145pfsz[7] BPM crmask3_obj145pfs_7.pl add- dele- verify- show+ update+ hedit obj145pfsz[8] BPM crmask3_obj145pfs_8.pl add- dele- verify- show+ update+ hedit obj147pfsz[1] BPM crmask3_obj147pfs_1.pl add- dele- verify- show+ update+ hedit obj147pfsz[2] BPM crmask3_obj147pfs_2.pl add- dele- verify- show+ update+ hedit obj147pfsz[3] BPM crmask3_obj147pfs_3.pl add- dele- verify- show+ update+ hedit obj147pfsz[4] BPM crmask3_obj147pfs_4.pl add- dele- verify- show+ update+ hedit obj147pfsz[5] BPM crmask3_obj147pfs_5.pl add- dele- verify- show+ update+ hedit obj147pfsz[6] BPM crmask3_obj147pfs_6.pl add- dele- verify- show+ update+ hedit obj147pfsz[7] BPM crmask3_obj147pfs_7.pl add- dele- verify- show+ update+ hedit obj147pfsz[8] BPM crmask3_obj147pfs_8.pl add- dele- verify- show+ update+ hedit obj148pfsz[1] BPM crmask3_obj147pfs_8.pl add- dele- verify- show+ update+ hedit obj148pfsz[2] BPM crmask3_obj147pfs_8.pl add- dele- verify- show+ update+ hedit obj148pfsz[3] BPM crmask3_obj147pfs_8.pl add- dele- verify- show+ update+ hedit obj148pfsz[4] BPM crmask3_obj147pfs_8.pl add- dele- verify- show+ update+ ...etc... 3. Generate Tangent Plane Projection a) First generate a reference image. This defines the reference point for the projection of individual to-be-stacked images to the tangent plane. Use wcsref, which is in the xwdred package (you'll have to "bye" from mscred): PACKAGE = xdwred TASK = wcsref image = Lockman_ref Reference image ra = 10.8786 RA (hours) dec = 57.4974 DEC (degrees) scale = 0.258 Scale (arcsec/pixel) rotation= 90. Rotation of DEC from N to E (degrees) (mode = ql) The RA and Dec are just the pointing for the Lockman Hole. The scale and rotation are what the DWFS guys are using. b) We're going to use mscimage to do the projecting, but mscimage has to know what to do with the gaps in the chips and with the space beyond the edges of the camera. We'll give it a mean sky value. This we feed to the parameter "blank" in mscimage below. This requires that we run iterstat on some central region of some central chip. We then run mscimage separately for each image, and we feed it the blank value that belongs to that image. (One could imagine deriving a more careful blank value -- by looking at more chips and taking a mean, for instance. Arjun claims that this doesn't help things, so you might as well not spend the time.) Note that iterstat wants an extension; it's not set up for MEF. c) Now we generate the tangent plane projections using the sky values we just got. The parameters for mscimage are below, but we'll run it with a script so as to feed it the proper values for blank and constant. (QUESTION: Arjun has blank=0 and constant=the.sky.value. The web notes have blank=the.sky.value and make not mention of constant. From the help page, it sounds like I should use the sky value for both. ANSWER: This was the right move. These pixels get masked out anyway, so blank and constant are really just cosmetic.) Here's a sample of the mscimage script. Notice the nomeclature is obj??> --> mos???. mscimage obj147pfsz mos147 reference=Lockman_ref.fits blank=10387. boundary=constant constant=10387. mscimage obj148pfsz mos148 reference=Lockman_ref.fits blank=12584. boundary=constant constant=12584. mscimage obj149pfsz mos149 reference=Lockman_ref.fits blank=11544. boundary=constant constant=11544. mscimage obj150pfsz mos150 reference=Lockman_ref.fits blank=12779. boundary=constant constant=12779. PACKAGE = mscred TASK = mscimage input = List of input mosaic exposures output = List of output images (referen= ) Reference image (pixmask= yes) Create pixel mask? (verbose= yes) Verbose output? # Resampling parmeters (blank = 0.) Blank value (interpo= sinc17) Interpolant for data (minterp= linear) Interpolant for mask (boundar= reflect) Boundary extension (constan= 0.) Constant boundary extension value (fluxcon= no) Preserve flux per unit area? (ntrim = 7) Edge trim in each extension (nxblock= 2048) X dimension of working block size in pixels (nyblock= 1024) Y dimension of working block size in pixels # Geometric mapping parameters (interac= no) Fit mapping interactively? (nx = 10) Number of x grid points (ny = 20) Number of y grid points (fitgeom= general) Fitting geometry (xxorder= 4) Order of x fit in x (xyorder= 4) Order of x fit in y (xxterms= half) X fit cross terms type (yxorder= 4) Order of y fit in x (yyorder= 4) Order of y fit in y (yxterms= half) Y fit cross terms type (fd_in = ) (fd_out = ) (fd_ext = ) (fd_coor= ) (mode = ql) 4. Remove any residual large scale gradients with mscskysub. PACKAGE = mscred TASK = mscskysub input = @list.mos Input images to be fit output = @list.mos//s Output images xorder = 2 Order of function in x yorder = 2 Order of function in y (type_ou= residual) Type of output (fit,residual,response,clean) (functio= leg) Function to be fit (legendre,chebyshev,spline3) (cross_t= yes) Include cross-terms for polynomials? (xmedian= 100) X length of median box (ymedian= 100) Y length of median box (median_= 50.) Minimum fraction of pixels in median box (lower = 0.) Lower limit for residuals (upper = 0.) Upper limit for residuals (ngrow = 0) Radius of region growing circle (niter = 0) Maximum number of rejection cycles (regions= mask) Good regions (all,rows,columns,border, (rows = *) Rows to be fit (columns= *) Columns to be fit (border = 50) Width of border to be fit (section= ) File name for sections list (circle = ) Circle specifications (mask = BPM) Mask (div_min= INDEF) Division minimum for response output (mode = ql) 5. Determine the scalings to be used for each image when combining. Here, we use mscgetcat to generate a catalog of stars. Then mscimatch goes and counts up the flux in the USNO A2 stars and figures out the image scales. You need to be a little careful, since, on the bright end, you don't want mscgetcat to feed you saturated stars, and on the faint end, you don't want to be considering patches of sky. So experiment with the magnitudes in in mscgetcat; use tvmark to see what you're actually generating. PACKAGE = mscred TASK = mscgetcatalog input = mos145s List of Mosaic files output = mos145s.usnoA2 Output file of sources (magmin = 17.) Minimum magnitude (magmax = 22.) Maximum magnitude (catalog= NOAO:USNO-A2) Catalog (rmin = 21.) Minimum radius (arcmin) (mode = ql) Now run mscimatch. It'll chug away and then bring up an interactive window. Check the fits, then quit out and see that the scales are mostly close to unity. If you accept them, the headers of the individual images will be updated with the scale values. PACKAGE = mscred TASK = mscimatch input = @list.moss List of images coords = mos145s.usnoA2 File of coordinates (bpm = BPM) List of bad pixel masks (measure= ) Measurment file (scale = yes) Determine scale? (zero = no) Determine zero offset? (box1 = 21) Inner box size for statistics (box2 = 51) Outer box size for statistics (lower = 100.) Lower limit for good data (upper = 25000.) Upper limit for good data (niterat= 3) Number of sigma clipping iterations (sigma = 3.) Sigma clipping factor (interac= yes) Interactive? (verbose= yes) Verbose? accept = yes Accept scaling and update images? (mode = ql) Incidentally, you should make a record of the FWHMs of some central star as it appears in each of the images. You may want to stack the images according to their seeing later on. 6. Your first stack! For the first composite imagine, do just a straight median (though you might leave out the bad seeing guys). PACKAGE = mscred TASK = mscstack input = @list.stack_1 List of images to combine output = Lockman_I_pass1a Output image (plfile = Lockman_I_pass1a_bpm.pl) List of output pixel list files (optional) (combine= median) Type of combine operation (median|average) (reject = none) Type of rejection (masktyp= goodvalue) Mask type (maskval= 0.) Mask value (blank = 50000.) Value if there are no pixels (scale = !mscscale) Image scaling (zero = !msczero) Image zero point offset (weight = none) Image weights (statsec= ) Image section for computing statistics (lthresh= 1.) Lower threshold (hthresh= INDEF) Upper threshold (nlow = 1) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg) (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 5.) Lower sigma clipping factor (hsigma = 5.) Upper sigma clipping factor (rdnoise= rdnoise) ccdclip: CCD readout noise (electrons) (gain = gain) ccdclip: CCD gain (electrons/DN) (snoise = 0.) ccdclip: Sensitivity noise (fraction) (sigscal= 0.1) Tolerance for sigma clipping scaling corrections (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0.) Radius (pixels) for neighbor rejection (mode = ql) 7. Now you want to set up to run diffdetect, which will compare the individual frames to the stacked frame and look for differences. Differences that are bigger than some parameter are flagged as artifacts and should included in a bad pixel mask. Before you can do this, though, you need to run imexpr. The bad pixel masks which get created by mscstack are of a different format than that which is expected by diffdetect (see the web-based notes for an explanation). In any case, run impexress on the BPB outputted by mscstack to make a revised version. Here, we took the output from mscstack "Lockman_I_stk1_bpm.pl" and turned it into "Lockman_I_diffstk1_bpm.pl". The numeral in the expr parameter is the number of frames that went into the stack. PACKAGE = imutil TASK = imexpr expr = (a != 11 ? 0 : a) expression output = Lockman_I_diffpass1a_bpm.pl output image (dims = auto) output image dimensions (intype = auto) minimum type for input operands (outtype= auto) output image pixel datatype (refim = auto) reference image for wcs etc (bwidth = 0) boundary extension width (btype = nearest) boundary extension type (bpixval= 0.) boundary pixel value (rangech= yes) perform range checking (verbose= yes) print informative messages (exprdb = none) expression database (lastout= Lockman_I_diffpass1a_bpm.pl) last output image a = Lockman_I_pass1a_bpm.pl operand a b = operand b c = operand c d = operand d e = operand e f = operand f g = operand g h = operand h i = operand i j = operand j k = operand k l = operand l m = operand m n = operand n o = operand o p = operand p q = operand q r = operand r s = operand s t = operand t u = operand u v = operand v w = operand w x = operand x y = operand y z = operand z (mode = ql) 8. Now you can run diffdetect, which is in the ace package. This will compare each frame that went into the stack to the composite image and flag pixels in the constituent images which deviate from those in the composite. The idea is to augment the BPMs which already contain info on bleed trails, CRs, and edges with pixels that differ substantially in one frame from the the same pixels in all the other frames. Some of the parameters need explaining. "Images" is the constituent image that you are comparing to the composite. "Rimages" is the reference image -- the composite, stacked image. "Objmasks" and "Catalogs" are the outputs of diffdetect -- the list of pixels in "images" which deviate from the composite in "Rimages". PACKAGE = ace TASK = diffdetect images = @list.diff_in List of images (masks = BPM) List of bad pixel masks (skys = SKYFIT) List of skys (sigmas = SKYSIG) List of sky sigmas (scales = !mscscale) List of image intensity scale factors # Reference Image(s) (rimages= Lockman_I_pass1a) List of reference images (rmasks = Lockman_I_diffpass1a_bpm.pl) List of reference bad pixel masks (rskys = SKYFIT) List of reference skys (rsigmas= SKYSIG) List of reference sky sigmas (rscales= 1) List of reference intensity scale factors # Output objmasks= diff145s List of output object masks catalogs= diff145s List of output catalogs (logfile= STDOUT) List of log files # Sky (newsky = yes) Determine new sky fit if one already exists? (nskylin= 100) Number of sky sample lines (skyblk1= 10) Sky block size for 1D sky estimation (skyhcli= 3.) High sky clipping during 1D sky estimation (skylcli= 3.) Low sky clippling during 1D sky estimation (skyxord= 4) Sky fitting x order (skyyord= 4) Sky fitting y order (skyxter= half) Sky fitting y order # Detection (convolv= block 3 3) Convolution kernel (hsigma = 2.5) Sigma threshold above sky (lsigma = 10.) Sigma threshold below sky (minpix = 10) Minimum number of pixels in detected objects (sigavg = 6.) Sigma of mean flux cutoff (sigmax = 6.) Sigma of maximum pixel (rfrac = 0.5) Minimum fraction of reference flux (bpval = 1) Output bad pixel value # Growing (ngrow = 4) Number of grow rings (mode = ql) Note that currently, Arjun runs the three steps above with a script: mscstack ("@11q2.list", "11q2Bw_1pass", plfile="11q2Bw_1pass_bpm.pl", combine="median", reject="none", masktype="goodvalue", maskvalue=0., blank=50000., scale="!mscscale", zero="!msczero", weight="none", statsec="", lthreshold=1., hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, lsigma=5., hsigma=5., rdnoise="rdnoise", gain="gain", snoise="0.", sigscale=0.1, pclip=-0.5, grow=0.) imexpr ("(a != 10 ? 0 : a)", "11q2Bw_diffpass_bpm.pl", dims="auto", intype="auto", outtype="auto", refim="auto", bwidth=0, btype="nearest", bpixval=0., rangech=yes, verbose=yes, exprdb="none", lastout="", a="11q2Bw_1pass_bpm.pl", b="") diffdetect ("@11q2.list", "@diff.out", "@diff.out", masks="BPM", skys="SKYFIT", sigmas="SKYSIG", scales="!mscscale", rimages="11q2Bw_1pass.fits", rmasks="11q2Bw_diffpass_bpm.pl", rskys="SKYFIT", rsigmas="SKYSIG", rscales="1.", logfiles="STDOUT", newsky=yes, nskylines=100, skyblk1d=10, skyhclip=3., skylclip=3., skyxorder=4, skyyorder=4, skyxterms="half", convolve="block 3 3", hsigma=2.5, lsigma=10., minpix=10, sigavg=6., sigmax=6., rfrac=0.5, bpval=1, ngrow=4) Pay special attention to how where the commands work on lists of files and where they call for single composite images or masks. Also note that this crashes at imexpr. 9. Diffdetect will return many many anomolous objects. (I'm not really sure how to diagnose what's going on yet, but I'm getting anywhere from a hundred to a couple hundred objects detected. Mike Brown says that a good number would be thirty objects.) Now you want to go back to the individual frames that went into the stack and use blkzap and satzap to remove obvious artifacts. Both these scripts operate on the object frame but modify the pixel mask named by the BPM keyword, e.g. satzap on mos145s modifies mos145_bpm.pl. Hence, it might be a good idea to copy all current mos???_bpm.pl bad pixel masks to a subdirectory (though unsatzap and unblkzap do exist to un-do mistakes). Also, it's probably a good idea to keep a record of what you mask. If you put all your zapping commands into a script, you can always run that script again if you have to back up a few steps, and you can add and comment out zap commands if you go through another iteration... A sample script is: # # mos145s # blkzap mos145s 19.69 7251.52 122.04 8219.72 # scattered light, upper left C5 blkzap mos145s 5507.54 7888.89 5634.44 8253.24 # scattered light, C6 blkzap mos145s 4965.09 8.44 5022.40 172.19 # scattered light, C3 blkzap mos145s 7.41 46.31 128.18 4109.48 # cool edge, left C1 blkzap mos145s 122.04 11.51 8340.49 91.34 # cool edge, bottom C1 - C4 blkzap mos145s 8264.76 79.06 8373.24 4111.52 # hot edge, right C4 # # mos232s # blkzap mos232s 16.68 7237.91 100.69 8221.94 # scattered light, upper left C5 blkzap mos232s 7.68 4159.82 85.69 5305.85 # cool edge, left C5 blkzap mos232s 4.68 40.70 133.69 4111.82 # cool edge, left C1 blkzap mos232s 136.69 7.70 6268.88 82.70 # cool edge, bottom C1 - C3 blkzap mos232s 6331.88 34.70 8263.94 100.70 # cool strip, bottom C4 blkzap mos232s 8257.94 43.70 8374.94 4111.82 # hot edge, right C4 blkzap mos232s 5278.85 7906.93 5347.85 8215.94 # scattered light, C6 satzap mos232s 110.44 1274.48 3440.54 42.95 hwidth=15 # sat, C1 - C2 satzap mos232s 4006.82 8251.94 4738.82 79.70 hwidth=15 # sat, C6 - C3 To execute this script, you need to be in stsdas.toolbox.imgtools. Load it up and let 'er rip. Incidentally, you could be using acedisplay and reviewproto at any time to display the mos???s files with (1) the BPM over-layed in red, and (2) the objects found by diffdetect marked with colored contours. Acedisplay shows the entire frame with all the diffdetect objects; run it on the mos???s file.. Reviewproto allows you to cycle through each object marked by diffdetect; run it on the diff???s.tab file. 10. Since diffdetect returned a lot of objects (10^2), we should re-stack the images with the updated mask, re-run imexpr on the new stacked mask, hide the old diff???s.tab and diff???s.pl files, and then re-run diffdetect, acedisplay, and reviewproto... Much better. Diffdetect is now returning 8734 objects detected 41 objects detected 49 objects detected 30 objects detected 24 objects detected 15 objects detected 131 objects detected 42 objects detected 42 objects detected 79 objects detected 29 objects detected What to do about that first guy? If the sky counts in that frame are a factor of 2 or 3 higher than in other frames: reject it. Same if the fluxes in a few faint objects are significantly lower than in other frames. (If either were the case, then you're just adding noise -- not signal -- to the composite image.) Otherwise, you could just crank up the detection sigma in diffdetect and run it again on that frame. 11. The diff???.pl files outputted by diffdetect are the final bad pixel masks. They contain the default mask, the masks generated by mkbleed4 and xzap, and the rejected objects marked by diffdetect. Now you want to make a new stack, but first you must change the mask named by the BPM keyword in the mos???s.fits files to reflect the diffdetect output. Do this with a simple script (notice that we don't have to deal with extensions any more): hedit mos232s BPM diff232.pl add- delete- verify- show+ update+ hedit mos233s BPM diff233.pl add- delete- verify- show+ update+ hedit mos238s BPM diff238.pl add- delete- verify- show+ update+ hedit mos239s BPM diff239.pl add- delete- verify- show+ update+ hedit mos240s BPM diff240.pl add- delete- verify- show+ update+ hedit mos633s BPM diff633.pl add- delete- verify- show+ update+ hedit mos634s BPM diff634.pl add- delete- verify- show+ update+ hedit mos635s BPM diff635.pl add- delete- verify- show+ update+ hedit mos636s BPM diff636.pl add- delete- verify- show+ update+ hedit mos637s BPM diff637.pl add- delete- verify- show+ update+ 12. Now make your final stack. Do a straight average with no rejection; ideally, the the rejection is handled by the updated bad pixel masks. PACKAGE = mscred TASK = mscstack input = @list.finalstack List of images to combine output = Lockman_I_H Output image (plfile = Lockman_I_H_bpm) List of output pixel list files (optional) (combine= average) Type of combine operation (median|average) (reject = none) Type of rejection (masktyp= goodvalue) Mask type (maskval= 0.) Mask value (blank = 50000.) Value if there are no pixels (scale = !mscscale) Image scaling (zero = !msczero) Image zero point offset (weight = none) Image weights (statsec= ) Image section for computing statistics (lthresh= 1.) Lower threshold (hthresh= INDEF) Upper threshold (nlow = 1) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg) (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 5.) Lower sigma clipping factor (hsigma = 5.) Upper sigma clipping factor (rdnoise= rdnoise) ccdclip: CCD readout noise (electrons) (gain = gain) ccdclip: CCD gain (electrons/DN) (snoise = 0.) ccdclip: Sensitivity noise (fraction) (sigscal= 0.1) Tolerance for sigma clipping scaling corrections (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0.) Radius (pixels) for neighbor rejection (mode = ql)