Pixel-based ACS/WFC CTE Forward Model

Pixel-based ACS/WFC CTE Forward Model#

This notebook demonstrates preparing data for input into the ACS/WFC pixel-based CTE forward model and running the model.

Table of Contents#

Introduction
0. Install/update stenv
1. Imports
2. Download data and reference files
3. Create an image of artificial stars

Option A: Start with an observed FLC image

4. Add artificial stars
5. Reverse the flat and dark correction
6. Run CTE forward model
7. (Optional) Run CTE correction
8. Apply flat and dark correction

Option B: Start with a synthetic image

4. Create your image
5. Run CTE forward model
6. (Optional) Run CTE correction
7. Apply flat and dark correction

About this Notebook

This notebook currently fails to execute, use as reference only

Introduction#

The charge transfer efficiency (CTE) of the Advanced Camera for Surveys (ACS) Wide Field Channel (WFC) has been decreasing over the lifetime of the instrument. Radiation damage from cosmic rays and other sources leads to charge traps within the detector. These traps remove electrons from charge packets as they are transferred between rows of the detector, and release the electrons in subsequent pixels. This causes flux to be removed from bright features and released into pixels behind the features (relative to the row closest to the amplifier), creating bright trails.

A pixel-based CTE correction model for the ACS/WFC detector is fully described in Anderson & Bedin (2010), and a recent update to the model is presented in ACS ISR 2018-04. The model is based on an empirical determination of the number and depth of charge traps distributed across the detector. It simulates detector readout of an input image, removes the result from the input, and iterates five times. In this way, a reverse model is successively approximated by the forward model. Electrons released in trails are removed and added back to the bright feature in which they originated.

The pixel-based correction was implemented in the calibration pipeline code for ACS (CALACS) in 2012 and the algorithm was updated and improved in 2018. The CTE correction step within CALACS runs on bias-corrected images, blv_tmp files, producing blc_tmp files, which lack the bright trails due to poor CTE. Further calibration, including dark correction and flat-fielding, produces flt and flc files from the blv_tmp and blc_tmps files, respectively. For more information on calibration of ACS/WFC data, see the ACS Data Handbook.

Users desiring to more fully understand the effects of pixel-based CTE correction on their science may wish to run the forward model (i.e., the detector readout simulation) on data containing artificial stars. Here we demonstrate two methods for running the CTE forward model. In Option A, we begin with an observed flc image, whereas in Option B, we begin with a raw image and generate synthetic data on which to run the forward model.

Note: The forward model, like the CTE correction step in CALACS, adds 10% of the difference between the input and output SCI extensions to the ERR extensions to account for uncertainty in the CTE model. Below, we provide guidance for properly repopulating the ERR extensions of forward-modeled data.

0. Install/update stenv#

*AstroConda is no longer supported and is superseded by stenv.#

stenv will include most packages from AstroConda and is recommended to process and analyze data from the Hubble Space Telescope (HST) and James Webb Space Telescope (JWST). To install and activate stenv, please refer to the documentation. NOTE: stenv requires Python 3.8 or greater.

If you already use stenv, make sure the versions of hstcal and acstools are both at least 2.1.0. The version of CALACS should be at least 10.1.0. Check the versions of all three on the command line:

$ conda list hstcal

$ conda list acstools

$ calacs.e --version

It is recommended, however, that you use the most up-to-date versions of these packages. To update these packages, run the following via the command line:

$ conda update hstcal

$ conda update acstools

1. Imports#

Start by importing several packages:

matplotlib notebook for creating interactive plots
os for setting environment variables
shutil for managing directories
numpy for math and array calculations
collections OrderedDict for making dictionaries easily
matplotlib pyplot for plotting
matplotlib.colors LogNorm for scaling images
astropy.io fits for working with FITS files
photutils datasets for creating synthetic stars and images
astroquery.mast Observations for downloading data from MAST
acstools acsccd for performing bias correction
acstools acscte for performing CTE correction (reversing CTE trailing)
acstools acs2d for performing dark correction and flat-fielding
acstools acscteforwardmodel for running CTE forward model (generating CTE trailing)

import os
import shutil
import numpy as np
from collections import OrderedDict
import matplotlib.pyplot as plt
from astropy.io import fits
from photutils import datasets
from astroquery.mast import Observations
from acstools import acsccd
from acstools import acscte
from acstools import acs2d
from acstools import acscteforwardmodel

Top of Page

2. Download data and reference files#

Full-frame, new-mode subarray, and 2K old-mode subarray ACS/WFC images can be run through the CTE forward model. New-mode subarrays were added to the HST flight software at the beginning of Cycle 24. These subarrays have APERTURE keywords of the type WFC1A-512, WFC1A-1K, WFC1A-2K, etc. Old-mode subarrays have APERTURE keywords of the type WFC1-512, WFC1-1K, WFC1-2K, etc. WFC apertures are also listed in Table 7.7 of the ACS IHB.

We recommend that the CTE forward model be run on data that has been bias-corrected, but not dark-corrected or flat-fielded. The flat and dark should be present in the image input into the CTE forward model because these features are present in the image when it is read out, and are therefore affected by CTE losses. The forward model can be run on flc files, but the results will technically be incorrect. Photometric tests of forward modeled data of both types show minor differences. Post-SM4 subarray data must be destriped with acs_destripe_plus, which will also perform the other calibration steps. Note: At this time, acs_destripe_plus only produces flt/flc images.

We download a full-frame 47 Tuc image, jd0q14ctq, from the ACS CCD Stability Monitor program (PI: Coe, 14402) from the Mikulski Archive for Space Telescopes (MAST) using astroquery. This image was taken in March 2016, and so it is strongly affected by CTE losses. We download the flc image for Option A and the raw image for Option B into the current directory.

obs_table = Observations.query_criteria(obs_id='jd0q14ctq')

dl_table = Observations.download_products(obs_table['obsid'], mrp_only=False,
                                          productSubGroupDescription=['FLC', 'RAW'])

for row in dl_table:
    oldfname = row['Local Path']
    unique_fname = np.unique(oldfname)
    newfname = os.path.basename(oldfname)
    print(row)
    os.rename(oldfname, newfname)
    
shutil.rmtree('mastDownload')

Next, update and download the correct flat and dark reference files for the jd0q14ctq dataset from the Calibration Reference Data System (CRDS). We use the CRDS command line tools to do this.

os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'
os.environ['CRDS_PATH'] = './crds_cache'
os.environ['jref'] = './crds_cache/references/hst/acs/'

!crds bestrefs --update-bestrefs --sync-references=1 --files jd0q14ctq_flc.fits

Next, we obtain the filenames of the flat (PFLTFILE), CTE-corrected dark (DRKCFILE), and superbias (BIASFILE) reference files from the image header. The flat will be used to add the effects of the flat field back into the image. If the data were post-flashed, then the flash file (FLSHFILE) is needed as well. This is shown in the commented line below. The CTE-corrected dark will be used to add dark current to the flc file and the synthetic data. The superbias will be used to repopulate the ERR extensions of the forward-modeled image.

hdr = fits.getheader('jd0q14ctq_flc.fits')

flat = hdr['PFLTFILE'].split('$')[-1]
dkc = hdr['DRKCFILE'].split('$')[-1]
bias = hdr['BIASFILE'].split('$')[-1]
# flash = hdr['FLSHFILE'].split('$')[-1]

We open the flat and dark images and obtain the SCI extensions of both CCDs, which are extension 1 for WFC2 and extension 4 for WFC1. We open the superbias image and obtain the ERR extensions of both CCDs, which are extension 2 for WFC2 and extension 5 for WFC2. If the flash file is needed, obtain the SCI extensions for both CCDs. This is shown in the commented out lines below.

# The jref environment variable gives the directory containing the reference files
flat_hdu = fits.open('{}/{}'.format(os.environ['jref'], flat))

flat_wfc1 = flat_hdu[4].data
flat_wfc2 = flat_hdu[1].data

dkc_hdu = fits.open('{}/{}'.format(os.environ['jref'], dkc))

dkc_wfc1 = dkc_hdu[4].data
dkc_wfc2 = dkc_hdu[1].data

# Darks can sometimes have negative pixels because they have been flash-corrected 
# and CTE-corrected. Set all negative pixels to zero
dkc_wfc1[dkc_wfc1 < 0.] = 0.
dkc_wfc2[dkc_wfc2 < 0.] = 0.

bias_hdu = fits.open('{}/{}'.format(os.environ['jref'], bias))

err_bias_wfc1 = bias_hdu[5].data
err_bias_wfc2 = bias_hdu[2].data

dq_bias_wfc1 = bias_hdu[6].data
dq_bias_wfc2 = bias_hdu[3].data

# flash_dhu = fits.open('{}/{}'.format(os.environ['jref'], flash))

# flash_wfc1 = flash_hdu[4].data
# flash_wfc2 = flash_hdu[1].data

Top of Page

3. Create image of artificial stars#

We will need an image containing artificial stars on zero background for both options presented below. Artificial stars are typically generated using models that do not include the flat field, or are produced from data that have been flat-fielded. If this is the case, artificial stars should be added to the image at the flc stage.

Users of this notebook may have a preferred method for generating artificial stars and adding them to data, so here we simply add several Gaussians to the image using utilities within photutils.datasets in astropy. These Gaussians are not representative of the true ACS/WFC PSF, and are added here for illustrative purposes only. Please note that artificial sources with peak values approaching or exceeding the WFC CCD full well value of about 80,000 electrons are not recommended for simulated data. Blooming of charge from saturated pixels is not implemented in this example.

There are many tools for generating artificial stars, including Tiny Tim, effective PSFs, or EPSFBuilder in photutils. A recent study of PSF models for ACS/WFC can be found here.

First, we generate a table of random Gaussian sources of typical brightness for our 47 Tuc field with $\mathrm{FWHM}\sim2.5$ pixels. Because $\mathrm{FWHM} = 2.355\sigma$, we will generate Gaussian sources with $\sigma \sim 1.06$ pixels in both $x$ and $y$.

n_sources = 300
param_ranges = [('amplitude', [500, 30000]),
                ('x_mean', [0, 4095]),
                ('y_mean', [0, 2047]),
                ('x_stddev', [1.05, 1.07]),
                ('y_stddev', [1.05, 1.07]),
                ('theta', [0, np.pi])]

param_ranges = OrderedDict(param_ranges)

sources = datasets.make_random_gaussians_table(n_sources, param_ranges, 
                                               seed=12345)

print(sources)

Next, we get the shape of one of the flc image SCI extensions and make an image from the table of Gaussian sources. Note that this step may take a few minutes. Finally, we run the synthetic image through a Poisson sampler in order to simulate the Poisson noise of the scene.

wfc2 = fits.getdata('jd0q14ctq_flc.fits', ext=1)

shape = wfc2.shape

synth_stars_image = datasets.make_gaussian_sources_image(shape, sources)

synth_stars_image = np.random.poisson(synth_stars_image)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(synth_stars_image, vmin=0, vmax=200, interpolation='nearest', 
          cmap='Greys_r', origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=10, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

Option A: Start with an observed FLC image#

In this Option, we add synthetic stars to the scene in an flc image and process it appropriately for use with the forward model. We use the flc image because it is the closest approximation to a pristine image of the sky. Below we plot a portion of the downloaded 47 Tuc image. Stars of various magnitudes are visible, as well as cosmic rays.

flc = fits.getdata('jd0q14ctq_flc.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(flc, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', origin='lower')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

Top of Page

4A. Add artificial stars#

Add the image of artificial stars generated above to both CCDs of the flc image, and save it as a new file.

hdu = fits.open('jd0q14ctq_flc.fits')

wfc1 = hdu[4].data
wfc2 = hdu[1].data

wfc1 += synth_stars_image
wfc2 += synth_stars_image

hdu.writeto('jd0q14ctq_stars_flc.fits', overwrite=True)

Below we plot a section of the flc image with a few artificial stars circled in red.

flc_stars = fits.getdata('jd0q14ctq_stars_flc.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(flc_stars, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', 
          origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=20, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

Top of Page

5A. Reverse the flat and dark correction#

First, calculate the total exposure time of the 47 Tuc image by combining the exposure time, flash time (if any), and 3 seconds of extra dark time to approximate instrument commanding overheads.

hdu = fits.open('jd0q14ctq_stars_flc.fits')

hdr = hdu[0].header

exptime = hdr['EXPTIME']
flashtime = hdr['FLASHDUR']
darktime = exptime + flashtime + 3

Next, we open the 47 Tuc image and obtain the SCI extensions of both CCDs. We multiply by the flat and scale the CTE-corrected dark by the total exposure time. We also run the scaled dark image through a Poisson sampler to include Poisson noise in the dark scene. We then add the dark current to the image. We save the result, which is now effectively a blc_tmp file. If the data are post-flashed, we also need to reverse the flash correction. We do this by multiplying the flash file by the flash duration, running it through a Poisson sampler, and adding it to the 47 Tuc image. The lines for this are commented out below. Note: It is not recommended to use a simulated exposure time that scales pixels in the dark or flash image to or above the full well depth of ~80,000 electrons.

wfc1 = hdu[4].data
wfc2 = hdu[1].data

wfc1 *= flat_wfc1
wfc2 *= flat_wfc2

wfc1 += np.random.poisson(dkc_wfc1*darktime)
wfc2 += np.random.poisson(dkc_wfc2*darktime)

# wfc1 += np.random.poisson(flash_wfc1*flashtime)
# wfc2 += np.random.poisson(flash_wfc2*flashtime)

hdu.writeto('jd0q14ctq_stars_pfl_dkc.fits', overwrite=True)

Finally, we update the header keyword PCTECORR to PERFORM, which is necessary for running the forward model.

fits.setval('jd0q14ctq_stars_pfl_dkc.fits', 'PCTECORR', value='PERFORM')

Top of Page

6A. Run CTE forward model#

We are now ready to run the CTE forward model, which simulates the effects of CTE losses while reading out the detector. In this example, we will use the acstools module acscteforwardmodel. Note that this step may take a few minutes. The resulting filename will be *_ctefmod.fits. We rename the file to have the suffix *_blv_tmp.fits so that it can be processed by CALACS in a later step.

acscteforwardmodel.acscteforwardmodel('jd0q14ctq_stars_pfl_dkc.fits')

os.rename('jd0q14ctq_stars_pfl_dkc_ctefmod.fits', 'jd0q14ctq_stars_ctefmod_blv_tmp.fits')

After the forward model is run, the SCI extensions of the image are equivalent to a blv_tmp file, in principle. However, the ERR extensions are the original flc ERR extensions plus 10% of the forward model correction. To ensure the ERR extensions are accurate for a blv_tmp file, we will set every pixel to zero and calculate new values for each pixel according to

$\mathrm{ERR} = \sqrt{\mathrm{SCI} + \mathrm{RN}^2 + (\mathrm{ERR}_{\mathrm{superbias}}g)^2}$,

where $\mathrm{SCI}$ is the pixel value in the SCI extension (all negative pixels are set to zero), $\mathrm{RN}$ is the readnoise, $\mathrm{ERR}_{\mathrm{superbias}}$ is the pixel value in the ERR extension of the superbias, and $g$ is the gain.

First, we access the header and SCI and ERR extensions of the forward-modeled data.

hdu = fits.open('jd0q14ctq_stars_ctefmod_blv_tmp.fits')

sci_wfc1 = hdu[4].data
sci_wfc2 = hdu[1].data

err_wfc1 = hdu[5].data
err_wfc2 = hdu[2].data

hdr = hdu[0].header

Next, we take the readnoise and gain values for each quadrant from the header.

rn_A = hdr['READNSEA']
rn_B = hdr['READNSEB']
rn_C = hdr['READNSEC']
rn_D = hdr['READNSED']

gain_A = hdr['ATODGNA']
gain_B = hdr['ATODGNB']
gain_C = hdr['ATODGNC']
gain_D = hdr['ATODGND']

Finally, we make copies of the SCI extensions in which to set all negative values to zero. We calculate the appropriate error for each quadrant and save them to the ERR extensions of the forward-modeled image.

sci_wfc1_pos = np.copy(sci_wfc1)
sci_wfc2_pos = np.copy(sci_wfc2)

sci_wfc1_pos[sci_wfc1_pos < 0] = 0
sci_wfc2_pos[sci_wfc2_pos < 0] = 0

err_A = np.sqrt(sci_wfc1_pos[:, :2048] + rn_A**2 + (err_bias_wfc1[20:, 24:2072]*gain_A)**2)
err_B = np.sqrt(sci_wfc1_pos[:, 2048:] + rn_B**2 + (err_bias_wfc1[20:, 2072:-24]*gain_B)**2)
err_C = np.sqrt(sci_wfc2_pos[:, :2048] + rn_C**2 + (err_bias_wfc2[:-20, 24:2072]*gain_C)**2)
err_D = np.sqrt(sci_wfc2_pos[:, 2048:] + rn_D**2 + (err_bias_wfc2[:-20, 2072:-24]*gain_D)**2)

err_wfc1[:] = np.hstack((err_A, err_B))
err_wfc2[:] = np.hstack((err_C, err_D))

hdu.writeto('jd0q14ctq_stars_ctefmod_blv_tmp.fits', overwrite=True)

Top of Page

(Optional) 7A. Run CTE correction#

If desired, we now CTE correct the forward-modeled image. To do this, we need to update the PCTECORR keyword to PERFORM again and update the NEXTEND keyword to 6, the number of extensions left after running the forward model. (This is because the forward model strips the distortion-related extensions from the input file, but does not update NEXTEND.) Finally, run acscte on the image. The resulting filename will be *_blc_tmp.fits.

fits.setval('jd0q14ctq_stars_ctefmod_blv_tmp.fits', 'PCTECORR', value='PERFORM')
fits.setval('jd0q14ctq_stars_ctefmod_blv_tmp.fits', 'NEXTEND', value=6)

acscte.acscte('jd0q14ctq_stars_ctefmod_blv_tmp.fits')

8A. Apply flat and dark correction#

Finally, we flat-field and dark-correct the forward-modeled image using acs2d to produce an flt-like image. First, we update the keywords DARKCORR and FLATCORR to PERFORM. All other relevant CALACS header keyword switches are set to COMPLETE because the original data was an flc image. The resulting filename will be *_flt.fits. Note: If the data were post-flashed, and the flash background was added back in during an earlier step, we must also set FLSHCORR equal to PERFORM. This option is shown below in a commented out line.

fits.setval('jd0q14ctq_stars_ctefmod_blv_tmp.fits', 'DARKCORR', value='PERFORM')
fits.setval('jd0q14ctq_stars_ctefmod_blv_tmp.fits', 'FLATCORR', value='PERFORM')
# fits.setval('jd0q14ctq_stars_ctefmod_blv_tmp.fits', 'FLSHCORR', value='PERFORM')

acs2d.acs2d('jd0q14ctq_stars_ctefmod_blv_tmp.fits')

If the forward-modeled image was CTE-corrected in Step 7A, we run acs2d on the CTE-corrected image. The resulting filename will be *_flc.fits.

fits.setval('jd0q14ctq_stars_ctefmod_blc_tmp.fits', 'DARKCORR', value='PERFORM')
fits.setval('jd0q14ctq_stars_ctefmod_blc_tmp.fits', 'FLATCORR', value='PERFORM')
# fits.setval('jd0q14ctq_stars_ctefmod_blv_tmp.fits', 'FLSHCORR', value='PERFORM')

acs2d.acs2d('jd0q14ctq_stars_ctefmod_blc_tmp.fits')

The 47 Tuc image(s) are now prepared for further analysis appropriate for the user’s science. The cells below plot a portion of the final images, the flt and, if produced, the flc.

flt_stars = fits.getdata('jd0q14ctq_stars_ctefmod_flt.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(flt_stars, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', 
          origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=20, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

flc_stars = fits.getdata('jd0q14ctq_stars_ctefmod_flc.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(flc_stars, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', 
          origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=20, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

Top of Page

Option B: Start with a synthetic image#

In this Option, we start with the raw file from the jd0q14ctq dataset and process the SCI extensions to make a completely synthetic dataset. We then process it appropriately for use in the forward model.

4B. Create a synthetic image#

We will create an image that is equivalent to a blc_tmp file. This means that the image is not flat-fielded and includes sky background, Poisson noise from the sky, artificial stars or other sources, Poisson noise from the sources, dark current, and dark noise.

We first process the data with acsccd within acstools to create a blv_tmp. This ensures that the error (ERR) and data quality (DQ) extensions are created and the header keywords for every extension are initially populated, which is necessary for the CTE forward model to run.

acsccd.acsccd('jd0q14ctq_raw.fits')

Next, we obtain the SCI extensions from the blv_tmp file, and set the pixels to zero.

hdu = fits.open('jd0q14ctq_blv_tmp.fits')
    
wfc1 = hdu[4].data
wfc2 = hdu[1].data

wfc1[:] = np.zeros(shape)
wfc2[:] = np.zeros(shape)

We then generate an image containing a user-selected sky background level (here we choose 40 electrons arbitrarily) and Poisson noise. Finally, we sum the noise image and the artificial star image, and save the result as a new file.

noise_image = datasets.make_noise_image(shape, distribution='poisson', mean=40, seed=12345)

wfc1 += noise_image + synth_stars_image
wfc2 += noise_image + synth_stars_image

hdu.writeto('synth.fits', overwrite=True)

Below we plot a portion of the synthetic image, with sources circled in red.

synth = fits.getdata('synth.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(synth, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', 
          origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=20, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

In order to properly remove the dark current after the forward model is run, we update the exposure time and flash duration header keywords, EXPTIME and FLASHDUR, with the desired simulated exposure time. The total exposure time used for the dark correction is the combination of exposure time, flash duration, and overhead (which depends on whether the data are post-flashed or not). Here we set EXPTIME to 300 seconds and FLASHDUR to 0 seconds. In principle one should retrieve the overheads from CCDTAB, but we provide a simpler approach here. Please note that the overhead times may change.

exptime = 300.
flashdur = 0.

hdu = fits.open('synth.fits')

hdr = hdu[0].header

hdr['EXPTIME'] = exptime
hdr['FLASHDUR'] = flashdur

if flashdur > 0:
    overhead = 2.43
else:
    overhead = 0.21

darktime = exptime + flashdur + overhead

We need to update the FITS header with the correct dark time for use in a later step. However, we leave off the overhead, because ACSCCD (run in a later step) will automatically append the overhead to the dark time in the header. We still use the darktime variable with overhead to manually add the dark current to our model below.

hdr['DARKTIME'] = darktime - overhead

Next, obtain the SCI extensions of both CCDs. We then multiply by the flat and scale the CTE-corrected dark by a chosen exposure time. We also run the scaled dark image through a Poisson sampler to include Poisson noise in the dark scene. We then add the dark current to the image. If post-flash is desired, multiply the flash reference file by the flash duration, run it through a Poisson sampler, and add to the synthetic data. This is shown in the commented out lines below. We save the result, which is now effectively a blc_tmp file. Note: It is not recommended to use a simulated exposure time that scales pixels in the dark or flash image to or above the full well depth of ~80,000 electrons.

Note that these reference files are specific to the anneal cycle in which these data were taken. If an observation date other than that listed in the DATE-OBS header keyword is desired for the synthetic data, different reference files will be needed. These can be found by updating the DATE-OBS header keyword in the synthetic image to the desired observation date, and rerunning the cell in Step 2 which uses CRDS bestrefs to download the correct reference files.

wfc1 = hdu[4].data
wfc2 = hdu[1].data

wfc1 *= flat_wfc1
wfc2 *= flat_wfc2

wfc1 += np.random.poisson(dkc_wfc1*darktime)
wfc2 += np.random.poisson(dkc_wfc2*darktime)

# wfc1 += np.random.poisson(flash_wfc1*flashdur)
# wfc2 += np.random.poisson(flash_wfc2*flashdur)

hdu.writeto('synth_blc_tmp.fits', overwrite=True)

Top of Page

5B. Run CTE forward model#

We are now ready to run the CTE forward model, which simulates the effects of CTE losses while reading out the detector. In this example, we will use the acstools module acscteforwardmodel. Note that this step may take a few minutes. The resulting filename will be *_ctefmod.fits. We rename this to *_blv_tmp.fits in order to ensure the correct behavior from acs2d in a later step.

acscteforwardmodel.acscteforwardmodel('synth_blc_tmp.fits')

os.rename('synth_blc_tmp_ctefmod.fits', 'synth_ctefmod_blv_tmp.fits')

At this point, we also add readnoise to the forward-modeled SCI extensions to complete the readout simulation. We find the readnoise values for each quadrant of the image from the header keywords READNSEA, READNSEB, etc. We make a noise image for each quadrant, concatenate quadrants A and B and quadrants C and D, and add them to the synthetic image.

rn_A = hdr['READNSEA']
rn_B = hdr['READNSEB']
rn_C = hdr['READNSEC']
rn_D = hdr['READNSED']

img_rn_A = datasets.make_noise_image((shape[0], int(shape[1]/2)), distribution='gaussian', 
                                     mean=0., stddev=rn_A)
img_rn_B = datasets.make_noise_image((shape[0], int(shape[1]/2)), distribution='gaussian', 
                                     mean=0., stddev=rn_B)
img_rn_C = datasets.make_noise_image((shape[0], int(shape[1]/2)), distribution='gaussian', 
                                     mean=0., stddev=rn_C)
img_rn_D = datasets.make_noise_image((shape[0], int(shape[1]/2)), distribution='gaussian', 
                                     mean=0., stddev=rn_D)

wfc1_rn = np.hstack((img_rn_A, img_rn_B))
wfc2_rn = np.hstack((img_rn_C, img_rn_D))

hdu = fits.open('synth_ctefmod_blv_tmp.fits')

wfc1 = hdu[4].data
wfc2 = hdu[1].data

wfc1 += wfc1_rn
wfc2 += wfc2_rn

hdu.writeto('synth_ctefmod_rn_blv_tmp.fits', overwrite=True)

The SCI extensions of this image are now equivalent to a blv_tmp file, in principle. However, the ERR extensions are the original blv_tmp ERR extensions plus 10% of the forward model correction. To ensure the ERR extensions are accurate, we will calculate new values for each pixel according to

$\mathrm{ERR} = \sqrt{\mathrm{SCI} + \mathrm{RN}^2 + (\mathrm{ERR}_{\mathrm{superbias}}g)^2}$,

where $\mathrm{SCI}$ is the pixel value in the SCI extension (all negative pixels are set to zero), $\mathrm{RN}$ is the readnoise, $\mathrm{ERR}_{\mathrm{superbias}}$ is the pixel value in the ERR extension of the superbias, and $g$ is the gain.

First, we access the header and SCI and ERR extensions of the forward-modeled data.

hdu = fits.open('synth_ctefmod_rn_blv_tmp.fits')

sci_wfc1 = hdu[4].data
sci_wfc2 = hdu[1].data

err_wfc1 = hdu[5].data
err_wfc2 = hdu[2].data

hdr = hdu[0].header

Next, we take the readnoise and gain values for each quadrant from the header.

rn_A = hdr['READNSEA']
rn_B = hdr['READNSEB']
rn_C = hdr['READNSEC']
rn_D = hdr['READNSED']

gain_A = hdr['ATODGNA']
gain_B = hdr['ATODGNB']
gain_C = hdr['ATODGNC']
gain_D = hdr['ATODGND']

Finally, we make copies of the SCI extensions in which to set all negative values to zero. We calculate the appropriate error for each quadrant and save them to the ERR extensions of the forward-modeled image.

sci_wfc1_pos = np.copy(sci_wfc1)
sci_wfc2_pos = np.copy(sci_wfc2)

sci_wfc1_pos[sci_wfc1_pos < 0] = 0
sci_wfc2_pos[sci_wfc2_pos < 0] = 0

# The superbias ERR arrays contain 20 rows of virtual overscan at the edge of 
# each CCD furthest from the amplifier and 24 columns of physical prescan on the
# left and right edges.
err_A = np.sqrt(sci_wfc1_pos[:, :2048] + rn_A**2 + (err_bias_wfc1[20:, 24:2072]*gain_A)**2)
err_B = np.sqrt(sci_wfc1_pos[:, 2048:] + rn_B**2 + (err_bias_wfc1[20:, 2072:-24]*gain_B)**2)
err_C = np.sqrt(sci_wfc2_pos[:, :2048] + rn_C**2 + (err_bias_wfc2[:-20, 24:2072]*gain_C)**2)
err_D = np.sqrt(sci_wfc2_pos[:, 2048:] + rn_D**2 + (err_bias_wfc2[:-20, 2072:-24]*gain_D)**2)

err_wfc1[:] = np.hstack((err_A, err_B))
err_wfc2[:] = np.hstack((err_C, err_D))

hdu.writeto('synth_ctefmod_rn_blv_tmp.fits', overwrite=True)

We also repopulate the DQ extensions because they reflect the processing of the original blv_tmp SCI extensions. To do this, we reprocess the data with DQICORR and SINKCORR, both within acsccd, and add in the DQ extensions of the appropriate superbias file. First, we rename the synthetic image to have a filename *_raw.fits, or acsccd will fail. Then, we update the DQICORR and SINKCORR header keywords to PERFORM and run acsccd. The output will be *_blv_tmp.fits again.

This step will also add the proper overhead to DARKTIME header keyword, but only if using CALACS v10.3.3 or later. We first check the version.

os.system('calacs.e --version')

If your version of CALCACS is prior to v10.3.3, please update. Now let’s run ACSCCD.

os.rename('synth_ctefmod_rn_blv_tmp.fits', 'synth_ctefmod_rn_raw.fits')

fits.setval('synth_ctefmod_rn_raw.fits', 'DQICORR', value='PERFORM')
fits.setval('synth_ctefmod_rn_raw.fits', 'SINKCORR', value='PERFORM')

acsccd.acsccd('synth_ctefmod_rn_raw.fits')

Next, obtain the DQ extensions of the superbias file, and add them to the new DQ extensions of the synthetic data with a bitwise_or operator.

hdu = fits.open('synth_ctefmod_rn_blv_tmp.fits')

dq_wfc1 = hdu[6].data
dq_wfc2 = hdu[3].data

dq_wfc1[:] = np.bitwise_or(dq_wfc1, dq_bias_wfc1[20:, 24:-24])
dq_wfc2[:] = np.bitwise_or(dq_wfc2, dq_bias_wfc2[:-20, 24:-24])

Also, acsccd assumes the input image is in DN and converts it to counts by multiplying by the gain. We defined our model in electrons though, so the acsccd conversion to electrons was unnecessary. We divide by the gains to fix this and replace the image.

sci_wfc1 = hdu[4].data
sci_wfc2 = hdu[1].data

sci_wfc1[:, :2048] /= gain_A
sci_wfc1[:, 2048:] /= gain_B
sci_wfc2[:, :2048] /= gain_C
sci_wfc2[:, 2048:] /= gain_D

err_wfc1 = hdu[5].data
err_wfc2 = hdu[2].data

err_wfc1[:, :2048] /= gain_A
err_wfc1[:, 2048:] /= gain_B
err_wfc2[:, :2048] /= gain_C
err_wfc2[:, 2048:] /= gain_D

hdu.writeto('synth_ctefmod_rn_blv_tmp.fits', overwrite=True)

Top of Page

6B. (Optional) Run CTE correction#

If desired, we now CTE correct the forward-modeled image. To do this, we need to update the PCTECORR keyword to PERFORM again, and run acscte on the image. The resulting filename will be *_blc_tmp.fits.

fits.setval('synth_ctefmod_rn_blv_tmp.fits', 'PCTECORR', value='PERFORM')

acscte.acscte('synth_ctefmod_rn_blv_tmp.fits')

7B. Apply flat and dark correction#

Finally, we flat-field and dark-correct the forward-modeled image using acs2d to produce an flt-like image. We first ensure that the DARKCORR, FLATCORR, and if necessary, FLSHCORR header keywords are set to PERFORM. The resulting filename will be *_flt.fits.

fits.setval('synth_ctefmod_rn_blv_tmp.fits', 'DARKCORR', value='PERFORM')
fits.setval('synth_ctefmod_rn_blv_tmp.fits', 'FLATCORR', value='PERFORM')
# fits.setval('synth_ctefmod_rn_blv_tmp.fits', 'FLSHCORR', value='PERFORM')

acs2d.acs2d('synth_ctefmod_rn_blv_tmp.fits')

If the forward-modeled image was CTE-corrected in Step 6, we run acs2d on the CTE-corrected image. The resulting filename will be *_flc.fits.

fits.setval('synth_ctefmod_rn_blc_tmp.fits', 'DARKCORR', value='PERFORM')
fits.setval('synth_ctefmod_rn_blc_tmp.fits', 'FLATCORR', value='PERFORM')
# fits.setval('synth_ctefmod_rn_blc_tmp.fits', 'FLSHCORR', value='PERFORM')

acs2d.acs2d('synth_ctefmod_rn_blc_tmp.fits')

The image(s) are now prepared for further analysis appropriate for the user’s science. The cells below plot a portion of the final images, the flt and, if produced, the flc.

synth_ctefmod = fits.getdata('synth_ctefmod_rn_flt.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(synth_ctefmod, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', 
          origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=20, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

synth_ctefmod_flc = fits.getdata('synth_ctefmod_rn_flc.fits', ext=1)

fig, ax = plt.subplots(1, 1, figsize=(9, 4))

ax.imshow(synth_ctefmod_flc, vmin=0, vmax=200, interpolation='nearest', cmap='Greys_r', 
          origin='lower')
ax.plot(sources['x_mean'], sources['y_mean'], marker='o', markersize=20, 
        markerfacecolor='none', markeredgecolor='red', linestyle='none')

ax.set_xlim(2000, 2800)
ax.set_ylim(1200, 1700)

For more help:#

More details may be found on the ACS website and in the ACS Instrument and Data Handbooks.

Please visit the HST Help Desk. Through the help desk portal, you can explore the HST Knowledge Base and request additional help from experts.

About this Notebook#

Author: Jenna Ryon, ACS Instrument Team
Updated On: 04/21/2022

Top of Page