corr_fid

1 #!/usr/bin/env python 2 """ 3 Apply a correction to ACIS fid light centroid values to compensate for 4 variations in the ACIS detector housing temperature. 5 6 Inputs 7 --------- 8 :: 9 10 --acisengfiles=ACISENGFILES 11 ACIS L0 engineering files with 1CBAT and 1CBBT 12 --acenfile=ACENFILE Aspect L1 ACEN file 13 --fidpropsfile=FIDPROPSFILE 14 Aspect L1 FIDPROPS file 15 --alg=ALG Centroiding algorithm 16 --nsmooth=NSMOOTH Median filter smoothing half-width 17 --tsample=TSAMPLE Sampling time for temperature correction (sec) 18 19 Outputs 20 --------- 21 :: 22 23 --acenfile=ACENFILE Aspect L1 ACEN file (ang_y_sm and ang_z_sm columns are updated 24 in place) 25 26 Source 27 ------ 28 `Prototype code <http://asc.harvard.edu/mta/ASPECT/corr_fid_cent/corr_fid_cent.py.txt>`_ 29 30 Description 31 ------------ 32 In order to maintain control margin to keep the ACIS focal plane (FP) at the 33 nominal setpoint of -119.7 C, the ACIS detector housing (DH) heater will be 34 disabled. The ACIS fid lights are mounted on the ACIS DH. Allowing the DH 35 temperature to vary has the undesired side effect of producing noticable shifts 36 (up to 0.25 arcsec on time scales of hours) in the ACIS fid light positions due 37 to thermal contraction. For the first 8 years of the mission the DH was controlled 38 at a temperature of -60 C. 39 40 During the nine 2007 Chandra Deep Field South observations the DH heater was 41 disabled at the start and then re-enabled 2 hours before the end of each 42 observation. These tests were performed using a variety of fid light 43 combinations and the shift in fid light positions were fully analyzed. All 44 the data were well-described by a simple physical model of the fid light 45 positions contracting or expanding radially about a single center point with a 46 CTE that corresponds roughly to that of aluminum. The numerical fit values 47 coming out of that analysis are captured in the 'Fid and DH calibration constants' 48 listed below. 49 50 The center of expansion was determined in ACA angular coordinates for ACIS-I 51 observations with the nominal SIM-Z position. The most rigorous (and 52 physically accurate) way to represent these values would be to translate all 53 the derived angular values into physical coordinates (mm) on the ACIS focal 54 plane. Then this correction tool would use the full transformation chain from 55 ACIS focal plane to ACA coordinates to calculate the needed corrections. 56 However, given the small size of the correction and fairly large uncertainties 57 in the derived calibration parameters, we simplify processing considerably by 58 assuming linearity throughout the system and just using offsets with respect to 59 the fid light angular positions used in the CDF-S calibration test 60 observations. These are captured in the fid_[yz]_ang_nom arrays. 61 62 In order to compensate for the temperature-dependent shift in aspect pipeline 63 processing an estimate of the DH temperature is needed. This is available in 64 ACIS engineering telemetry as MSIDs 1CBAT and 1CBBT. Unfortunately the 65 digitization of these values is approximately 2.5 C and using them directly 66 would produce undesirable discontinuities in the aspect solution SIM DY, DZ and 67 DTHETA values. To avoid this we take advantage of the physical characteristic 68 of LEDs that they get brighter as the temperature is lowered. The test data 69 were used to calibrate a linear relationship between changes in DH temperature 70 and fid light brightness (counts). 71 72 Although the fid light brightness are fairly stable over the mission, we cannot 73 depend on stability to the level required for an absolute translation from fid 74 brightness to DH temperature. Instead a hybrid approach is taken -- read in 75 both the ACIS engineering temperature data and the fid data, filter to be using 76 only good quality values, and then match with a simple average over the aspect 77 time interval. At the same time apply the empirically determined scaling 78 factor to convert delta counts into delta degC. This results in a high-resolution 79 estimate of the actual DH temperature during the interval. 80 81 Using the DH temperature and the derived expansion coefficients, the expected 82 radial expansion or contraction is then removed from the observed (smoothed) 83 fid light centroids that are used in constructing the SIM displacement. 84 Correcting the centroids in place has the advantage that downstream tools (in 85 particular V&V) need no modification. The original (unmodified) centroid 86 values are still available in the unsmoothed angular centroids as well as the 87 original CCD pixel centroids. 88 89 Prototype performance 90 --------------------- 91 The plot below shows an example of the original (blue) and corrected (red) 92 centroids for obsid 8595. The key feature to note is a smooth correction 93 during the abrupt DH warmup near the end when the heater was re-enabled during 94 the test. This indicates that the model is doing a good job for a fast 95 temperature shift of about 8 C. The very bottom sub-plot shows the raw ACIS DH 96 temperature telemetry (blue) and the hybrid estimate using fid light counts. 97 98 .. image:: obs8595.jpg 99 100 In other cases the model does not do quite as well and one see residuals of up 101 to about 0.1 arcsec. This is expected to be acceptable in operations because 102 the actual temperature deltas should be less than a few degrees. Nevertheless 103 some additional effort will go into improving the calibration values. 104 105 .. image:: obs8592.jpg 106 107 CALDB update 108 ------------ 109 The constants defined in the prototype code could logically be put in a new 110 binary table HDU in the CALDB PCAD CALALIGN product with an extension name 111 FIDACISDH (or FID_ACIS_DH?). 112 113 References 114 ------------------ 115 - `Initial fid cooling investigations <http://occweb.cfa.harvard.edu/twiki/Aspect/FidLightCooling>`_ 116 - `CDF-S fid cooling test analysis <http://occweb.cfa.harvard.edu/twiki/Aspect/DHcoolingFidCalibration>`_ 117 118 Analysis and the source for this document are in the directory: 119 /proj/sot/ska/analysis/fid_cooling 120 """ 121 __docformat__ = "restructuredtext" 122 123 import pyfits 124 import numpy as np 125 import sys 126 from glob import glob 127 128 # Fid and DH Calibration constants 129 dh_temp_base = -60.0 # Baseline ACIS DH temperature degC 130 degC_per_count = -7.0 / 26000. # Degrees per fid light count (integrated e-) 131 fid_CTE = 16.3e-6 # Fid light mount coeff. of thermal expansion 132 fid_y_center_nom = -65 / 3600. # Fid light Y-center of expansion (deg) 133 fid_z_center_nom = 1395 / 3600. # Fid light Z-center of expansion (deg) 134 135 # Fid light nominal Y and Z positions relative to expansion center (deg) 136 # These were the fid positions for the CDF-S calibration observations used 137 # to derive the fid_[yz]_center_nom values. 138 fid_y_ang_nom = np.array([ 928, -769, 46, 2145, -1820, 391]) / 3600. 139 fid_z_ang_nom = np.array([-837, 845, -970, 1061, 1060, 1703]) / 3600. 140

141 -def main():

142 opt, args = get_options() 143 144 # Read the input data files. 145 acis0eng = get_acis0eng(opt.acisengfiles).data 146 fidprops = pyfits.open(glob(opt.fidpropsfile)[0])[1].data 147 acen_hdus = pyfits.open(glob(opt.acenfile)[0]) 148 acen = acen_hdus[1].data 149 150 # Make sure there is at least one good fid to correct 151 if (fidprops.field('id_status') == 'GOOD').shape[0] == 0: 152 print 'No good fids - skipping corr_fid_cent' 153 sys.exit(0) 154 155 # Calculate detector housing temperature using ACIS eng. telemetry and fid counts 156 times, dh_temp = calc_dh_temp(acis0eng, fidprops, acen, opt) 157 158 # Apply centroid correction to acen values in-place 159 apply_cent_corr(times, dh_temp, acen, fidprops, opt) 160 161 # For development, optionally make a plot showing results of correction 162 if opt.plotfile: 163 make_plot(opt.plotfile, times, dh_temp, acen, fidprops, acis0eng)

164 165 # Write updated values back to original file, OR do whatever 166 # makes sense within pipeline processing. 167

168 -def get_options():

169 """Get program options.""" 170 from optparse import OptionParser 171 global opt, args 172 parser = OptionParser() 173 parser.set_defaults() 174 parser.add_option("--acisengfiles", 175 default='acisf*_2_eng0.fits*', 176 help="ACIS L0 engineering files with 1CBAT and 1CBBT", 177 ) 178 parser.add_option("--acenfile", 179 default='pcadf*_acen1.fits*', 180 help="Aspect L1 ACEN file", 181 ) 182 parser.add_option("--fidpropsfile", 183 default='pcadf*_fidpr1.fits*', 184 help="Aspect L1 FIDPROPS file", 185 ) 186 parser.add_option("--alg", 187 default=8, 188 type='int', 189 help="Centroiding algorithm", 190 ) 191 parser.add_option("--nsmooth", 192 default=5, 193 type='int', 194 help="Median filter smoothing half-width", 195 ) 196 parser.add_option("--tsample", 197 default=32.8, 198 type='float', 199 help="Sampling time for temperature correction (sec)", 200 ) 201 parser.add_option("--plotfile", 202 default='corr_plot.jpg', 203 help="Output file for diagnostic plot of fid centroid adjustments", 204 ) 205 # args = file names 206 (opt, args) = parser.parse_args() 207 208 return (opt, args)

209

210 -def get_acis0eng(files):

211 """Read ACIS engineering telemetry opt.acisengfiles and return a pyfits 212 bintable HDU with time and detector housing temperature MSID columns.""" 213 214 # IMPORTANT NOTE: The current code does not do quality filtering on the 215 # input acis0eng files. Final code should select on good quality for 216 # the acis0eng_col_names columns. 217 218 # Convert from a glob or .lis file to a list of file names 219 # (Allowing an input glob is useful for testing, but not required in the final tool). 220 if files.startswith('@'): 221 files = [x.strip() for x in open(files[1:]).readlines()] 222 else: 223 files = glob(files) 224 225 # Read each of the input ACIS L0 eng. FITS files and keep track of rows 226 nrows = 0 227 hdulist = {} 228 nrow = {} 229 for f in files: 230 hdulist[f] = pyfits.open(f) 231 nrow[f] = hdulist[f][1].header['naxis2'] 232 nrows += nrow[f] 233 234 # Select the 3 cols we care about about and concatentate all the input data sets 235 # into a new (in-memory) pyfits FITS HDU bintable. [Could probably do this with 236 # a numpy recarray but I found an example in the pyfits docs first]. 237 acis0eng_col_names = ['TIME', '1CBAT', '1CBBT'] 238 acis0eng_cols = [x for x in hdulist[files[0]][1].columns if x.name in acis0eng_col_names] 239 acis0eng_hdu = pyfits.new_table(acis0eng_cols, nrows=nrows) 240 i0 = 0 241 for f in files: 242 i1 = i0 + nrow[f] 243 for name in acis0eng_col_names: 244 acis0eng_hdu.data.field(name)[i0:i1] = hdulist[f][1].data.field(name) 245 i0 = i1 246 247 return acis0eng_hdu

248

249 -def calc_dh_temp(acis0eng, fidprops, acen, opt):

250 """Use empirically determined fid light brightness vs. temperature scaling relation 251 and the coarsely measured ACIS detector housing temperature telemetry to get a 252 good estimate of the DH temperature during observation. 253 Return values:: 254 255 times : Sample times (sec) [numpy.ndarray] 256 dh_temp : ACIS DH average temperature at 'times' (degC) [numpy.ndarray] 257 """ 258 259 # Because temperatures vary slowly we can do this processing at a relatively 260 # low sample rate, probably once per 32.8 sec). Set up an array 261 # of sample 'times'. 262 t0 = min(acen.field('time')) - opt.tsample 263 t1 = max(acen.field('time')) + opt.tsample 264 times = np.arange(t0,t1,opt.tsample) 265 266 # Use only acis0eng values within the aspect interval 267 acis0eng = acis0eng[ np.logical_and(acis0eng.field('time') >= t0, 268 acis0eng.field('time') <= t1) ] 269 270 # Select only the good fid lights 271 fidprops = fidprops[fidprops.field('id_status') == 'GOOD'] 272 273 # Select good centroid data with correct algorithm 274 # [Are other status values acceptable here?] 275 acen = acen[np.logical_and(acen.field('status') == 0, acen.field('alg') == opt.alg)] 276 277 # Create an array for the smoothed counts values 278 fid_counts = np.zeros([fidprops.shape[0], len(times)]) 279 280 for i_fid, fid in enumerate(fidprops): 281 # Select centroids from the fid slot 282 cen = acen[acen.field('slot') == fid['slot']] 283 n_cen = cen.shape[0] 284 285 # Determine cen indexes corresponding to the desired output 'times'. 286 # Then for each time stamp calculate the median filtered 'counts' 287 # value. Note that time gaps in centroid values are not an issue since 288 # we just take a certain number of samples backward and forward in time 289 i_cens = cen.field('time').searchsorted(times) 290 for i_time, i_cen in enumerate(i_cens): 291 i0 = i_cen - opt.nsmooth 292 i1 = i_cen + opt.nsmooth 293 if i0 < 0: 294 i0 = 0 295 if i1 > n_cen: 296 i1 = n_cen 297 fid_counts[i_fid, i_time] = np.median(cen[i0:i1].field('counts')) 298 299 # Take the average counts value over the GOOD fid lights at each time sample 300 fid_counts = fid_counts.mean(0) 301 fid_counts_avg = fid_counts.mean() 302 303 # Calc DH temperature (degC) from average of thermistor MSIDs 1CBAT and 1CBBT 304 dh_temp_avg = np.mean(acis0eng.field('1cbat') + acis0eng.field('1cbbt')) / 2.0 - 273.15 305 306 # DH temperature calculated by assuming that the average DH temp. from 307 # acis0eng telemetry corresponds to the average fid counts and then using 308 # the empirically determined scaling of degC per fid count. Note that 309 # degC_per_count is negative -- the fid gets brighter when it is cooled. 310 dh_temp = dh_temp_avg + (fid_counts - fid_counts_avg) * degC_per_count 311 312 return times, dh_temp

313

314 -def calc_fid_center(acen, fidprops, opt):

315 """Calculate fid light center of expansion by shifting the nominal fid center 316 based on the difference between the current observed fid positions and those 317 from the nominal (calibration) fid positions. 318 319 Return value:: 320 321 fid_center: Dict with 'y' and 'z' values [degrees] 322 """ 323 fidprops = fidprops[fidprops.field('id_status') == 'GOOD'] 324 n_fid = fidprops.shape[0] 325 326 # Select good centroid data with correct algorithm 327 acen = acen[np.logical_and(acen.field('status') == 0, acen.field('alg') == opt.alg)] 328 329 # Calculate the average y and z offsets. 330 y_offset_avg = 0.0 331 z_offset_avg = 0.0 332 for i, fid in enumerate(fidprops): 333 cen = acen[acen.field('slot') == fid['slot']] 334 id_num = fid['id_num'] 335 y_offset_avg += cen.field('ang_y_sm').mean() - fid_y_ang_nom[id_num-1] 336 z_offset_avg += cen.field('ang_z_sm').mean() - fid_z_ang_nom[id_num-1] 337 y_offset_avg /= n_fid 338 z_offset_avg /= n_fid 339 340 # Return the adjusted fid center dict 341 fid_center = {'y': (fid_y_center_nom + y_offset_avg), 342 'z': (fid_z_center_nom + z_offset_avg)} 343 return fid_center

344 345

346 -def apply_cent_corr(times, dh_temp, acen, fidprops, opt):

347 """ Remove (subtract) the predicted amount of expansion based on the DH temperature 348 relative to the baseline temperature at which the fid positions were calibrated. 349 The 'ang_y_sm' and 'ang_z_sm' columns of the 'acen' bintable HDU are modified in place. 350 """ 351 # In this routine we do NOT select on status or alg -- all fid centroids get corrected. 352 353 # Calculate the fid_center values that correspond to the SIM offset for this interval. 354 # The fid_center_nom values were determined using the 2007 CDF-S observations on ACIS-I. 355 fid_center = calc_fid_center(acen, fidprops, opt) 356 357 for fid in fidprops: 358 ok = np.where(acen.field('slot') == fid['slot']) 359 ang_y_sm = acen.field('ang_y_sm') 360 ang_z_sm = acen.field('ang_z_sm') 361 362 # Find indexes into 'times' of the acen time stamps. Temps vary slowly so don't worry 363 # about interpolating values, just use picks from searchsorted. 364 i_dh = times.searchsorted(acen.field('time')[ok]) 365 366 # If making a plot (for TEST ONLY) then copy the smoothed values into the unsmoothed cols 367 if opt.plotfile: 368 acen.field('ang_y')[ok] = ang_y_sm[ok] 369 acen.field('ang_z')[ok] = ang_z_sm[ok] 370 371 # Make the actual centroid corrections for the y and z axes 372 cor = (dh_temp[i_dh] - dh_temp_base) * fid_CTE 373 ang_y_sm[ok] -= cor * (ang_y_sm[ok] - fid_center['y']) 374 ang_z_sm[ok] -= cor * (ang_z_sm[ok] - fid_center['z'])

375

376 -def make_plot(plotfile, times, dh_temp, acen, fidprops, acis0eng):

377 """ Make a plot of the before (using the unsmoothed centroids) and after centroids 378 for each fid light.""" 379 from matplotlib import pylab as pl 380 381 t0 = acen[0].field('time') 382 n_fid = fidprops.shape[0] 383 pl.figure(1, figsize=(7,8)) 384 pl.clf() 385 iplot = 1 386 387 for fid in fidprops: 388 ok = np.logical_and(acen.field('alg')==8, acen.field('slot')==fid.field('slot')) 389 cen = acen[ok] 390 t = (cen.field('time') - t0)/1000. 391 392 pl.subplot(n_fid+1, 2, iplot) 393 pl.plot(t, cen.field('ang_y')*3600) 394 pl.plot(t, cen.field('ang_y_sm')*3600, 'r') 395 if iplot==1: 396 pl.title('Y angle (arcsec)') 397 if iplot==5: 398 pl.xlabel('Time (ksec)') 399 pl.ylabel(fid.field('id_string')) 400 401 pl.subplot(n_fid+1, 2, iplot+1) 402 pl.plot(t, cen.field('ang_z')*3600) 403 pl.plot(t, cen.field('ang_z_sm')*3600, 'r') 404 if iplot==1: 405 pl.title('Z angle (arcsec)') 406 if iplot==5: 407 pl.xlabel('Time (ksec)') 408 409 iplot += 2 410 acis_eng_dh_temp = (acis0eng.field('1cbat') + acis0eng.field('1cbbt')) / 2.0 - 273.15 411 pl.subplot(n_fid+1, 2, iplot) 412 pl.plot((acis0eng.field('time') - t0)/1000, acis_eng_dh_temp) 413 pl.plot((times - t0)/1000, dh_temp) 414 pl.xlabel('Time (ksec)') 415 pl.ylabel('Temp (degC)') 416 pl.savefig(plotfile) 417 pl.show()

418 419 420 if __name__ == '__main__': 421 main() 422

Source Code for Module corr_fid_cent