;+ ; NAME: ra_mike_pipeline.pro ; ; PURPOSE: ; Run the Mike pipeline (Rosetta Alice) to process engineering level ; data into science level data ; ; CATEGORY: data procesing ; ; CALLING SEQUENCE: ; ra_mike_pipeline, FileList, alldata = alldata, outdir = outdir, Verbose=Verbose, ; deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag, ; wavefile = wavefile, flatfile = flatfile, aefffile = aefffile, ; badflag = badflag, badfile = badfile, linflag = linflag, logfile = logfile, ; kernelfile = kernelfile, loadkernels = loadkernels, movefiles = movefiles ; ; OPTIONAL INPUTS: ; filelist -- if supplied, a string or string array giving the ; filenames of the files to be processed. The string ; may be a valid UNIX wildcard. ; flatfile -- filepath to flatfield file ; aefffile -- filepath to effective area calibration file ; ; KEYWORD PARAMETERS: ; ; deadflag -- Default behavior is to apply deadtime correction to ; all data. If set to 0, deadtime correction is not ; performed. ; ; darkflag -- Default behavior is to apply dark counts correction to ; all data. If set to 0, dark count correction is not ; performed. ; ; flatflag -- Default behavior is to apply flatfield correction to ; all data. If set to 0, flatfield correction is not ; performed. ; ; aeffflag -- Default behavior is to apply effective area ; calibration to all data. If set to 0, the effective ; area calibration is not performed. ; ; badflag -- Set this keyword to apply the bad pixel ; mask. CURRENTLY, THIS KEYWORD IS NOT OPERATIONAL ; ; linflag -- Set this keyword to apply the linearization correction ; to the data. ; ; alldata -- set this keyword to process all eng files in the ; rosetta data directory. ; ; ra_data_dir -- set this keyword to the fully-qualified path to the ; rosetta data directory. If not provided, the ; current directory is used. ; ; outdir -- set this keyword to the directory in which to output ; files. ; ; movefiles -- If this keyword is provided, data are written ; to the /data/sci/YYYY/MM/ directory, where YYYY and MM ; are the appropriate Year and Month, respectively, of the ; data file. ; ; OUTPUTS: The pipeline code produces science level data files ; ; COMMON BLOCKS: ; mike_data -- contains a structure with various flags for processing ; ; SIDE EFFECTS: ; none ; ; LIST OF POSSIBLE DATA QUALITY FLAGS: ; 0: Data are nominal ; 1: Data are missing packets ; 2: Aperture door motion event during exposure ; 4: Stim position/Lyman alpha wavelength correction not available, ; wavelength errors of 5-10 A are possible. ; ; PROCEDURE: ; ; ; ; EXAMPLE: ; ; ; ; MODIFICATION HISTORY: ; The code has been completely rewritten to simplify the Rosetta Alice ; pipeline and significantly improve performance. After this rewrite, (v4) ; THE PIPELINE IS NOT COMPATIBLE WITH PREVIOUS VERSIONS. ; Version Date: ; v4. 2007 Feb 22 AJS ; v5. 2007 Jun 13 AJS Added et column to pixel list table ; v6. 2008 Jan 10 AJS Use HK Event file to check aperture door state ; v7. 2008 Jan 23 AJS Use non-paralyzable deadtime correction ; v8. 2008 Feb 11 AJS Use proper HV-level for dark subtraction ; v9. 2008 Sep 02 AJS - Default to using post PC8 effective ; area. Correction for 3.8 kV level in Commissioning ; and Mars data. New dark count subtraction ; algorithm. Stim pixel wavelength corrections ; v10. 2008 Nov 25 AJS - Change error handling to use Poisson error bars when ; raw image counts are lower than 2500. Use Gaussian ; (sqrt) error bars when counts are greater than ; 2500. Also, run deadtime correction before calculating ; error bars. ; v11. 2009 Feb 08 AJS - Use high-order polynomial wavelength solution, ; RA_AEFF_007.TAB, hybrid model deadtime. Create ; linearized data products. Handle deadtime ; associated errors like the ideal non-paralyzable ; case. Correct Countrate files for dark counts ; v12. 2009 Feb 22 AJS - Update version flag in data filenames and fixed bug ; that prevented Aeff correction on some files. ; v13. 2009 Apr 14 AJS - Fixed bug that prevented the LIN file header from ; being updated. ; v14. 2009 Nov 10 AJS - Update to generalize how wavelength solution is ; calculated and make code more machine independent ; v15. 2009 Dec 14 AJS - Use PC10 Aeff and newly derived corrections for ; other epochs, correct count rate extension ; interval, add wavefile flag ; v16. 2009 Dec 18 AJS - Fixed bug with LIN files, wherein the derived flux ; levels were off by the ratio of the pixel size in ; the "natural" frame to the pixel size in the linear frame. ; v17. 2010 Feb 10 AJS - Update to produce "Version 3" data for all data ; through ESB3 ; v18. 2010 Feb 15 AJS - Fix bug that prevented the headers for the primary ; data unit and error image extension from being ; written ; v19. 2010 Mar 26 AJS - Update ra_mike_add_geometry.pro to use in-flight ; pixel look vectors ; v20. 2010 Mar 27 AJS - Added comment to wavelength extension ; ; v21. 2010 Jul 08 AJS - Changed how effective area correction is applied. ; ; v22. 2010 Jul 13 AJS - Fixed bug that resulted in certain files being ; reported erroneously as darks. Removed correction ; for CRINTVC error that has been fixed in v17 of ; LIMA. ; v23. 2010 Jul 15 AJS - Added optional sky subtraction pro ra_mike_pipeline, infiles, alldata = alldata, outdir = outdir, Verbose=Verbose, $ deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag, $ wavefile = wavefile, flatfile = flatfile, skyflag = skyflag, skyfile = skyfile, skyindex = skyindex, $ badflag = badflag, badfile = badfile, linflag = linflag, logfile = logfile, $ kernelfile = kernelfile, loadkernels = loadkernels, movefiles = movefiles, $ hvflag = hvflag, hvfile38 = hvfile38, hvfile37 = hvfile37, no_linearize = no_linearize, $ version = version COMMON MIKE_DATA, mikeflags, logmessages IF n_params() EQ 0 AND NOT keyword_set(alldata) THEN BEGIN print, 'ra_mike_pipeline, infiles, alldata = alldata, outdir = outdir, Verbose=Verbose, $' print, ' deadflag = deadflag, darkflag = darkflag, flatflag = flatflag, aeffflag = aeffflag, $' print, ' wavefile = wavefile, flatfile = flatfile, skyflag = skyflag, skyfile = skyfile, skyindex = skyindex, $' print, ' badflag = badflag, badfile = badfile, linflag = linflag, logfile = logfile, $' print, ' kernelfile = kernelfile, loadkernels = loadkernels, movefiles = movefiles, $' print, ' hvflag = hvflag, hvfile38 = hvfile38, hvfile37 = hvfile37, no_linearize = no_linearize, $' print, ' version = version' RETURN ENDIF ; This is the calibration version number, not the pipeline version number IF NOT keyword_set(version) THEN version = '3' ELSE version = strtrim(version, 2) ; Default behavior is to apply corrections to the data files IF n_elements(deadflag) EQ 0 THEN deadflag = 1 IF n_elements(darkflag) EQ 0 THEN darkflag = 1 IF n_elements(aeffflag) EQ 0 THEN aeffflag = 1 IF n_elements(hvflag) EQ 0 THEN hvflag = 0 IF n_elements(flatflag) EQ 0 THEN flatflag = 0 ; DEFAULT IS NOW NO FLATFIELDING! IF n_elements(badflag) EQ 0 THEN badflag = 0 IF n_elements(rectflag) EQ 0 THEN rectflag = 0 IF n_elements(skyflag) EQ 0 THEN skyflag = 0 ; What machine are we on? spawn, 'hostname', Host Host = strsplit(host, '.', /extract) ; These directories and file references will need to be modified to ; match the directory structure on the local machine. CASE strlowcase(Host[0]) OF 'argonath': BEGIN calibration_dir = '/Volumes/d0/joel/projects/rosetta/ADA/data/cal/' data_dir = '/Volumes/d0/joel/projects/rosetta/ADA/data/' IF NOT keyword_set(kernelfile) THEN $ kernelfile = '/Volumes/d0/joel/projects/rosetta/ADA/data/spice/rosetta_kernels_argonath.txt' END 'ariel': BEGIN calibration_dir = '/Users/joel/projects/rosetta/ADA/data/cal/' data_dir = '/User/joel/projects/rosetta/ADA/data/' IF NOT keyword_set(kernelfile) THEN $ kernelfile = '/Volumes/d0/joel/projects/rosetta/ADA/data/spice/rosetta_kernels_argonath.txt' END 'rioja': BEGIN calibration_dir = '/Users/steffl/data/ralice/data/cal/' data_dir = '/Users/steffl/data/ralice/data/' IF NOT keyword_set(kernelfile) THEN $ kernelfile = '/Users/steffl/kernels/rosetta/rosetta_kernels.txt' END 'rosetta': BEGIN calibration_dir = '/home/soc/rdat/data/cal/' data_dir = '/home/soc/rdat/data/' IF NOT keyword_set(kernelfile) THEN $ kernelfile = '/home/soc/rdat/data/spice/rosetta_kernels_rosetta.txt' END ELSE: BEGIN print, '% RA_MIKE_PIPELINE: Unknown host name. Unable to determine required filepaths.' print, '% RA_MIKE_PIPELINE: Returning.' return END ENDCASE ; Wavelength Calibration file IF NOT keyword_set(wavefile) THEN wavefile = calibration_dir + $ 'wavelength/RA_WAVE_006.FIT' ; Flatfield Correction file IF NOT keyword_set(flatfile) THEN flatfile = calibration_dir + $ 'flat/RA_FLAT_002.FIT' IF keyword_set(flatflag) THEN flat = readfits(flatfile, flathdr, /silent) ; Sky Background file IF NOT keyword_set(skyfile) THEN skyfile = calibration_dir + $ 'sky/RA_SKY_3.9KV_001.FIT' IF keyword_set(skyflag) THEN BEGIN rdfits_struct, skyfile, sky, /silent skyerr = float(sky.im1) sky = float(sky.im0) IF NOT keyword_set(skyindex) THEN skyindex = [6, 7, 8, 20, 21, 22] ENDIF ;;; ; Set up the various keywords and flags. ; IF n_elements(mikeflags) EQ 0 THEN $ MikeFlags = {Version:'23 [2010 Jul 15]', $ Mode:'', $ Verbose: keyword_set(verbose) ? verbose : 0 , $ Log:1, $ Error:0, $ StartTime:systime(1), $ kernel:0, $ nogeometry:0, $ apdoor:0} Time0 = systime(1) ; check to see that the CSPICE ICY DLM has been registered with IDL ; either at startup (the RSI preferred way of doing this) or via the ; DLM_REGISTER procedure. help, /dlm, /brief, out = dlm_help foo = where(strpos(dlm_help, 'ICY') GE 0, icy_loaded) IF icy_loaded EQ 0 THEN BEGIN mike_log, 'The CSPICE ICY DLM has not been registered with IDL, ' + $ 'so geometrical information cannot be determined.', errnum = 1 mikeflags.nogeometry = 1 ENDIF ; The kernels required for SPICE to work have not been loaded, so load ; them now. This is accomplished via a SPICE metakernel file IF NOT keyword_set(mikeflags.kernel) OR keyword_set(loadkernels) THEN BEGIN CSPICE_FURNSH, kernelfile mikeflags.kernel = 1 ENDIF ; Start the log file mike_log, mess, /start, /init, /time, verbose = 1 ; Get the files to process. ; Use FILE_SEARCH to allow the user to enter wildcards in the filelist string IF keyword_set(alldata) THEN BEGIN filelist = file_search(data_dir + 'eng/', 'ra_*eng.fit') filelist = filelist[where(strpos(filelist, 'ra_temp') EQ -1)] nfiles = n_elements(filelist) ENDIF ELSE BEGIN nfilelist = n_elements(infiles) nfiles = 0 filelist = '' FOR i = 0, nfilelist-1 DO BEGIN IF file_test(infiles[i]) THEN BEGIN filelist_i = infiles[i] nfiles_i = 1 ENDIF ELSE filelist_i = file_search(data_dir + 'eng/', infiles[i], count=nfiles_i) filelist = [filelist, filelist_i] nfiles += nfiles_i ENDFOR IF nfiles GT 0 THEN filelist = filelist[1:*] ENDELSE IF (nfiles EQ 0) THEN BEGIN mike_log, err=2, 'No files to process.' RETURN ENDIF ; Find the (chronologically) first and last of the files that will be ; processed, in order to find the range of data directories spanned by ; the data. Since it is theoretically possible that a data file could ; span multiple months, add one month before and one month after this ; range. file_basenames = file_basename(filelist) file_basenames = file_basenames[sort(file_basenames)] first = file_basenames[0] last = file_basenames[nfiles-1] first_year = fix('20' + strmid(first, 3, 2)) first_month = fix(strmid(first, 5, 2)) - 1 IF first_month LT 1 THEN BEGIN first_year -= 1 first_month = 12 ENDIF last_year = fix('20' + strmid(last, 3, 2)) last_month = fix(strmid(last, 5, 2)) + 1 IF last_month GT 12 THEN BEGIN last_year += 1 last_month = 1 ENDIF ; Find all of the available HK Event files in the time span of the data mike_log, 'Searching for valid HK Event files' hkeventfiles = '' nhkeventfiles = 0l FOR year = first_year, last_year DO BEGIN yearstr = strtrim(year, 2) FOR month = 1, 12 DO BEGIN IF NOT ((year EQ first_year AND month LT first_month) OR $ (year EQ last_year AND month GT last_month)) THEN BEGIN monthstr = month LT 10 ? '0'+strtrim(month,2) : strtrim(month,2) hkevfiles = file_search(data_dir + 'eng/' + yearstr + '/' + monthstr + '/', $ '*evnt_eng.txt', count = nhkevfiles) IF nhkevfiles GT 0 THEN BEGIN hkeventfiles = [hkeventfiles, hkevfiles] nhkeventfiles += nhkevfiles ENDIF ENDIF ENDFOR ENDFOR IF n_elements(hkeventfiles) GT 1 THEN hkeventfiles = hkeventfiles[1:*] ;Initialize variables hkformat = "(F15.3,A24,2F15.3,I2,I6,A7,A60)" line_scetr = 0d line_scetc = '' line_insttime = 0d line_toffset = 0d line_instsync = 0b line_eidsubt = 0l line_ident = '' line_descript = '' ; Apdoor contains flags for the aperture door motion event at time, ; scetc scetc = '' apdoor = 0 valid_hkevtimes = [0d, 0d] ; Loop over the HK event files FOR i = 0, nhkeventfiles-1 DO BEGIN mike_log, 'Now reading ' + hkeventfiles[i] openr, lun, hkeventfiles[i], /get first_time = '' last_time = '' dataflag = 0 WHILE NOT eof(lun) DO BEGIN ; Read in the header info line = '' WHILE dataflag EQ 0 DO BEGIN readf, lun, line IF line EQ 'START DATA' THEN dataflag = 1 ENDWHILE readf, lun, line_Scetr, line_ScetC, line_InstTime, line_TOffset, line_InstSync, $ line_EidSubt, line_Ident, line_Descript, format = hkformat IF first_time EQ '' THEN first_time = line_scetc last_time = line_scetc ; Search for Aperture Door Motion events. If found, determine what ; kind of motion event it is. IF line_ident EQ ' ADMOT:' THEN BEGIN scetc = [scetc, line_scetc] CASE 1 OF strpos(line_descript, 'closed') GE 0 : line_apdoor = 0 strpos(line_descript, 'open') GE 0 : line_apdoor = 1 strpos(line_descript, 'between') GE 0 : line_apdoor = 2 strpos(line_descript, 'error') GE 0 : line_apdoor = 3 ENDCASE apdoor = [apdoor, line_apdoor] ENDIF ENDWHILE ; Get the first and last times covered by this HK event file CSPICE_UTC2ET, first_time, first_et CSPICE_UTC2ET, last_time, last_et valid_hkevtimes = [[valid_hkevtimes], [first_et, last_et]] free_lun, lun ENDFOR IF nhkeventfiles GT 0 THEN BEGIN scetc = scetc[1:*] apdoor = apdoor[1:*] valid_hkevtimes = valid_hkevtimes[*, 1:*] ENDIF et_hkev = dblarr(n_elements(scetc)) FOR i = 0, n_elements(scetc)-1 do BEGIN CSPICE_UTC2ET, scetc[i], et et_hkev[i] = et ENDFOR ; Create the 2D array of row offsets in pixel units and initialize the ; wavelength solution coefficients. rdfits_struct, wavefile, wave, /silent wavesoln = wave.im0 poly2 = wave.im1 poly1 = wave.im2 row_offset = wave.im3 row_offset = (fltarr(1024) + 1.) # row_offset wave_comments = [$ 'Note that the actual wavelength range covered by the Alice UV ', $ 'spectrometer is 700-2050 Angstroms. When a charge pulse strikes the ', $ 'readout anodes, the detector electronics assigns it a 10-bit X value', $ 'and a 5-bit Y-value. The mapping of physical space on the detector to', $ 'data space does not use the full 10-bit X range. Instead, events on the', $ 'detector are mapped into columns 130-900 (approximately) in data space.', $ 'The regions in data space on either side of this range do not have any', $ 'physical meaning, and are conceptually similar to overscan areas found', $ 'in CCD data.'] ; Define the nominal linear wavelength scale xindex = findgen(1024) lin_wavelen = float(poly(xindex, poly1)) lin_wave_image = lin_wavelen # (fltarr(32)+1.) yindex = (fltarr(1024) + 1.) # findgen(32) ; Call poisson_errorbars.pro to get vector of Poisson errors poisson_errors = dblarr(2501) FOR i = 0, 2500 DO poisson_errors[i] = poisson_errorbar(i, /double) poisson_errors = FLOAT(poisson_errors) ;;; ; Loop over the files IF (n_elements(NumIms) NE 0) THEN nfiles = nfiles < NumIms mike_log, string(nfiles)+' files to process', Verb=1 FOR F = 0,(nfiles-1) DO BEGIN Filename = Filelist[F] Filenum = string(F+1,strtrim(nfiles,2),format='(I5,"/",A)') IF keyword_set(verbose) THEN print, '% RA_MIKE_PIPELINE: Now processing file ' + $ file_basename(filename) + ' ' + filenum IF strpos(strupcase(filename), 'SCI') GE 0 THEN BEGIN mike_log, err = 2, 'File has already been processed!' GOTO, error ENDIF oldversion = strmid(filename, 8, 1, /reverse_offset) outfile = repstr(file_basename(filename), oldversion + '_eng', version + '_sci') linfile = repstr(file_basename(filename), oldversion + '_eng', version + '_lin') ; Test for existence of file exists = file_test(filename) IF exists EQ 0 THEN BEGIN mike_log, 'File ' + filename + ' not found.', err = 2, verbose = verbose GOTO, error ENDIF ; Make sure we are reading a FITS file A = bytarr(10) openr, Lun, Filename, /get_lun readu, Lun, A free_lun, Lun IF (string(A) NE 'SIMPLE = ') THEN BEGIN mike_log, err = 2, 'The file: ' + filename + ' is not a valid FITS file' GOTO, error ENDIF ; Read in the data rdfits_struct, filename, data, /silent primary_header = data.hdr0 sxaddpar, primary_header, 'BITPIX', -32, /savecomment exptime = float(sxpar(primary_header, 'EXPTIME')) hvlevel = sxpar(primary_header, 'HVLEVELR') ; Initialize the MIKE header block ra_mike_header_init, filename, outfile, primary_header ; Add the wavefile sxaddpar, primary_header, 'WAVEFILE', file_basename(wavefile), /savecomment ; Add the geometry keywords to the header ra_mike_add_geometry, primary_header ; Are there are missing packets in this file? lostpackets = sxpar(primary_header, 'LOSTPKTS') IF lostpackets GT 0 THEN BEGIN data_quality = sxpar(primary_header, 'DATAQUAL') data_quality += 1 sxaddpar, primary_header, 'DATAQUAL', data_quality, /savecomment ENDIF ; Is the aperture door open? If the door is closed for the whole time, ; the image is a dark, and it's not appropriate to apply dark current, ; effective area and flatfield corrections. Attempt to determine the ; aperture door state from the HK event files. CSPICE_UTC2ET, sxpar(primary_header, 'STRTSCET'), et_start CSPICE_UTC2ET, sxpar(primary_header, 'STOPSCET'), et_stop ; Initialize the open_time and closed_time variables open_time = 0. closed_time = 0. hkevind = where(valid_hkevtimes[0,*] LE et_start AND valid_hkevtimes[1,*] GE et_start AND $ valid_hkevtimes[0,*] LE et_stop AND valid_hkevtimes[1,*] GE et_stop, $ hkevfound) IF hkevfound EQ 1 THEN BEGIN ; The start and stop times of this exposure fall into the time span ; covered by one of the HK event files. Thus, the first aperture door ; motion event prior to the start time should give the initial state ; of the aperture door. hkev_index = max(where(et_hkev-et_start LT 0)) ; The exposure happened before the first door motion event, so the ; door state cannot be determined from the event file. IF hkev_index EQ -1 THEN GOTO, nodoorinfo mikeflags.apdoor = apdoor[hkev_index] ; If the door changed during the exposure, or if the error condition ; was set, set the data quality flag and deal with it appropriately. events = where(et_hkev GE et_start AND et_hkev LT et_stop, nevents) IF nevents GT 0 THEN BEGIN apdoor_events = [mikeflags.apdoor, apdoor[events], apdoor[events[nevents-1]]] event_et = [et_start, et_hkev[events], et_stop] ; Although the door open/close times are of the order of 55 msec., the ; door event sampling rate is usually 1 Hz, so assume if an event ; comes through, that the door position instantly changed at the event ; time. IF mikeflags.apdoor EQ 1 THEN open = 1 ELSE open = 0 FOR i = 0, n_elements(apdoor_events)-2 DO BEGIN IF open EQ 1 THEN BEGIN open_time += event_et[i+1] - event_et[i] IF apdoor_events[i+1] EQ 0 OR apdoor_events[i+1] GT 1 THEN open = 0 ENDIF ELSE BEGIN closed_time += event_et[i+1] - event_et[i] IF apdoor_events[i+1] EQ 1 OR apdoor_events[i+1] EQ 2 THEN open = 1 ENDELSE ENDFOR open_time = float(open_time) closed_time = float(closed_time) ; If there was an aperture door motion error event, all bets are off, ; so don't do any more processing. IF max(apdoor_events GT 2) THEN mikeflags.apdoor = 3 ; Write warnings to the log file and adjust the data quality flag mike_log, err = 1, strtrim(nevents, 2) + $ ' aperture door motion event during exposure!' data_quality = sxpar(primary_header, 'DATAQUAL') data_quality += 2 sxaddpar, primary_header, 'DATAQUAL', data_quality, /savecomment ; We have HK event info and there was no aperture door motion event ; during the exposure. ENDIF ELSE CASE mikeflags.apdoor OF 0: closed_time = exptime 1: open_time = exptime ELSE: ENDCASE sxaddpar, primary_header, 'APDRSRC', 'HK event file', $ ' Source of aperture door status', after = 'DATAQUAL' ENDIF ELSE BEGIN ; We don't have any reliable information from the HK event file, so ; just trust the value in the data header. nodoorinfo:apdoor_state = strupcase(sxpar(primary_header, 'APDOOR')) CASE 1 OF apdoor_state EQ 'CLOSED': BEGIN mikeflags.apdoor = 0 closed_time = exptime END apdoor_state EQ 'OPEN': BEGIN mikeflags.apdoor = 1 open_time = exptime END apdoor_state EQ 'BETWEEN': mikeflags.apdoor = 2 apdoor_state EQ 'ERROR': mikeflags.apdoor = 3 ENDCASE sxaddpar, primary_header, 'APDRSRC', 'Data header', $ ' Source of aperture door status', after = 'DATAQUAL' ENDELSE ; Add the OPENTIME and CLOSTIME keywords to the primary header sxaddpar, primary_header, 'OPENTIME', open_time, ' Time aperture door was open', $ after = 'EXPTIME' sxaddpar, primary_header, 'CLOSTIME', closed_time, ' Time aperture door was closed', $ after = 'OPENTIME' ; If the opened time is less than the total exptime, treat as a science image, ; if the open time is equal to the exptime, it is, by definition, a dark IF closed_time EQ exptime THEN isdark = 1 ELSE isdark = 0 ; Right now skip data that have been binned or windowed. If this ; becomes a popular thing to do, this will have to be revised. IF sxpar(primary_header, 'WILOSPEC') NE 0 OR $ (sxpar(primary_header, 'WIHISPEC') NE 0 AND $ sxpar(primary_header, 'WIHISPEC') NE 1023) OR $ sxpar(primary_header, 'WILOSPAT') NE 0 OR $ (sxpar(primary_header, 'WIHISPAT') NE 0 AND sxpar(primary_header, 'WIHISPAT') NE 31) OR $ (sxpar(primary_header, 'WICOSPEC') NE 0 AND sxpar(primary_header, 'WICOSPEC') NE 1) OR $ (sxpar(primary_header, 'WICOSPAT') NE 0 AND sxpar(primary_header, 'WICOSPAT') NE 1) $ THEN BEGIN mike_log, err = 2, 'Windowed and/or binned data currently not allowed' GOTO, error ENDIF mikeflags.mode = strupcase(sxpar(primary_header, 'ACQMODE')) ; High voltage set point hvsetc = sxpar(primary_header, 'hvsetc') ; Split the pipeline here based on the acquisition mode of the ; observation. CASE mikeflags.mode OF ;===================================== ; HISTOGRAM observations ;===================================== 'HISTOGRAM': BEGIN image = data.im0 ; Create error image. Use the Gaussian approximation where image > 2500 ; counts, otherwise use Poisson errors. error_image = FLTARR(SIZE(image, /DIM)) over2500 = WHERE(image GT 2500, nover2500, COMPLEMENT = under2500, NCOMPLEMENT = nunder2500) IF nover2500 GT 0 THEN error_image[over2500] = SQRT(image[over2500]) IF nunder2500 GT 0 THEN error_image[under2500] = poisson_errors[image[under2500]] mkhdr, error_header, error_image, /image sxaddpar, error_header, 'EXTNAME', 'Errors' ; Create the wavelength image stim_correction = ra_stim_lya_correction(image, primary_header, stimflag = stimflag, $ lyaflag = lyaflag, /inverse) IF n_elements(stim_correction) EQ 1024 THEN BEGIN wavelength_indices = stim_correction wavelength_indices = rebin(wavelength_indices, 1024, 32) mike_log, ' Stim pixel correction = YES' sxaddpar, primary_header, 'WAVEFLAG', 'T', /savecomment ENDIF ELSE wavelength_indices = rebin(findgen(1024), 1024, 32) wavelength_indices_no_offset = wavelength_indices wavelength_indices -= row_offset wave_image = ra_wavelen(wavesoln, poly2, index = wavelength_indices, wavefile = wavefile) ; The format of the images as it comes out of LIMA contains the PHD in ; the first 16 elemnts of the image array. This is unneccesary since ; this information is also contained in the PHD extension, and it also ; messes up automatic image scaling routines, since the values of PHD ; bins will be significantly greater than actual image values, ; reducing the available dynamic range. So, set the PHD portion of the ; image to 0. image[0:15, 0] = 0 ; Dead time correction IF keyword_set(deadflag) THEN BEGIN countrate = sxpar(primary_header, 'EVENTS')/exptime ; If the (digital) countrate is > 50 kHz, the file is probbaly a test pattern ; or something really weird is going on. Either way, we probbaly shouldn't ; process it. IF countrate LT 5e4 THEN BEGIN deadtime_correction = ra_mike_deadtime(countrate) image *= deadtime_correction error_image *= deadtime_correction sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?' ENDIF ELSE BEGIN mike_log, err = 2, 'File has unrealistically high countrate and will not be processed.' mike_log, err = 2, 'This may be because the file is a test pattern and not real science data.' GOTO, error ENDELSE ENDIF ; If the HISTOGRAM image has saturated, then replace these pixel ; values with NaNs (Not A Numbers). saturated = where(image EQ 65535, nsaturated) IF nsaturated GT 0 THEN image[saturated] = !values.f_nan ; Correction for dark counts. IF keyword_set(darkflag) AND isdark EQ 0 THEN BEGIN ra_darksubtract, calibration_dir, primary_header, image, error_image, wavelength_indices_no_offset mike_log, ' dark subtraction = YES' ENDIF ; Subtract off sky background. IF keyword_set(skyflag) AND isdark EQ 0 THEN BEGIN ra_mike_skysubtract, primary_header, image, error_image, sky, skyerr, skyfile, $ wavelength_indices_no_offset, skyindex = skyindex mike_log, ' sky subtraction = YES' ENDIF ; Divide by the time the aperture door was open, or the exposure time ; if this is a dark. IF isdark EQ 0 THEN BEGIN image /= open_time error_image /= open_time ENDIF ELSE BEGIN image /= exptime error_image /= exptime ENDELSE sxaddpar, primary_header, 'BUNIT', 'counts s**-1', /savecomment sxaddpar, error_header, 'BUNIT', 'counts s**-1', $ 'Units of data in the error extension' ; Flatfield Correction IF keyword_set(flatflag) AND isdark EQ 0 THEN BEGIN mike_log, ' flatfield correction = YES' image /= flat error_image /= flat sxaddpar, primary_header, 'FLATFLAG', 'T', ' Flatfield correction applied?' sxaddpar, primary_header, 'FLATFILE', /savecomment, file_basename(flatfile) ENDIF ; Effective Area Calibration IF keyword_set(aeffflag) AND isdark EQ 0 THEN BEGIN aeff_image = ra_aeff_norm(primary_header, wave_image, calibration_dir = calibration_dir) image /= aeff_image error_image /= aeff_image sxaddpar, primary_header, 'BUNIT', 'photons cm**-2 s**-1' sxaddpar, error_header, 'BUNIT', 'photons cm**-2 s**-1' ENDIF ; Now write the FITS file and its extensions sxaddpar, primary_header, 'BZERO', 0, ' Data are 32 bit floating point numbers' writefits, OutFile, Image, primary_header writefits, outfile, error_image, error_header, /app mkhdr, Hwave, wave_image, /image sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image' sxaddpar, Hwave, '','' sxaddpar, Hwave, 'COMMENT', 'Alice Wavelength Calibration 2-D Image' IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $ 'Stim pixel wavelength correction: APPLIED' ELSE BEGIN sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED' sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms' ENDELSE sxaddpar, Hwave, 'EXTNAME','Wavelengths' sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array' sxaddhist, wave_comments, Hwave, /comment writefits, OutFile, Wave_Image, Hwave, /app HDRtmp = data.hdr1 sxaddpar, HDRtmp, 'EXTNAME','Pulse Height Distribution' writefits, OutFile, data.im1, HDRtmp, /app HDRtmp = data.hdr2 sxaddpar, data.hdr2, 'EXTNAME','Count Rate' writefits, OutFile, data.im2, HDRtmp, /app ; Now create the linearized image and write it and its extensions to a FITS file IF NOT keyword_set(no_linearize) THEN BEGIN ; Divide by the natural wavelength scale, such that the units are: ; photons / s / cm^2 / Angstrom. ; This prevents the integrated flux level from changing when you interpolate. dwave_image = wave_image - shift(wave_image,-1,0) image_lin = image / dwave_image error_image_lin = error_image / dwave_image lin_wavelength_indices = fltarr(1024, 32) FOR jj = 0, 31 DO lin_wavelength_indices[*, jj] = interpol(xindex, wave_image[*, jj], lin_wavelen) lin_wavelength_indices_no_offset = rebin(lin_wavelength_indices[*, 15], 1024, 32) lin_normalization = total(image_lin[80:950, *]) ; Interpolate the image onto the new Linear scale image_lin = interpolate(image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5) error_image_lin = interpolate(error_image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5) ; Normalize so that we haven't lost any flux due to interpolation errors lin_normalization /= total(image_lin[80:950, *]) mike_log, ' Image Linearization = YES' sxaddpar, primary_header, 'FILE_OUT', strupcase(linfile), /savecomment sxaddpar, primary_header, 'LINFLAG', 'T', /savecomment sxaddpar, primary_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1' sxaddpar, error_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1' writefits, linFile, Image_lin, primary_header writefits, linfile, error_image_lin, error_header, /app mkhdr, Hwave, lin_wavelen, /image sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image' sxaddpar, Hwave, '','' sxaddpar, Hwave, 'COMMENT', 'R-Alice liearized wavelength solution' IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $ 'Stim pixel wavelength correction: APPLIED' ELSE BEGIN sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED' sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms' ENDELSE sxaddpar, Hwave, 'EXTNAME','Wavelength' sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array' sxaddhist, wave_comments, Hwave, /comment writefits, linFile, lin_wavelen, Hwave, /app HDRtmp = data.hdr1 sxaddpar, HDRtmp, 'EXTNAME','Pulse Height Distribution' writefits, linfile, data.im1, HDRtmp, /app HDRtmp = data.hdr2 sxaddpar, HDRtmp, 'EXTNAME','Count Rate' writefits, linfile, data.im2, HDRtmp, /app ENDIF END ;===================================== 'PIXELLIST': BEGIN ;===================================== image = data.im0 ; Create error image. Use the Gaussian approximation where image > 2500 ; counts, otherwise use Poisson errors. error_image = FLTARR(SIZE(image, /DIM)) over2500 = WHERE(image GT 2500, nover2500, COMPLEMENT = under2500, NCOMPLEMENT = nunder2500) IF nover2500 GT 0 THEN error_image[over2500] = SQRT(image[over2500]) IF nunder2500 GT 0 THEN error_image[under2500] = poisson_errors[image[under2500]] mkhdr, error_header, error_image, /image sxaddpar, error_header, 'EXTNAME', 'Errors' ; Create the wavelength image stim_correction = ra_stim_lya_correction(image, primary_header, stimflag = stimflag, $ lyaflag = lyaflag, /inverse) IF n_elements(stim_correction) EQ 1024 THEN BEGIN wavelength_indices = stim_correction wavelength_indices = rebin(wavelength_indices, 1024, 32) mike_log, ' Stim pixel correction = YES' sxaddpar, primary_header, 'WAVEFLAG', 'T', ' Stim pixel correction applied to wavelength solution?' ENDIF ELSE wavelength_indices = rebin(findgen(1024), 1024, 32) wavelength_indices_no_offset = wavelength_indices wavelength_indices -= row_offset wave_image = ra_wavelen(wavesoln, poly2, index = wavelength_indices, wavefile = wavefile) ; Call the routine deform_pixellist to get information about each event pixtable = pixellist_table(data.im1, wave_image = wave_image) ; Add a column containing the ephemeris time of the photon event to ; the pixel list table. time = pixtable[3,*] * sxpar(data.hdr0, 'TIMEHCKC') CSPICE_UTC2ET, sxpar(data.hdr0, 'STRTSCET'), starttime time += starttime pixtable = [pixtable, time] ; Dead time correction IF keyword_set(deadflag) THEN BEGIN orig_image = image image = fltarr(1024, 32) cntrate = histogram(pixtable[3,*], binsize = 1, reverse_indices = R) hacktime = 2.^(-8 + sxpar(primary_Header,'TIMEHCKR')) deadtime_correction = ra_mike_deadtime(cntrate/hacktime) ; Use the arcane and powerful reverse indices to find the number of ; observed photon events between consecutive timehacks. FOR i = 0l, n_elements(cntrate)-1 DO BEGIN IF r[i] NE r[i+1] THEN BEGIN inbin = r[r[i]:r[i+1]-1] ; There is probably a better way to do this without using a for ; loop,but I can't think of one right now, so just plow on ahead with ; brutre force... FOR j = 0, n_elements(inbin)-1 DO $ image[pixtable[0,inbin[j]], pixtable[1,inbin[j]]] += $ deadtime_correction[i] ENDIF ENDFOR notzero = where(orig_image NE 0, count) IF count GT 0 THEN factor = image[notzero]/orig_image[notzero] sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?' ENDIF ; Correction for dark counts. IF keyword_set(darkflag) AND mikeflags.apdoor THEN BEGIN ra_darksubtract, calibration_dir, primary_header, image, error_image, wavelength_indices_no_offset mike_log, ' dark subtraction = YES' ENDIF ; Subtract off sky background. IF keyword_set(skyflag) AND isdark EQ 0 THEN BEGIN ra_mike_skysubtract, primary_header, image, error_image, sky, skyerr, skyfile, $ wavelength_indices_no_offset, skyindex = skyindex mike_log, ' sky subtraction = YES' ENDIF ; Divide by the time the aperture door was open, or the exposure time ; if this is a dark. IF isdark EQ 0 THEN BEGIN image /= open_time error_image /= open_time ENDIF ELSE BEGIN image /= exptime error_image /= exptime ENDELSE sxaddpar, primary_header, 'BUNIT', 'counts s**-1' sxaddpar, error_header, 'BUNIT', 'counts s**-1', $ 'Units of data in the error extension' ; Flatfield Correction IF keyword_set(flatflag) AND mikeflags.apdoor THEN BEGIN mike_log, ' flatfield correction = YES' image /= flat error_image /= flat sxaddpar, primary_header, 'FLATFLAG', 'T', ' Flatfield correction applied?' sxaddpar, primary_header, 'FLATFILE', /savecomment, file_basename(flatfile) ENDIF ; Effective Area Calibration IF keyword_set(aeffflag) AND isdark EQ 0 THEN BEGIN aeff_image = ra_aeff_norm(primary_header, wave_image, calibration_dir = calibration_dir) image /= aeff_image error_image /= aeff_image sxaddpar, primary_header, 'BUNIT', 'photons cm**-2 s**-1' sxaddpar, error_header, 'BUNIT', 'photons cm**-2 s**-1' ENDIF ; Now write the FITS file and its extensions sxaddpar, primary_header, 'BZERO', 0, ' Data are 32 bit floating point numbers' writefits, OutFile, Image, primary_Header writefits, outfile, error_image, error_header, /app mkhdr, Hwave, wave_image, /image sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image' sxaddpar, Hwave, '','' sxaddpar, Hwave, 'COMMENT', 'Alice Wavelength Calibration 2-D Image' IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $ 'Stim pixel wavelength correction: APPLIED' ELSE BEGIN sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED' sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms' ENDELSE sxaddpar, Hwave, 'EXTNAME','Wavelengths' sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array' sxaddhist, wave_comments, Hwave, /comment writefits, OutFile, Wave_Image, Hwave, /app cols = 5 rows = n_elements(pixtable)/cols ftcreate, cols, rows, tabhdr, tab sxaddpar, tabhdr, 'EXTNAME','Pixel List Event Table' ftput, tabhdr, tab, 'X Position', 0, string(reform(pixtable[0,*]), format = '(I4)') ftput, tabhdr, tab, 'Y Position', 0, string(reform(pixtable[1,*]), format = '(I3)') ftput, tabhdr, tab, 'Wavelength', 0, string(reform(pixtable[2,*]), format = '(F8.2)') ftput, tabhdr, tab, 'Timestep', 0, string(reform(pixtable[3,*]), format = '(I6)') ftput, tabhdr, tab, 'Time (et)', 0, string(reform(pixtable[4,*]), format = '(F16.3)') sxaddpar, tabhdr, 'TFORM1', 'I4 ', /savecomment sxaddpar, tabhdr, 'TFORM2', 'I3 ', /savecomment sxaddpar, tabhdr, 'TFORM3', 'F8.2 ', /savecomment sxaddpar, tabhdr, 'TFORM4', 'I6 ', /savecomment sxaddpar, tabhdr, 'TFORM5', 'I ', /savecomment sxaddpar, tabhdr, 'TUNIT1', 'pixels', after = 'TFORM1', /savecomment sxaddpar, tabhdr, 'TUNIT2', 'pixels', after = 'TFORM2', /savecomment sxaddpar, tabhdr, 'TUNIT3', 'Angstroms', after = 'TFORM3', /savecomment sxaddpar, tabhdr, 'TUNIT4', 'time steps', after = 'TFORM4', /savecomment sxaddpar, tabhdr, 'TUNIT5', 'seconds past J2000', after = 'TFORM5', /savecomment ; Add a pair to the end of each ascii table row. These are ; "stealth" characters in that they are not described by the FITS ; header as being in columns. cr = bytarr(1, rows) + 13b lf = bytarr(1, rows) + 10b tab = [tab, cr, lf] sxaddpar, tabhdr, 'NAXIS1', n_elements(tab[*,0]) , /savecomment writefits, outfile, tab, tabhdr, /app HDRtmp = data.hdr2 sxaddpar, HDRtmp, 'EXTNAME','Count Rate' writefits, OutFile, data.im2, HDRtmp, /app ; Now create the linearized image and write it and its extensions to a FITS file IF NOT keyword_set(no_linearize) THEN BEGIN ; Divide by the natural wavelength scale, such that the units are counts / ; angstrom. This prevents the integrated flux level from changing when you interpolate. dwave_image = wave_image - shift(wave_image,-1,0) image_lin = image / dwave_image error_image_lin = error_image / dwave_image lin_wavelength_indices = fltarr(1024, 32) FOR jj = 0, 31 DO lin_wavelength_indices[*, jj] = interpol(xindex, wave_image[*, jj], lin_wavelen) lin_wavelength_indices_no_offset = rebin(lin_wavelength_indices[*, 15], 1024, 32) lin_normalization = total(image_lin[80:950, *]) image_lin = interpolate(image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5) error_image_lin = interpolate(error_image/dwave_image, lin_wavelength_indices, yindex, cubic = -0.5) ; Normalize so that we haven't lost any flux due to interpolation errors lin_normalization /= total(image_lin[80:950, *]) mike_log, ' Image Linearization = YES' sxaddpar, primary_header, 'FILE_OUT', strupcase(linfile), /savecomment sxaddpar, primary_header, 'LINFLAG', 'T', /savecomment sxaddpar, primary_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1' sxaddpar, error_header, 'BUNIT', 'photons cm^-2 s^-1 Angstrom^-1' writefits, linfile, Image_lin, primary_header writefits, linfile, error_image_lin, error_header, /app mkhdr, Hwave, lin_wavelen, /image sxaddpar, Hwave, 'XTENSION', 'IMAGE', 'Image extension: Wavelength Image' sxaddpar, Hwave, '','' sxaddpar, Hwave, 'COMMENT', 'R-Alice liearized wavelength solution' IF n_elements(stimcorrection) EQ 1024 THEN sxaddpar, Hwave, 'COMMENT', $ 'Stim pixel wavelength correction: APPLIED' ELSE BEGIN sxaddpar, hwave, 'COMMENT', 'Stim pixel wavelength correction: NOT APPLIED' sxaddpar, hwave, 'COMMENT', 'wavelength solution may have innaccuracies of 5-10 Angstoms' ENDELSE sxaddpar, Hwave, 'EXTNAME','Wavelength' sxaddpar, Hwave, 'BUNIT', 'Angstroms', 'Units for the wavelength array' sxaddhist, wave_comments, Hwave, /comment writefits, linfile, lin_wavelen, Hwave, /app writefits, linfile, tab, tabhdr, /app HDRtmp = data.hdr2 sxaddpar, HDRtmp, 'EXTNAME','Count Rate' writefits, linfile, data.im2, HDRtmp, /app ENDIF END ;===================================== 'COUNTRATE':BEGIN ;===================================== cnt = data.im0 cnterr = sqrt(cnt) mkhdr, error_header, cnterr, /image sxaddpar, error_header, 'EXTNAME','Errors' ; Divide by the exposure time cntinterval = sxpar(primary_header, 'crintvc') cnt /= cntinterval cnterr /= cntinterval sxaddpar, primary_header, 'BUNIT', 'counts s**-1', /savecomment sxaddpar, error_header, 'BUNIT', 'counts s**-1', 'Units of data in the error extension' ; Dead time correction IF keyword_set(deadflag) THEN BEGIN deadtime_correction = ra_mike_deadtime(cnt) cnt *= deadtime_correction cnterr *= deadtime_correction sxaddpar, primary_header, 'DEADFLAG', 'T', ' Dead-time correction applied?' ENDIF ; Correction for dark counts. IF keyword_set(darkflag) AND isdark EQ 0 THEN BEGIN ra_darksubtract, calibration_dir, primary_header, cnt mike_log, ' dark subtraction = YES' ENDIF sxaddpar, primary_header, 'EXTEND', 'T', ' file may contain extensions', $ before = 'ORIGIN' sxaddpar, primary_header, 'BZERO', 0, ' Data are 32 bit floating point numbers' writefits, OutFile, cnt, primary_Header writefits, outfile, cnterr, error_header, /append END ENDCASE mike_log, 'Processed data written to file: ' + OutFile IF keyword_set(movefiles) THEN BEGIN year = '20' + strmid(filename, strpos(filename, 'ra_') + 3, 2) mon = strmid(filename, strpos(filename, 'ra_') + 5, 2) outdir = data_dir + 'sci/' + year + '/' + mon + '/' outdir_lin = repstr(outdir, 'sci', 'lin') ENDIF IF keyword_set(outdir) THEN BEGIN file_move, outfile, outdir, /overwrite IF NOT keyword_set(no_linearize) AND mikeflags.mode NE 'COUNTRATE' THEN file_move, linfile, outdir_lin, /overwrite ENDIF ; The actual error handling routine.If we ever get here, some error ; occurred and the calibrated file could not be created. ; STILL NEEDS WORK! ; catch, error_status ; IF error_status NE 0 THEN BEGIN ; help, calls = calls ; mike_log, err = 2, '%% Fatal Error at ' + calls ; mike_log, !error_state.msg, errnum = 2 ; error: mike_log, 'File ' + filename + ' not processed!', errnum = 2 ; message, /reset_error_state ; catch, /cancel ; ENDIF error:foo = 0 ENDFOR dT = systime(1) - Time0 dT = strtrim(string(dT / 60., format='(F20.2)'),2) mike_log, 'Total processing time = ' + dT + ' minutes', verb=1 IF keyword_set(logfile) THEN BEGIN openw, lun, logfile, /get printf, lun, logmessages free_lun, lun ENDIF END