PRO psu_view, template, NO_XY_BLUR = no_xy_blur, NO_PHA_BLUR=no_pha_blur, $ NO_PLOT = no_plot, BIG_SUCKER = big_sucker, SMALL_DUDE=small_dude, $ PSU = psu, L1 = l1, SAVE_IT = save_it, PROCESS_ALL=process_all, $ FITS_OUT=fits_out ; Procedure to read a PSU event file, convert to quasi-level1 product, ; and create a useful plot (color .gif). ; ; ; Input Arguments/KEYWORDS: ; ; template The trw_id to read/plot: 'H-HAS-EA-2.001' ; It can also be wild-carded: 'H-HAS-PI-*' to do a bunch ; of plots, however only the last is returned in PSU and L1. ; /NO_XY_BLUR Normally the event positions are uniformly ; blurred by +/- 0.5 pixel in making the Level 1 DETX, DETY ; (in mm from nominal aim point); this keyword turns ; off the blurring. ; /NO_PHA_BLUR Likewise, the PHA is blurred by +/-0.5 ADU before ; making ENERGY. ; /NO_PLOT Don't make a plot ; /BIG_SUCKER Use this keyword to make one big ; plot of the colored events in X,Y ; (for John Houck's poster) ; Output KEYWORDS: ; ; PSU This is the structure from the psu FITS file and ; contains tags: ; TIME EXPOSURE CCD_ID AMP_ID DETX DETY PHAS QUALCODE ; or for on-board graded data: ; TIME EXPOSURE CCD_ID AMP_ID DETX DETY PHA GRADE CNRMEAN [CHIPX CHIPY] ; Note: psu.DETX and psu.DETY here are really TDETX and TDETY ; ; L1 This is the structure containing the value-added ; "ASC Level 1" like tags: ; DETX DETY ENERGY GRADE FLTGRADE see below for def.s ; ; /SAVE_IT If this KEYWORD is set the L1 structure is saved ; into an IDL save file in save_dir (below) with the ; name: trw_id_l1.idlsav ; /PROCESS_ALL Generally if a wild-card template gives multiple files ; then the list is printed and execution stops. ; To allow multiple file processing PROCESS_ALL loops ; through all files found. (Useful for a one-time ; creation of .gif and .idlsav file set.) ; /FITS_OUT Create a FITS file for the output data ; ; "ASC Level 1" output includes the following columns ; many are available in the PSU dataset and others (above) ; are created here. ; TIME - same as PSU TIME ; CCD_ID - same as PSU CCD_ID ; EXPNO - same as PSU EXPOSURE ; CHIPX - (future PSU will include this) ; CHIPY - (future PSU will include this) ; TDETX - same as PSU DETX (which will change to TDETX) ; TDETY - same as PSU DETY (which will change to TDETY) ; DETX - NEW: Value is in real mm from nominal aim point ; DETY - NEW: Value is in real mm from nominal aim point ; These are randomized with +/- 0.5 pixel unless ; the /NO_XY_BLUR keyword is used. ; X - --- sky coord, value not created here ; Y - --- sky coord, value not created here ; PHAS - same as PSU PHAS (if available) ; PHA - NEW: calculated from 3x3 islands and simple grading ;or ; PHA - same as PSU PHA ; CORN_PHA - same as PSU CNRMEAN ; PI - (not provided, = FIX(ENERGY*1000.) ) ; ENERGY - NEW: real energy in keV calculated with gain/offsets from PHA ; FLTGRADE - NEW: created from event islands OR psu.GRADE ; GRADE - NEW: Returns grade values of 0, 2, or 7. ; STATUS - same as PSU QUALCODE (?) Note that QUALCODE is used ; to assign events a grade of 7 if appropriate. ; ; ; Tried to make it standalone except for: ; mrdfits from astrolib (might have to type astrolib to IDL ; prompt to setup stuff.) ; pointer to acis_gain.fits is needed for gain values ; asca_of_flight - (~dd/idl/xrcf/) returns ASCA grade array ; given ACIS/FLT grade array ; dd_load_ct - (from ~dd/idl/useful/) loads color table for event plots ; ;--------------------------------------------------------------- ; dd 11/21/97 ; dd 12/4-5/97 Modify to work on non-grating files too (e.g., change plot ; ranges, etc.). ; dd 12/8/97 Modify to read acis_gain.fits to fill the gain and ; offset arrays. Use real ENERGY. (Thanks Joel!) ; dd 5/13/98 Add TIME and PHA to L1 save file, Add FITS_OUT ;--------------------------------------------------------------- ; . ... . . ... ... ; Some "constant" directories: ; psu data at MIT: psu_data_dir = '/nfs/opabinia/d5/PHASE_H_DATA/trw_ids/' ; ACIS gain file: acis_gain_file = '/nfs/spectra/d6/CIP/acis_gain.fits' ; what directory to put the resulting .gif plot... ;;gif_dir = '' ; current directory gif_dir = '/nfs/spectra/d6/ACIS_anal/' ; Dan's ACIS_anal directory ; and IDL save file ;;save_dir = '' ; current directory save_dir = '/nfs/spectra/d6/ACIS_anal/' ; Dan's ACIS_anal directory ; . ... . . ... ... ; Some numerical constants, etc. ; pixel size in mm pix_size = 0.024 ; mm ; split event threshold: split_dn = 13 ; e.g., dn=13 is above threshold ; CCD/amp gain and offsets ; E_eV = gain * PHA + offset ; Use the PSU amp_id and ccd_id to select gain/offset amp_gain = FLTARR(10,4) ; ccd_id(0-9), amp_id(0-3) amp_offset = FLTARR(10,4) ; Fill these from FITS file! gains = mrdfits(acis_gain_file, 1) for il=0,n_elements(gains)-1 do begin ccd = gains(il).ccd_id amp = gains(il).ccdnode amp_gain(ccd, amp) = gains(il).gain ; eV/ADU amp_offset(ccd, amp) = gains(il).offset ; eV end ; nominal aim points for I and S arrays in TDET coordinates: ; ACIS-S (on S3) ; 6.04 mm across S3 from S2-S3 edge: s_aim_tdetx = 3918 + FIX(6.04/pix_size) ; Centered in S3: s_aim_tdety =1703 + 512 ; ACIS-I (on I3) i_aim_tdetx = 4108 + FIX(1.46/pix_size) i_aim_tdety = 4085 - FIX(1.46/pix_size) ; . ... . . ... ... ; Get the file name(s) ; If more than one file matches, list them and return... ; (unless /PROCESS_ALL!) ; get the filename - same as the TRW-ID ; get around the file structure with a FINDFILE with wild cards (*): all_found = findfile(psu_data_dir+'*/*/'+template, COUNT=nfound) if nfound EQ 0 then begin print, 'psu_view: non-existant TRW-ID' RETURN end if nfound GT 1 AND NOT(KEYWORD_SET(PROCESS_ALL)) then begin print, 'psu_view: found '+STRING(nfound,FORMAT='(I4)')+' files:' for ifile = 0, n_elements(all_found)-1 do begin file_name = all_found(ifile) ; get the TRW ID from the filename... pieces = STR_SEP(file_name, '/') trw_id = pieces(n_elements(pieces)-1) print, trw_id end print, 'psu_view: please select one of the above' RETURN end print, '' ; . ... . . ... ... ; Process the file(s) for ifile=0, n_elements(all_found)-1 do begin file_name = all_found(ifile) ; get the TRW ID from the filename... pieces = STR_SEP(file_name, '/') trw_id = pieces(n_elements(pieces)-1) print, trw_id ; Use TRW ID to determine I or S array: s_array = 0 if STRPOS(trw_id,'S') EQ 4 then s_array = 1 if s_array then print, 'ACIS-S' else print, 'ACIS-I' ; read the event file in ; and headers primary = mrdfits(file_name, 0, primary_hdr) psu = mrdfits(file_name, 1, extension_hdr) print, tag_names(psu) ; The file H-HAS-MC-3.010 with ~3.5 M events busts ; my little SPARC 20 - limit it to 2 million events: if trw_id EQ 'H-HAS-MC-3.010' then begin psu = psu(0:1999999) end ; OK now for some quasi Level0 to 1 stuff before plot(s) ; ----------------------------------------------------- ;# Bias subtract: the PSU data are already bias subtracted. ;# Event Grading, Pulse height summation, Attribute flagging ; If the PHAS tag is available then grade them etc. ; otherwise PHA and (ACIS/FLT)GRADE will be available... grade_em = 0 if TOTAL(STRPOS(tag_names(psu),'PHAS') GE 0) GT 0. then begin grade_em = 1 print, ' ... grading PHAS, size(PHAS) = ',size(psu.phas) end ; Define the PHA and energy arrays acis_pha = INTARR(n_elements(psu)) acis_e = FLTARR(n_elements(psu)) ; and the (ASCA) grade array grade = 2 + INTARR(n_elements(psu)) ; and flight grade array: flt_grade = INTARR(n_elements(psu)) ; use while loop so ie can be a LONG if grade_em GE 1 then begin ie = LONG(0) while ie LE n_elements(psu)-1 do begin ; E: sum all the above-split PHA values ; to get the photon energy ; This is not quite right for isolated corners? above = psu(ie).PHAS GE split_dn ; Total PHA: acis_pha(ie) = TOTAL(psu(ie).PHAS * above) ; Flight grade: weights = [1,2,4,8,16,32,64,128] flt_grade(ie) = TOTAL(above(1:8)*weights) ; Use QUALCODE to veto events through grading: if psu(ie).QUALCODE NE 0 then flt_grade(ie) = 255 ie = ie + 1 end end else begin ; on-board grading gives PHA and (ACIS/FLT)GRADE acis_pha = psu.PHA ; Flight grade: flt_grade = psu.GRADE end ; In any case need to assign "ASCA" grades to the events print, ' ... creating GRADE from FLTGRADE ...' grade = asca_of_flight(flt_grade) ;# Gain correction ; convert to keV: if KEYWORD_SET(NO_PHA_BLUR) then begin ie = LONG(0) while ie LE n_elements(psu)-1 do begin acis_e(ie) = 0.001 * ( acis_pha(ie) * $ amp_gain(psu(ie).ccd_id, psu(ie).amp_id) + $ amp_offset(psu(ie).ccd_id, psu(ie).amp_id) ) ie = ie + 1 end end else begin ; Blur values pha_blur = RANDOMU(SEED, n_elements(psu)) - 0.5 ie = LONG(0) while ie LE n_elements(psu)-1 do begin acis_e(ie) = 0.001 * ( (acis_pha(ie)+pha_blur(ie)) * $ amp_gain(psu(ie).ccd_id, psu(ie).amp_id) + $ amp_offset(psu(ie).ccd_id, psu(ie).amp_id) ) ie = ie + 1 end end ;# Coordinate transformation: ; PSU tiled coords have 18 pixel gap in S array (0.43 mm) - good. ; Create floating values: acis_x = FLOAT(psu.DETX) acis_y = FLOAT(psu.DETY) ; Randomize for better event X,Y plots (also simulates dither and aspect ; correction effects...) ; e.g., blur them by a "pixel" if NOT(KEYWORD_SET(NO_XY_BLUR)) then begin acis_x = acis_x + RANDOMU(SEED,n_elements(acis_x)) - 0.5 acis_y = acis_y + RANDOMU(SEED,n_elements(acis_y)) - 0.5 end ; Keep acis_X and acis_Y around for the plotting (below) ; Create the DETX and DETY in mm from nominal aim point: if s_array then begin aim_tdetx = s_aim_tdetx aim_tdety = s_aim_tdety end else begin aim_tdetx = i_aim_tdetx aim_tdety = i_aim_tdety end acis_detx = pix_size * (acis_x - aim_tdetx) acis_dety = pix_size * (acis_y - aim_tdety) ; - - - - - - - - - - - - - - - - - ; end of "Level 1" processing ; Create the L1 structure to return ; L1 has the tags: DETX, DETY, ENERGY, GRADE, FLTGRADE, PHA, TIME one_event = {detx:0.0, dety:0.0, energy:0.0, grade:0, fltgrade:0, $ pha:0, time:DOUBLE(0.0)} L1 = REPLICATE(one_event, n_elements(psu)) ; Fill it L1.detx = acis_detx L1.dety = acis_dety L1.energy = acis_e L1.grade = grade L1.fltgrade = flt_grade L1.pha = acis_pha L1.time = PSU.time If KEYWORD_SET(save_it) then begin SAVE, L1, FILENAME=save_dir+trw_id+'_l1.idlsav' end ; Create a FITS file for world If KEYWORD_SET(FITS_OUT) then begin ; write headers out to a new file fits_file = save_dir+trw_id+'_dd_evt.fits' fxwrite, fits_file, primary_hdr, /NOUPDATE fxbhmake, l1_extension_hdr, n_elements(l1) ; Create columns in the header for all the junk... fxbaddcol, 1, l1_extension_hdr, L1(0).time, 'TIME', 'PSU.TIME' fxbaddcol, 2, l1_extension_hdr, PSU(0).detx, 'TDETX', $ 'Det X in PSU tiled pixels' fxbaddcol, 3, l1_extension_hdr, PSU(0).dety, 'TDETY', $ 'Det Y in PSU tiled pixels' fxbaddcol, 4, l1_extension_hdr, acis_pha(0), 'PHA', $ 'Sum of PHAS above threshold or PSU.PHA' fxbaddcol, 5, l1_extension_hdr, L1(0).energy, 'ENERGY', $ 'Energy from gains, offsets' fxbaddcol, 6, l1_extension_hdr, L1(0).grade, 'GRADE', $ 'ASCA grade' fxbaddcol, 7, l1_extension_hdr, L1(0).fltgrade, 'FLTGRADE', $ 'ACIS flight grade' ; Now make it on disk... fxbcreate, fxunit, fits_file, l1_extension_hdr ; and write the data... print, ' Starting FITS writing...' for irow = 0L, n_elements(L1)-1 do begin FXBWRITE, FXUNIT, L1(irow).time, 1, irow+1 FXBWRITE, FXUNIT, PSU(irow).detx, 2, irow+1 FXBWRITE, FXUNIT, PSU(irow).dety, 3, irow+1 FXBWRITE, FXUNIT, acis_pha(irow), 4, irow+1 FXBWRITE, FXUNIT, L1(irow).energy, 5, irow+1 FXBWRITE, FXUNIT, L1(irow).grade, 6, irow+1 FXBWRITE, FXUNIT, L1(irow).fltgrade, 7, irow+1 end fxbfinish, fxunit close, fxunit free_lun, fxunit print, ' FITS complete: created '+fits_file print, '' end ; of FITS ; All done with L1 ;- - - - - - - - - - - - - - - - - - - - - - - - - ; Make pretty pictures... if NOT(KEYWORD_SET(NO_PLOT)) then begin ; Calc a floating time since beginning of data set ; in seconds ftime = psu.TIME ftime = ftime - MIN(ftime) ; time order the events and arrays ; (in case they are not, e.g. came from qpoe somewhere) time_order = SORT(ftime) ftime = ftime(time_order) psu = psu(time_order) acis_x = acis_x(time_order) acis_y = acis_y(time_order) acis_pha = acis_pha(time_order) acis_e = acis_e(time_order) grade = grade(time_order) flt_grade = flt_grade(time_order) ; PLOTTING section ; Select events with: ; energy above 10 keV and less than 0.3 keV ; ASCA grade of 0,2,3,4, or 6 ok_es = where(acis_e GT 0.3 AND acis_e LT 10.0 AND $ (grade EQ 0 OR grade EQ 2 OR grade EQ 3 $ OR grade EQ 4 OR grade EQ 6), nok) if nok GE 1 then begin ftime = ftime(ok_es) acis_x = acis_x(ok_es) acis_y = acis_y(ok_es) acis_e = acis_e(ok_es) print, ' Total events accepted ( 0.3 keVcolor coding dd_load_ct ; Colors for the events based on energy acis_colors = 3+FIX(252.*((acis_E - 0.3)/(10.0-0.3))) ; Plot the events in physical space with a cross-hair ; at the nominal aim point location if s_array then begin plot, acis_x(pes), acis_y(pes), PSYM=3, BACK=2, COLOR=1, $ XRANGE=[500.,7000.], XTICKS=13, $ YRANGE = [1600.,2800.], YTICKS = 6, /NODATA, CHARSIZE = char_size, $ XTITLE = 'TDETX (pixels)', YTITLE='TDETY (pixels)', $ TITLE = trw_id plots, acis_x(pes), acis_y(pes), COLOR=acis_colors(pes), PSYM=3 print, ' nominal aim point (TDETX,TDETY) = '+$ STRING(s_aim_tdetx)+STRING(s_aim_tdety) oplot, s_aim_tdetx+[0.,0.], s_aim_tdety+[-100.,100.], COLOR=1 oplot, s_aim_tdetx+[-100.,100.], s_aim_tdety+[0.,0.], COLOR=1 end else begin plot, acis_x(pes), acis_y(pes), PSYM=3, BACK=2, COLOR=1, $ XRANGE=[3000.,5200.], XSTYLE=1, $ YRANGE = [3000.,5200.], YSTYLE=1, /NODATA, CHARSIZE = char_size, $ XTITLE = 'TDETX (pixels)', YTITLE='TDETY (pixels)', $ TITLE = trw_id plots, acis_x(pes), acis_y(pes), COLOR=acis_colors(pes), PSYM=3 print, ' nominal aim point (TDETX,TDETY) = '+$ STRING(i_aim_tdetx)+STRING(i_aim_tdety) oplot, i_aim_tdetx+[0.,0.], i_aim_tdety+[-100.,100.], COLOR=1 oplot, i_aim_tdetx+[-100.,100.], i_aim_tdety+[0.,0.], COLOR=1 end if NOT(KEYWORD_SET(BIG_SUCKER) OR KEYWORD_SET(SMALL_DUDE)) then begin ; I-array has last 3 plots smaller if NOT(s_array) then !p.multi = [3,1,6] ; ACIS PHA histogram of all events ; pha_bin_size = 0.01 ; pha_hist = histogram(acis_e, MIN=0., MAX=10.0, BINSIZE=pha_bin_size) ; pha_es = pha_bin_size*indgen(n_elements(pha_hist)) ; keV ; 1/26/98 dd replaced the above const E bins with const log bins: ; log_hist is in ~dd/idl/useful. pha_bin_EdE = 100.0 log_hist, acis_e, 1.+1./pha_bin_EdE, pha_es, pha_hist plot_oo, pha_es, pha_hist, PSYM=10, BACK=2, COLOR=1, $ XRANGE=[0.3, 10.], XSTYLE=1, YRANGE=[1.,(MAX(pha_hist) > 1.E4)], $ TITLE = 'ACIS PHA: .../'+trw_id, CHARSIZE = char_size, $ XTITLE = 'Energy (keV)', YTITLE='counts (bin = E/'+$ STRCOMPRESS(STRING(pha_bin_EdE,FORMAT='(F8.1)'),/REMOVE_ALL)+')' pha_colors = 3+FIX(252.*((pha_es - 0.3)/(10.0-0.3))) plots, pha_es, pha_hist, PSYM=1, COLOR=pha_colors ; Plot the Energy vs X if s_array then begin plot_io, acis_x(pes), acis_e(pes), PSYM=3, BACK=2, COLOR=1, $ XRANGE=[500.,7000.], XTICKS=13, $ YRANGE = [0.3, 10.], /NODATA, CHARSIZE = char_size, $ XTITLE = 'TDETX (pixels)', $ YTITLE='Energy (keV)', YSTYLE=1 plots, acis_x(pes), acis_e(pes), COLOR=acis_colors(pes), PSYM=3 end else begin plot_io, acis_x(pes), acis_e(pes), PSYM=3, BACK=2, COLOR=1, $ XRANGE=[3000.,5200.], XSTYLE=1, $ YRANGE = [0.3, 10.], /NODATA, CHARSIZE = char_size, $ XTITLE = 'TDETX (pixels)', $ YTITLE='Energy (keV)', YSTYLE=1 plots, acis_x(pes), acis_e(pes), COLOR=acis_colors(pes), PSYM=3 end ; Plot the events vs time of arrival - see dropped frames etc! y_rang = [3000.,5200.] if s_array then y_rang = [0.,8000.] plot, fuzzy_time, acis_x, PSYM=3, BACK=2, COLOR=1, /NODATA, $ YRANGE=y_rang, YSTYLE=1, CHARSIZE = char_size, $ XTITLE = 'Time (TIME - MIN(TIME))', YTITLE='TDETX (pixels)' plots, fuzzy_time, acis_x, PSYM=3, COLOR=acis_colors oplot, ftime(min_pes)*[1.,1.], [500.,7500.], COLOR=1 oplot, ftime(max_pes)*[1.,1.], [500.,7500.], COLOR=1 end ; Finish the .gif output wshow,2,iconic=0 image = tvrd() tvlct, red, green, blue, /GET if KEYWORD_SET(BIG_SUCKER) OR KEYWORD_SET(SMALL_DUDE)then begin if KEYWORD_SET(BIG_SUCKER) then begin write_gif, gif_dir+trw_id+'_big.gif', image, red, green, blue end else begin write_gif, gif_dir+trw_id+'_small.gif', image, red, green, blue end end else begin write_gif, gif_dir+trw_id+'_anal.gif', image, red, green, blue end !p.multi = 0 end ; of NO_PLOT ;- - - - - - - - - - - - - - - - - - - - - - - - - end ; of ifile loop over files RETURN END