FUNCTION load_events, filename, fileroot=fileroot, t0=t0, texp=texp, $ chan_lim=chan_lim, gti=gti, SILENT=SILENT, $ pi=pi, PATT_LIM=patt_lim ; ---------------------------------------------------------- ;+ ; NAME: ; LOAD_EVENTS ; ; PURPOSE: ; Read XMM EPIC data from filtered events file (FITS) ; ; AUTHOR: ; Simon Vaughan (U.Leicester) ; ; CALLING SEQUENCE: ; t = LOAD_EVENTS("mydata.fits") ; ; INPUTS: ; filename - (string) file name ; ; OPTIONAL INPUTS: ; silent - (logical) switch off output to screen ; chan_lim - (array) set lower and upper PI channel boundaries ; patt_lim - (array) set lower and upper PATTERN ranges ; t_clip - (float) length of time (sec) to 'clip' from the ; end of PN GTI. ; ; OUTPUTS: ; t - vector of event arrival times ; ; OPTIONAL OUTPUTS: ; t0 - (double) time (MET sec) of first event ; texp - (float) Expsure time (from LIVETIME keyword) ; gti - (array) two column table of good start/stop times ; pi - (vector) list of event PI (energy in eV) ; ; PROCEDURES CALLED: ; READFITS, TBGET, SXPAR, FXPOSIT, FXMOVE ; ; DETAILS: ; Read columns of data from a FITS file. The input is a filtered ; XMM EPIC event list. This output is a list of event arrive ; times that can be used to generate a light curve. Optional ; outputs include a list (same length as main output) of the ; event PIs (energy in eV), which can be used to filter on ; evenet energy. ; A GTI (Good Time Interval) table is also generated. This shows ; the start and stop times of each GTI (to the nearest sec). The ; MET (Mission Elapsed Time) in sec of the first event is given ; in T0 - the main output sets the start time to zero, i.e. t[0] ; = 0. The "live" exposure time is also passed via TEXP. ; ; The SILENT keyword switches of the screen output showing the ; data information. ; ; The CHAN_LIM input can be used to specify the lower and upper ; PI channel ranges that are to be collected. Events with PI ; valuess outside this range will be ignored. (The limit is ; inclusive of the end points.) ; ; The PATT_LIM input can be used to specify the lower and upper ; event pattern ranges that are to be collected. Single pixel ; events are PATTERN=0, double pixel events are 1-4, triples are ; 5-8 and quadruples are 9-12. By default we use all PATTERNS in ; the input files (typically 0-12). But to use only singles and ; doubles set PATT_LIM=[0,4], which uses only events for which ; PATTERN is in the range 0-4. TO use only single events set ; PATT_LIM=[0,0]. ; ; T_CLIP is used to curtail the end of the pn exposure. By ; default T_CLIP = 100 sec, which means the GTI table for the pn ; data is modified so that the final interval ends 100s earlier ; than specified in the file. This removes problems caused by ; the exposure ending before the end of the GTI as tabulated, ; which seems to happen (by a few 10's of sec) for the pn (but ; not MOS as far as I know). ; ; (Added 19/01/2011) Automatically collect the GTI table ; from each event file. For the pn this is assumed to be for ; CCDNR=4 (counting 1-12) and for both MOS this is assumed to be ; for CCDNR=1 (counting 1-7). These tables are found in the ; filtered event files. The FXPOSIT routine is used to find the ; extension number in each case. ; ; If no time filtering has been applied there should be a ; STDGTIxx extension in each event file, where xx is the CCDNR ; starting from 01. So STDGTI04 and STDGTI01 are the extensions ; for the aimpoint chips of the pn and MOS, respectively. ; But if time filtering has been applied to the event file there ; will be a second GTI, called GTI0yy05, where yy is the CCDNR ; now starting from 00. So GTI00305 and GTI00005 are the ; extensions for the aimpoint chips of the pn and MOS, ; respectively. Where the latter is available we shall us it, ; since it includes the additional time filtering. (Who decided ; to change the CCD number scheme part-way through the analysis?) ; ; ; Example calls: ; ; IDL> t = load_events("/data/49/sav2/xmm/ngc4051/events/rev0263_pn_src.ds", $ ; chan_lim = [300,1000], t0=t0) ; IDL> plot, histogram(t, binsize=100.0), psym=10 ; HISTORY: ; 23/04/10 - v1.0 - First basic version ; 09/07/10 - v1.1 - Added PATT_LIM keyword for PATTERN ; selection ; 01/10/10 - v1.2 - Added T_CLIP keyword ; 19/01/11 - v1.3 - Automatically find GTI extension ; event files (using FXPOSIT routine) which ; resulted in a large reduction in the amount ; of code needed for the MOS GTI ; collection. Removed redundant GT_FILE ; keyword. Use double precision ; arrays for time series output. ;- ; ---------------------------------------------------------- ; options for compilation (recommended by RSI) COMPILE_OPT idl2 ; watch out for errors on_error, 2 ; ---------------------------------------------------------- ; Check the arguments ; is the file name defined? if (n_elements(filename) eq 0) then begin filename = '' read,'-- Enter file name (ENTER to list current directory): ',filename if (filename eq '') then begin list = findfile() print, list read,'-- Enter file name: ',filename endif endif if NOT KEYWORD_SET(t_clip) then t_clip = 100.0 ; ---------------------------------------------------------- ; call READFITS to read the data from FITS file if KEYWORD_SET(fileroot) then begin file = fileroot+filename endif else begin file = filename endelse filedata = READFITS(file, htab, EXTEN_NO=1, /SILENT) ; check for errors from READFITS if (N_ELEMENTS(filedata) le 1) then begin if (filedata eq -1) then begin print,'** Error calling READFITS in LOAD_EVENTS' print,'** '+!ERROR_STATE.MSG return, -1 endif endif ; now convert three columns - ; TIME (col 1) ; PI (energy; col 9) ; PATTERN (col 11) ; - from binary to floating point time = TBGET(htab, filedata, 1, /NOSCALE) pi = TBGET(htab, filedata, 9) pattern = TBGET(htab, filedata, 11) ; save, and then subtract, the time of first event: T0 t0 = time[0] time = time - t0 ; extract header information target_name = SXPAR(htab, "OBJECT") t_duration = SXPAR(htab, "TELAPSE") t_ontime = SXPAR(htab, "ONTIME") t_exp = SXPAR(htab, "LIVETIME") name_tscope = SXPAR(htab, "TELESCOP") name_inst = strtrim(SXPAR(htab, "INSTRUME")) n_rev = SXPAR(htab, "REVOLUT") datamode = SXPAR(htab, "SUBMODE") pi_name = SXPAR(htab, "OBSERVER") pi_0 = SXPAR(htab, "TLMIN9") pi_1 = SXPAR(htab, "TLMAX9") pi_lim = [pi_0, pi_1] ; ---------------------------------------------------------- ; Find a GTI table in the event file header if (name_inst eq "EPN") then begin gti_name1 = 'STDGTI04' gti_name0 = 'GTI00305' endif else begin gti_name1 = 'STDGTI01' gti_name0 = 'GTI00005' endelse ; use FXPOSIT to find the extension number of the ; GTI table in the fits file. unit = FXPOSIT(file, gti_name0, EXTNUM=extnum, /SILENT, ERRMSG=errmsg) ; if cannot find GTI0yy05 then try STDGTIxx if (unit eq -1) then begin unit = FXPOSIT(file, gti_name1, EXTNUM=extnum, /SILENT, ERRMSG=errmsg) ; not found either GTI if (unit eq -1) then begin print, '** Error using FXPOSIT in LOAD_EVENTS on:',file print, '** ',errmsg return, -1 endif endif ; This is necessary because FXPOSIT does not close the open file ; (unit) even when it returns an error.ro5e_bud FREE_LUN, unit, /FORCE ; call READFITS to read the data from FITS file filedata = READFITS(file, htab, EXTEN_NO=extnum, /SILENT) ; check for errors from READFITS if (N_ELEMENTS(filedata) le 1) then begin if (filedata eq -1) then begin print,'** Error calling READFITS in LOAD_EVENTS' print,'** '+!ERROR_STATE.MSG return, 0 endif endif ; how many rows are there in the file? n_row = SXPAR(htab, "NAXIS2") ; prepare data array for output gti = MAKE_ARRAY(n_row, 2, /DOUBLE) ; now convert two columns - start, stop times - ; from binary to floating point x = TBGET(htab, filedata, 1) gti[0,0] = x[*] x = TBGET(htab, filedata, 2) gti[0,1] = x[*] x = 0 ; subtract zero time T0 gti = gti - t0 ; extract header information t_frame = SXPAR(htab, "FRMTIME") / 1000.0D ; If this is the PN be sure to remove an extra T_CLIP sec from the end of ; the data. For some reason the PN often shuts down before the end of ; the its GTI (01/10/2010) if (name_inst eq "EPN") then begin gti[n_row-1, 1] = MAX([gti[n_row-1, 0], gti[n_row-1, 1] - t_clip]) endif ; ---------------------------------------------------------- ; apply the channel limits if required if KEYWORD_SET(chan_lim) then begin mask = WHERE(pi ge chan_lim[0] and pi le chan_lim[1], count) time = time[mask] pi = pi[mask] pattern = pattern[mask] endif ; ---------------------------------------------------------- ; apply the pattern limits if required if KEYWORD_SET(patt_lim) then begin mask = WHERE(pattern ge patt_lim[0] and pattern le patt_lim[1], count) time = time[mask] pi = pi[mask] endif ; ---------------------------------------------------------- ; print header information IF not KEYWORD_SET(silent) then begin print,"--------------------------------" print,"-- FITS file:",filename print,"-- Observation details:" print,"-- Target: ", target_name print,"-- Mission: ", name_tscope print,"-- Instrument: ", name_inst print,"-- Data mode: ", datamode print,"-- Observer: ", pi_name print,"-- Duration: ", t_duration print,"-- Exposure: ", t_exp print,"-- Revolution: ", n_rev print,"-- Frame time: ", t_frame print,"-- PI range ", pi_lim print,"-- GTI range ", MIN(gti), MAX(gti) print,"--------------------------------" endif ; ---------------------------------------------------------- return, time END