#!/usr/bin/env python """ Apply a correction to ACIS fid light centroid values to compensate for variations in the ACIS detector housing temperature. Inputs --------- :: --acisengfiles=ACISENGFILES ACIS L0 engineering files with 1CBAT and 1CBBT --acenfile=ACENFILE Aspect L1 ACEN file --fidpropsfile=FIDPROPSFILE Aspect L1 FIDPROPS file --alg=ALG Centroiding algorithm --nsmooth=NSMOOTH Median filter smoothing half-width --tsample=TSAMPLE Sampling time for temperature correction (sec) Outputs --------- :: --acenfile=ACENFILE Aspect L1 ACEN file (ang_y_sm and ang_z_sm columns are updated in place) Source ------ `Prototype code ` Description ------------ In order to maintain control margin to keep the ACIS focal plane (FP) at the nominal setpoint of -119.7 C, the ACIS detector housing (DH) heater will be disabled. The ACIS fid lights are mounted on the ACIS DH. Allowing the DH temperature to vary has the undesired side effect of producing noticable shifts (up to 0.25 arcsec on time scales of hours) in the ACIS fid light positions due to thermal contraction. For the first 8 years of the mission the DH was controlled at a temperature of -60 C. During the nine 2007 Chandra Deep Field South observations the DH heater was disabled at the start and then re-enabled 2 hours before the end of each observation. These tests were performed using a variety of fid light combinations and the the shift in fid light positions were fully analyzed. All the data were well-described by a simple physical model of the fid light positions contracting or expanding radially about a single center point with a CTE that corresponds roughly to that of aluminum. The numerical fit values coming out of that analysis are captured in the 'Fid and DH calibration constants' listed below. The center of expansion was determined in ACA angular coordinates for ACIS-I observations with the nominal SIM-Z position. The most rigorous (and physically accurate) way to represent these values would be to translate all the derived angular values into physical coordinates (mm) on the ACIS focal plane. Then this correction tool would use the full transformation chain from ACIS focal plane to ACA coordinates to calculate the needed corrections. However, given the small size of the correction and fairly large uncertainties in the derived calibration parameters, we simplify processing considerably by assuming linearity throughout the system and just using offsets with respect to the fid light angular positions used in the CDF-S calibration test observations. These are captured in the fid_[yz]_ang_nom arrays. In order to compensate for the temperature-dependent shift in aspect pipeline processing an estimate of the DH temperature is needed. This is available in ACIS engineering telemetry as MSIDs 1CBAT and 1CBBT. Unfortunately the digitization of these values is approximately 2.5 C and using them directly would produce undesirable discontinuities in the aspect solution SIM DY, DZ and DTHETA values. To avoid this we take advantage of the physical characteristic of LEDs that they get brighter as the temperature is lowered. The test data were used to calibrate a linear relationship between changes in DH temperature and fid light brightness (counts). Although the fid light brightness are fairly stable over the mission, we cannot depend on stability to the level required for an absolute translation from fid brightness to DH temperature. Instead a hybrid approach is taken -- read in both the ACIS engineering temperature data and the fid data, filter to be using only good quality values, and then match with a simple average over the aspect time interval. At the same time apply the empirically determined scaling factor to convert delta counts into delta degC. This results in a high-resolution estimate of the actual DH temperature during the interval. Using the DH temperature and the derive expansion coefficients, the expected radial expansion or contraction is then removed from the observed (smoothed) fid light centroids that are used in constructing the SIM displacement. Correcting the centroids in place has the advantage that downstream tools (in particular V&V) need no modification. The original (unmodified) centroid values are still available in the unsmoothed angular centroids as well as the original CCD pixel centroids. The plot below shows an example of the original (blue) and corrected (red) centroids for obsid 8595. The key feature to note is a smooth correction during the abrupt DH warmup near the end when the heater was re-enabled during the test. This indicates that the model is doing a good job for a fast temperature shift of about 8 C. The very bottom sub-plot shows the raw ACIS DH temperature telemetry (blue) and the hybrid estimate using fid light counts. .. image:: obs8595.jpg In other cases the model does not do quite as well and one see residuals of up to about 0.1 arcsec. This is expected to be acceptable in operations because the actual temperature deltas should be less than a few degrees. Nevertheless some additional effort will go into improving the calibration values. .. image:: obs8592.jpg The constants defined in the prototype code could logically be put in a new binary table HDU in the CALDB PCAD CALALIGN product with an extension name FIDACISDH (or FID_ACIS_DH?). References ------------------ - `Initial fid cooling investigations `_ - `CDF-S fid cooling test analysis `_ Analysis and the source for this document are in the directory: /proj/sot/ska/analysis/fid_cooling """ __docformat__ = "restructuredtext" import pyfits import numpy as np import sys from glob import glob # Fid and DH Calibration constants dh_temp_base = -60.0 # Baseline ACIS DH temperature degC degC_per_count = -7.0 / 26000. # Degrees per fid light count (integrated e-) fid_CTE = 16.3e-6 # Fid light mount coeff. of thermal expansion fid_y_center_nom = -65 / 3600. # Fid light Y-center of expansion (deg) fid_z_center_nom = 1395 / 3600. # Fid light Z-center of expansion (deg) # Fid light nominal Y and Z positions relative to expansion center (deg) # These were the fid positions for the CDF-S calibration observations used # to derive the fid_[yz]_center_nom values. fid_y_ang_nom = np.array([ 928, -769, 46, 2145, -1820, 391]) / 3600. fid_z_ang_nom = np.array([-837, 845, -970, 1061, 1060, 1703]) / 3600. def main(): opt, args = get_options() # Read the input data files. acis0eng = get_acis0eng(opt.acisengfiles).data fidprops = pyfits.open(glob(opt.fidpropsfile)[0])[1].data acen_hdus = pyfits.open(glob(opt.acenfile)[0]) acen = acen_hdus[1].data # Make sure there is at least one good fid to correct if (fidprops.field('id_status') == 'GOOD').shape[0] == 0: print 'No good fids - skipping corr_fid_cent' sys.exit(0) # Calculate detector housing temperature using ACIS eng. telemetry and fid counts times, dh_temp = calc_dh_temp(acis0eng, fidprops, acen, opt) # Apply centroid correction to acen values in-place apply_cent_corr(times, dh_temp, acen, fidprops, opt) # For development, optionally make a plot showing results of correction if opt.plotfile: make_plot(opt.plotfile, times, dh_temp, acen, fidprops, acis0eng) # Write updated values back to original file, OR do whatever # makes sense within pipeline processing. def get_options(): """Get program options.""" from optparse import OptionParser global opt, args parser = OptionParser() parser.set_defaults() parser.add_option("--acisengfiles", default='acisf*_2_eng0.fits*', help="ACIS L0 engineering files with 1CBAT and 1CBBT", ) parser.add_option("--acenfile", default='pcadf*_acen1.fits*', help="Aspect L1 ACEN file", ) parser.add_option("--fidpropsfile", default='pcadf*_fidpr1.fits*', help="Aspect L1 FIDPROPS file", ) parser.add_option("--alg", default=8, type='int', help="Centroiding algorithm", ) parser.add_option("--nsmooth", default=5, type='int', help="Median filter smoothing half-width", ) parser.add_option("--tsample", default=32.8, type='float', help="Sampling time for temperature correction (sec)", ) parser.add_option("--plotfile", default='corr_plot.jpg', help="Output file for diagnostic plot of fid centroid adjustments", ) # args = file names (opt, args) = parser.parse_args() return (opt, args) def get_acis0eng(files): """Read ACIS engineering telemetry opt.acisengfiles and return a pyfits bintable HDU with time and detector housing temperature MSID columns.""" # IMPORTANT NOTE: The current code does not do quality filtering on the # input acis0eng files. Final code should select on good quality for # the acis0eng_col_names columns. # Convert from a glob or .lis file to a list of file names # (Allowing an input glob is useful for testing, but not required in the final tool). if files.startswith('@'): files = [x.strip() for x in open(files[1:]).readlines()] else: files = glob(files) # Read each of the input ACIS L0 eng. FITS files and keep track of rows nrows = 0 hdulist = {} nrow = {} for f in files: hdulist[f] = pyfits.open(f) nrow[f] = hdulist[f][1].header['naxis2'] nrows += nrow[f] # Select the 3 cols we care about about and concatentate all the input data sets # into a new (in-memory) pyfits FITS HDU bintable. [Could probably do this with # a numpy recarray but I found an example in the pyfits docs first]. acis0eng_col_names = ['TIME', '1CBAT', '1CBBT'] acis0eng_cols = [x for x in hdulist[files[0]][1].columns if x.name in acis0eng_col_names] acis0eng_hdu = pyfits.new_table(acis0eng_cols, nrows=nrows) i0 = 0 for f in files: i1 = i0 + nrow[f] for name in acis0eng_col_names: acis0eng_hdu.data.field(name)[i0:i1] = hdulist[f][1].data.field(name) i0 = i1 return acis0eng_hdu def calc_dh_temp(acis0eng, fidprops, acen, opt): """Use empirically determined fid light brightness vs. temperature scaling relation and the coarsely measured ACIS detector housing temperature telemetry to get a good estimate of the DH temperature during observation. Return values:: times : Sample times (sec) [numpy.ndarray] dh_temp : ACIS DH average temperature at 'times' (degC) [numpy.ndarray] """ # Because temperatures vary slowly we can do this processing at a relatively # low sample rate, probably once per 32.8 sec). Set up an array # of sample 'times'. t0 = min(acen.field('time')) - opt.tsample t1 = max(acen.field('time')) + opt.tsample times = np.arange(t0,t1,opt.tsample) # Use only acis0eng values within the aspect interval acis0eng = acis0eng[ np.logical_and(acis0eng.field('time') >= t0, acis0eng.field('time') <= t1) ] # Select only the good fid lights fidprops = fidprops[fidprops.field('id_status') == 'GOOD'] # Select good centroid data with correct algorithm # [Are other status values acceptable here?] acen = acen[np.logical_and(acen.field('status') == 0, acen.field('alg') == opt.alg)] # Create an array for the smoothed counts values fid_counts = np.zeros([fidprops.shape[0], len(times)]) for i_fid, fid in enumerate(fidprops): # Select centroids from the fid slot cen = acen[acen.field('slot') == fid['slot']] n_cen = cen.shape[0] # Determine cen indexes corresponding to the desired output 'times'. # Then for each time stamp calculate the median filtered 'counts' # value. Note that time gaps in centroid values are not an issue since # we just take a certain number of samples backward and forward in time i_cens = cen.field('time').searchsorted(times) for i_time, i_cen in enumerate(i_cens): i0 = i_cen - opt.nsmooth i1 = i_cen + opt.nsmooth if i0 < 0: i0 = 0 if i1 > n_cen: i1 = n_cen fid_counts[i_fid, i_time] = np.median(cen[i0:i1].field('counts')) # Take the average counts value over the GOOD fid lights at each time sample fid_counts = fid_counts.mean(0) fid_counts_avg = fid_counts.mean() # Calc DH temperature (degC) from average of thermistor MSIDs 1CBAT and 1CBBT dh_temp_avg = np.mean(acis0eng.field('1cbat') + acis0eng.field('1cbbt')) / 2.0 - 273.15 # DH temperature calculated by assuming that the average DH temp. from # acis0eng telemetry corresponds to the average fid counts and then using # the empirically determined scaling of degC per fid count. Note that # degC_per_count is negative -- the fid gets brighter when it is cooled. dh_temp = dh_temp_avg + (fid_counts - fid_counts_avg) * degC_per_count return times, dh_temp def calc_fid_center(acen, fidprops, opt): """Calculate fid light center of expansion by shifting the nominal fid center based on the difference between the current observed fid positions and those from the nominal (calibration) fid positions. Return value:: fid_center: Dict with 'y' and 'z' values [degrees] """ fidprops = fidprops[fidprops.field('id_status') == 'GOOD'] n_fid = fidprops.shape[0] # Select good centroid data with correct algorithm acen = acen[np.logical_and(acen.field('status') == 0, acen.field('alg') == opt.alg)] # Calculate the average y and z offsets. y_offset_avg = 0.0 z_offset_avg = 0.0 for i, fid in enumerate(fidprops): cen = acen[acen.field('slot') == fid['slot']] id_num = fid['id_num'] y_offset_avg += cen.field('ang_y_sm').mean() - fid_y_ang_nom[id_num-1] z_offset_avg += cen.field('ang_z_sm').mean() - fid_z_ang_nom[id_num-1] y_offset_avg /= n_fid z_offset_avg /= n_fid # Return the adjusted fid center dict fid_center = {'y': (fid_y_center_nom + y_offset_avg), 'z': (fid_z_center_nom + z_offset_avg)} return fid_center def apply_cent_corr(times, dh_temp, acen, fidprops, opt): """ Remove (subtract) the predicted amount of expansion based on the DH temperature relative to the baseline temperature at which the fid positions were calibrated. The 'ang_y_sm' and 'ang_z_sm' columns of the 'acen' bintable HDU are modified in place. """ # In this routine we do NOT select on status or alg -- all fid centroids get corrected. # Calculate the fid_center values that correspond to the SIM offset for this interval. # The fid_center_nom values were determined using the 2007 CDF-S observations on ACIS-I. fid_center = calc_fid_center(acen, fidprops, opt) for fid in fidprops: ok = np.where(acen.field('slot') == fid['slot']) ang_y_sm = acen.field('ang_y_sm') ang_z_sm = acen.field('ang_z_sm') # Find indexes into 'times' of the acen time stamps. Temps vary slowly so don't worry # about interpolating values, just use picks from searchsorted. i_dh = times.searchsorted(acen.field('time')[ok]) # If making a plot (for TEST ONLY) then copy the smoothed values into the unsmoothed cols if opt.plotfile: acen.field('ang_y')[ok] = ang_y_sm[ok] acen.field('ang_z')[ok] = ang_z_sm[ok] # Make the actual centroid corrections for the y and z axes cor = (dh_temp[i_dh] - dh_temp_base) * fid_CTE ang_y_sm[ok] -= cor * (ang_y_sm[ok] - fid_center['y']) ang_z_sm[ok] -= cor * (ang_z_sm[ok] - fid_center['z']) def make_plot(plotfile, times, dh_temp, acen, fidprops, acis0eng): """ Make a plot of the before (using the unsmoothed centroids) and after centroids for each fid light.""" from matplotlib import pylab as pl t0 = acen[0].field('time') n_fid = fidprops.shape[0] pl.figure(1, figsize=(7,8)) pl.clf() iplot = 1 for fid in fidprops: ok = np.logical_and(acen.field('alg')==8, acen.field('slot')==fid.field('slot')) cen = acen[ok] t = (cen.field('time') - t0)/1000. pl.subplot(n_fid+1, 2, iplot) pl.plot(t, cen.field('ang_y')*3600) pl.plot(t, cen.field('ang_y_sm')*3600, 'r') if iplot==1: pl.title('Y angle (arcsec)') if iplot==5: pl.xlabel('Time (ksec)') pl.ylabel(fid.field('id_string')) pl.subplot(n_fid+1, 2, iplot+1) pl.plot(t, cen.field('ang_z')*3600) pl.plot(t, cen.field('ang_z_sm')*3600, 'r') if iplot==1: pl.title('Z angle (arcsec)') if iplot==5: pl.xlabel('Time (ksec)') iplot += 2 acis_eng_dh_temp = (acis0eng.field('1cbat') + acis0eng.field('1cbbt')) / 2.0 - 273.15 pl.subplot(n_fid+1, 2, iplot) pl.plot((acis0eng.field('time') - t0)/1000, acis_eng_dh_temp) pl.plot((times - t0)/1000, dh_temp) pl.xlabel('Time (ksec)') pl.ylabel('Temp (degC)') pl.savefig(plotfile) pl.show() if __name__ == '__main__': main()