FUNCTION EPIC_LOAD_LC, ROOTPATH=rootpath, DT=dt, N_CHAN=N_chan, OBS_LIST=obs_list, $ DATA_T=data_t, ERROR=error, SEG_LIST=seg_list, T_START=t_start, $ DATA_B=data_b, CHAN_LIST=chan_list, OBS_NAME=obs_name, $ INST_MASK=inst_mask, PATT_LIM=patt_lim, PI_LIST=pi_list, $ SEED=seed, MINF=minf, FRAC=frac, QUIET=quiet ; ---------------------------------------------------- ;+ ; NAME: ; EPIC_LOAD_LC ; ; PURPOSE: ; Load into an array multiple XMM/EPIC time series ; ; AUTHOR: ; Simon Vaughan (U.Leicester) ; ; CALLING SEQUENCE: ; result = EPIC_LOAD_LC() ; ; INPUTS: ; ROOTPATH - (string) path to the event and aux files ; ; OPTIONAL INPUTS: ; DT - (float) time bin size ; N_CHAN - (integer) Number of (log spaced) frequency bands ; OBS_LIST - (vector) which observations to process ; INST_MASK - (vector) use to ignore PN, M1, M2 cameras ; PATT_LIM - (array) set lower and upper PATTERN ranges ; PI_LIST - (array) list of PI energy intervals to use ; SEED - (integer) optional seed for random number generation ; QUIET - (logical) more or less feedback to screen? ; ; OUTPUTS: ; result - (array) combined, net EPIC count rate in ; [N_chan, N_time] ; ; OPTIONAL OUTPUTS: ; DATA_T - (vector) start time of each output bin (sec) ; ERROR - (array) error on each output count rate ; SEG_LIST - (array) [n_seg, 2] array listing first and last ; element numbers of the combined time series ; T_START - (vector) list of the start times of each obs ; DATA_B - (array) background count rate in each ; time/energy bin ; CHAN_LIST - (array) the list of PI channel ranges ; OBS_NAME - (vector) list of the observation names ; FRAC - (array) combined fractional exposure per bin ; ; DETAILS: ; This is essentially a wrapper routine for MAKE_EPIC_LC, which ; loads EPIC data from event files (and auxillary files), ; combined data over all instruments (if specified) correcting ; for GTI losses and background (if specified) with appropriate ; scaling (contained in auxillary files). ; ; ROOTPATH must specify path (directories) where the necessary ; files can be found: event lists and index files. ; Input are the time and energy resolution for the output time ; series. Time bins size DT is in sec. The energy resolution N_CHAN ; is number of logarithmically space channels between PI ; 200-10000 (i.e. 0.2-10 keV). ; ; Output is an array of dimensions [N_chan, N_time] where N_time ; is the total number of time bins and contains the count rates ; in each energy band for each time bin. These are combined over ; all available EPIC cameras (PN, M1, M2) and background ; subtracted and corrected for GTI losses. The actual time of each bin ; is contained in DATA_T, and SEG_LIST records the start:stop ; elements which define each time series. ERROR has the same ; dimensions as RESULT and contains the uncertainties on the ; count rates. See MAKE_EPIC_LC, LOAD_EVENTS and GTI_FIX for ; more details of the extraction and correction procedures. The ; output array T_START list the start times of each observation ; (time for first bin) in absolute (spacecraft seconds) units. ; ; If you wish to ignore a camera completely use the INST_MASK ; vector. By default this is set to [1,1,1] meaning pn, M1 and ; M2 are all used. Set a value to 0 to ignore than camera ; completely. For example [1,0,0] will use only pn data, and ; [0,1,1] will use only M1+M2 data. ; ; The files necessary (contained in the directory ROOTPATH) are ; as follows: ; ; scaling-pn.txt ; scaling-M1.txt ; scaling-M2.txt ; ; These are ASCII files that have a one header line at the top ; and then three columns listing for each of the N_OBS ; observations an identifer name, and the source and background ; region areas (e.g. from the BACKSCAL keywords). For example: ; ; src bkg ; rev1721 2225606 16356161 ; rev1722 2254404 16224966 ; ... ... ... ; ; And for each observation identifier (e.g. "rev1721") we need ; some event files to process. These should be of the form: ; ; __.ds ; ; where is the observation identifier. is the ; instrument ("pn", "M1" or "M2"). is either "src" or ; "bkg" to indicate source or background events. Asuming all ; three cameras are available for "rev1721" we expect to find ; the following files: ; ; rev1721_M1_bkg.ds ; rev1721_M1_src.ds ; rev1721_M2_bkg.ds ; rev1721_M2_src.ds ; rev1721_pn_bkg.ds ; rev1721_pn_src.ds ; ; If there are no data for a given camera we set the area ; scalings (in the "scaling-.txt" file) to zero. This ; will then be skipped over. ; ; For each obsid listed in the "scaling-.txt" files we ; therefore expect 6 files (3 src events, 3 bkg events) assuming ; all three cameras were used. ; ; The vector OBS_LIST can be used to select only a few observations. ; For example, if 5 observations are listed in the index files ; (__.ds) the default value is ; OBS_LIST = [1,2,3,4,5] (note: start counting from 1). ; But if you only wish to take data from observation 4 and 5 set ; OBS_LIST = [4,5]. ; ; The FRAC output keyword shows the fractional exposure per bin ; (0.0-1.0) averaged over the cameras used. ; ; By default EPIC_LOAD_LC prints a one-line summary of each observation ; as it runs. To request more information on each event file set ; QUIET=0 (default QUIET=1) ; ; EXAMPLE USAGE: ; rootpath = '/data/49/sav2/xmm/ngc4051/events/' ; x = EPIC_LOAD_LC(rootpath=ROOTPATH, N_chan=5, data_t=time) ; plot, time, x[0,*], psym=10 ; ; HISTORY: ; 18/06/2010 - v1.0 - first working version ; 23/06/2010 - v1.1 - added OBS_NAME output ; 24/06/2010 - v1.2 - replaced FLAG with SEG_LIST output ; 09/07/2010 - v1.3 - added INST_MASK keyword ; 16/07/2010 - V1.4 - Added PATT_LIM keyword ; 10/11/2010 - v1.5 - Added PI_LIST keyword ; 25/04/2012 - v1.6 - added SEED input/output keyword, for ; controlling randomisation step in ; GTI_FIX, useful for debugging ; Added MINF input keyword and FRAC ; output keyword. Use double precision ; arrays for time series output and ; exposure calculations ; 26/04/2012 - v1.7 - Changed default values for ROOTPATH; ; N_OBS now defines itself based on ; the length of the event index list. ; New OBS input can use used to select ; specific observations. Added QUIET ; keyword ; ; PROCEDURES CALLED: ; ENERGY_BANDS, READ_TABLE, ; MAKE_EPIC_LC [LOAD_EVENTS, GTI_FIX] ; ; NOTES: ;- ; ---------------------------------------------------------- ; options for compilation (recommended by RSI) COMPILE_OPT idl2 ; watch out for errors on_error, 2 ; --------------------------------------------------------- ; Check arguments ; set path root if NOT KEYWORD_SET(rootpath) then begin print,'** ROOTHPATH must be set in call to EPIC_LOAD_LC' print,'** Example: rootpath = "/data/49/sav2/xmm/ngc4051/events/"' print,'** x = EPIC_LOAD_LC(ROOTPATH=rootpath, ...)' RETURN, -1 endif ; set channel number if (N_ELEMENTS(N_chan) ne 1) THEN N_chan = 1 ; set time bin size to default if (N_ELEMENTS(dt) ne 1) THEN dt = 100.0 ; set INST_MASK to default [1,0,0], i.e. only pn if NOT KEYWORD_SET(inst_mask) then inst_mask = [1,0,0] ; set PATT_LIM (PATTERN selection) to default (0-4) if not given on ; input IF (N_ELEMENTS(patt_lim) eq 0) THEN patt_lim = [0, 4] ; check the SEED keyword, if not set assign to !NULL and it will be ; automatically defined when the ranndom number generator is first ; called IF NOT KEYWORD_SET(seed) THEN seed = !NULL ; have we been given minf? if (N_ELEMENTS(minf) ne 1) then minf = 0.3D ; define N_OBS is we have OBS_LIST specified if KEYWORD_SET(obs_list) then N_obs = N_ELEMENTS(obs_list) ; define the level out output to screen if NOT KEYWORD_SET(quiet) then quiet = 1 ; --------------------------------------------------------- ; Begin main routine ; make list of N_CHAN log-spaced energy bands if (N_ELEMENTS(pi_list) eq 0) then begin pi_list = (ENERGY_BANDS(N_chan)) endif xpi_list = TRANSPOSE(pi_list) ; prepare to load the tables listing and src/bkg scaling for ; each instrument N_inst = 3 scale_name = ['scaling-pn.txt', 'scaling-M1.txt', 'scaling-M2.txt'] ; ---------------------------------------------------- ; loop over all instruments and load the src/bkg scaling factors ; load the source and background regions, and calculate the src/bkg ; scaling factor, SCALE. Also load the root name of the dataset, ; OBS_NAME. ; If INST_MASK[i] = 0 then we skip this camera completely. counter = 0 for i = 0, N_inst - 1 do begin if (inst_mask[i] eq 0) then continue scale_data = READ_TABLE(rootpath+scale_name[i], head=1, /TEXT) ; if this is the first index file to be read, use it's size to ; define N_obs if not alredy defined if (counter eq 0) then begin N_eventfiles = (size(scale_data, /DIMENSION))[1] if NOT KEYWORD_SET(N_obs) then N_obs = N_eventfiles if NOT KEYWORD_SET(obs_list) then obs_list = INDGEN(N_obs)+1 scale = MAKE_ARRAY(N_obs, N_inst, /DOUBLE) counter = 1 endif src = FLOAT(scale_data[1, obs_list-1]) bkg = FLOAT(scale_data[2, obs_list-1]) mask = WHERE(src gt 0.0 AND bkg gt 0.0, count) scale[mask, i] = src[mask] / bkg[mask] obs_name = scale_data[0, obs_list-1] endfor ; Trim the array (list of dataset names) to one dimension obs_name = REFORM(obs_name) ; ---------------------------------------------------- ; now load the event files for each dataset ; for each dataset call MAKE_EPIC_LC to produce a merged EPIC time ; series from the relevant event files. Record the start time and ; store the count rate and error in an array. Then concatinate these ; arrays over all datasets. Note that the arrays contain ; [energy, time], i.e. they contain time series in N_chan energy bins ; and N_time time intervals in an array of dimensions [N_chan, ; N_time]. Use the FLAG array to keep track of which elements in the ; output array are from which dataset, by listing the elements of the ; output are the start:stop times of each observation seg_list = MAKE_ARRAY(N_obs, 2, /LONG) t_start = MAKE_ARRAY(N_obs, /DOUBLE) str = '-- EPIC_LOAD_LC processing'+STRING(N_obs)+' observations using PATTERNs:'+STRING(PATT_LIM[0])+"-"+STRING(PATT_LIM[1]) print, STRCOMPRESS(str) for j=0, N_obs-1 do begin scale_j = REFORM(scale[j, *]) file_j = obs_name[j] qinst = '' if (scale_j[0] ne 0) then qinst = qinst + ' pn' if (scale_j[1] ne 0) then qinst = qinst + ' M1' if (scale_j[2] ne 0) then qinst = qinst + ' M2' x = MAKE_EPIC_LC(file_j, scale_j, ROOTPATH=rootpath, T0=t0, DT=dt, /POIS, $ PI_LIST=xpi_list, TIME=time, ERR=err, SILENT=quiet, BKG=bkg, $ PATT_LIM=PATT_LIM, FRAC_I=frac_i, SEED=seed, MINF=minf) ; sum over all energy ranges, and compute mean of this total rate mean_x = MEAN(TOTAL(x, 1)) mean_bkg = MEAN(TOTAL(bkg, 1)) mask = WHERE(inst_mask eq 1) loss = WHERE(frac_i[mask,*] lt minf, count) loss = FLOAT(count)/N_ELEMENTS(frac_i[mask,*]) frac_i = TOTAL(frac_i[mask,*], 1) / N_ELEMENTS(mask) if (j eq 0) then begin print,'-- Name Inst. Duration Start time src rate bkg rate intp frac intp no scale (pn, MOS1, MOS2)' endif print, '-- ', STRTRIM(file_j,2), ' ', STRTRIM(qinst,2), max(time), t0, mean_x, mean_bkg, loss, count, scale_j[*] t_start[j] = t0 N = N_ELEMENTS(time) if (j eq 0) then begin data_x = x data_t = time error = err seg_list[0,*] = [0, N-1] data_b = bkg frac = frac_i endif else begin t_offset = t_start[j] - t_start[0] seg_list[j,*] = [0, N-1] + N_ELEMENTS(data_t) data_x = [[data_x], [x]] data_t = [data_t, time+t_offset] error = [[error], [err]] data_b = [[data_b], [bkg]] frac = [[frac], frac_i] endelse endfor ; --------------------------------------------------------- ; Return to user chan_list = pi_list return, data_x END