; this idl program extracts pha spectra ; from selected regions of a user provided list ; of Chandra events files ;----------------------------------------------------------------- file = par_get("xtract.par","file",status1) sourcefile = par_get("xtract.par","sourcefile",status1) path = par_get("xtract.par","path",status1) nreg = fix(par_get("xtract.par","nreg",status1)) nevt = fix(par_get("xtract.par","nevt",status1)) grpnum = fix(par_get("xtract.par","grpnum",status1)) nchan = fix(par_get("xtract.par","nchan",status1)) dn2evfile = '/bulk/pkg/xray/idl_lib/pro/acis.dn2ev' split_threshold = 14 response = 'none' arf = 'none' ; define program arrays rebin_spec = dblarr(nchan) rebin_chn = dblarr(nchan) counts = dblarr(nchan) bkg_counts = dblarr(nchan) area = dblarr(nreg) bkg_area = dblarr(nreg) specsum_reg = dblarr(nreg,nchan) bkg_specsum_reg = dblarr(nreg,nchan) src_in_rad = dblarr(nreg) src_out_rad = dblarr(nreg) bkg_in_rad = dblarr(nreg) bkg_out_rad = dblarr(nreg) xref = dblarr(nevt) yref = dblarr(nevt) src_name = strarr(nreg) x_c = dblarr(nreg) y_c = dblarr(nreg) src_in_rad(*) = 0 src_out_rad(*) = 0 bkg_in_rad(*) = 0 bkg_out_rad(*) = 0 x_c(*) = 0 y_c(*) = 0 evtfile='' specfile='' sname='' fl = '' header ='' struc = {TYPES,A:1D,B:1D,C:1D,D:'string'} ; open filelist containing evt files get_lun, lun openr, lun, file point_lun, lun, 0 nl = numlines(file) spec_all = dblarr(nl+1,nreg,nchan) bkg_spec_all = dblarr(nl+1,nreg,nchan) tsum = 0.0 ; readf, lun,header ; ====================== ; Loop around event files for ievt = 0,nevt-1 do begin ; readf, lun, format = '(a30,I6)',evtfile,tm,xr,yr readf, lun,struc tm = struc.A xr = struc.B yr = struc.C evtfile = strtrim(struc.D,2) print,evtfile,tm,xr,yr xref(ievt) = xr yref(ievt) = yr time = tm tsum = tsum + time ; ccd= fix(ccd(0)) evtfile = path + strtrim(evtfile) ; Read in .evt file ReadEvents, evtfile ; Filter and extract info from .evt file no_g255 = 0 ontime=dblarr(3) ontime(1) = 1 ; frt ontime(2) = 0; nexp_thresh ; Now extract the spectrum out of each desired extraction region ; ; Open source file get_lun, lun1 openr, lun1, sourcefile point_lun, lun1, 0 nl = numlines(sourcefile) fl = '' for ireg = 0,nreg-1 do begin readf, lun1,xc,yc,srcinrad,srcoutrad,bkginrad,bkgoutrad,sname src_name(ireg) = sname xc = xc + (xref(ievt) - xref(0)) yc = yc + (yref(ievt) - yref(0)) src_annu =[xc,yc,srcinrad,srcoutrad] bkg_annu =[xc,yc,bkginrad,bkgoutrad] x_c(ireg) = xc y_c(ireg) = yc src_in_rad(ireg) = srcinrad bkg_in_rad(ireg) = bkginrad src_out_rad(ireg) = srcoutrad bkg_out_rad(ireg) = bkgoutrad specfile = 'tmp.pha' src_file = "src.fits" ; Create Source Spectrum ; ---------------------------------------------------------------------------- AcisFilter, specfile, split_threshold, dn2evfile, status,events_in_region, $ ccd_id=[0,1,2,3,4,5,6,7,8,9],AMP=[0,1,2,3], TIME=0,$ G0=1, G1=0, G2=1, G3=1, G4=1, G5=0, G6=1, G7=0, $ NO_G255=no_g255, CIRCLE=0, RECTANGLE=0, $ ONTIME=ontime, USE_QUAL=0, /DN ,ANNULUS=src_annu result = mrdfits(specfile,1,h) chn = result.channel spec = result.counts print,'total COUNTS', total(spec) ; Filtered Area area(ireg) = (!pi)*(((src_out_rad(ireg))^2.) - ((src_in_rad(ireg))^2.)) ; REBIN spec which has 1-4096 CHN from the PSU simulator ; to 1-nchan CHN of "standard" ASC rmf files binfactor = 4096/nchan for i = 0,nchan-1 do begin rebin_spec(i) = total(spec(binfactor*i:binfactor*i+1)) rebin_chn(i) = i+1 endfor ; ; now store each spectrum from each obsid and each region spec_all(ievt,ireg,*) = rebin_spec ; Create Background Spectrum ; ---------------------------------------------------------------------------- AcisFilter, specfile, split_threshold, dn2evfile, status,events_in_region, $ ccd_id=[0,1,2,3,4,5,6,7,8,9],AMP=[0,1,2,3], TIME=0,$ G0=1, G1=0, G2=1, G3=1, G4=1, G5=0, G6=1, G7=0, $ NO_G255=no_g255, CIRCLE=0, RECTANGLE=0, $ ONTIME=ontime, USE_QUAL=0, /DN ,ANNULUS=bkg_annu result = mrdfits(specfile,1,h) chn = result.channel spec = result.counts print,'total background COUNTS', total(spec) ; Filtered Area bkg_area(ireg) = (!pi)*(((bkg_out_rad(ireg))^2.) - ((bkg_in_rad(ireg))^2.)) ; REBIN spec which has 1-4096 CHN from the PSU simulator ; to 1-nchan CHN of "standard" ASC rmf files binfactor = 4096/nchan for i = 0,nchan-1 do begin rebin_spec(i) = total(spec(binfactor*i:binfactor*i+1)) rebin_CHN(i) = i+1 endfor ; ; now store each spectrum from each obsid and each region bkg_spec_all(ievt,ireg,*) = rebin_spec endfor endfor free_lun,lun close,lun ; Calculate summed spectra and corresponding .pha files for i_reg = 0,nreg-1 do begin for i_evt = 0,nevt-1 do begin specsum_reg(i_reg,*) = 0 bkg_specsum_reg(i_reg,*) = 0 endfor endfor for i_reg = 0,nreg-1 do begin for i_evt = 0,nevt-1 do begin specsum_reg(i_reg,*) = specsum_reg(i_reg,*) + spec_all(i_evt,i_reg,*) bkg_specsum_reg(i_reg,*) = bkg_specsum_reg(i_reg,*) + bkg_spec_all(i_evt,i_reg,*) endfor endfor for i_reg = 0,nreg-1 do begin counts(*) = specsum_reg(i_reg,*) bkg_counts(*) = bkg_specsum_reg(i_reg,*) ; Write out XSPEC compatible source file ;-------------------------------------------------------------------- pha_str = replicate({CHANNEL:0.0,COUNTS:0.0},nchan) pha_str.CHANNEL = rebin_chn pha_str.COUNTS = counts day = STRMID(SYSTIME(0), 8, 2) mon = STRMID(SYSTIME(0), 4, 3) year = STRMID(SYSTIME(0), 22, 2) date = strtrim(day) +'-'+ strtrim(mon) +'-'+ strtrim(year) theader = HEADFITS('/home/bondi/chartas/acis/pileup/spec6.pha',exten=1) ; name of source pha files phaname_src = strtrim(src_name(i_reg),2) + '_spec_sum_reg' + strtrim(i_reg,2) + '.pha' phaname_src_gr = strtrim(src_name(i_reg),2) + '_spec_sum_reg' + strtrim(i_reg,2) + '_gr' + strtrim(grpnum,2) + '.pha' ; name of background pha files bkg_phaname_src = 'bkg' + strtrim(i_reg,2) + '.pha' fxaddpar, theader, 'EXTNAME', 'SPECTRUM' fxaddpar, theader, 'HDUCLASS', 'OGIP' fxaddpar, theader, 'HDUCLAS1', 'SPECTRUM' fxaddpar, theader, 'HDUCLAS2', 'TOTAL' fxaddpar, theader, 'HDUCLAS3', 'COUNT' fxaddpar, theader, 'HDUVERS1', '1.1.0' fxaddpar, theader, 'FILENAME', phaname_src fxaddpar, theader, 'DATE', date fxaddpar, theader, 'TELESCOP', 'CHANDRA' fxaddpar, theader, 'INSTRUME', 'ACIS' fxaddpar, theader, 'FILTER', 'none' fxaddpar, theader, 'EXPOSURE', tsum fxaddpar, theader, 'AREASCAL', 1.0000E+00 fxaddpar, theader, 'BACKSCAL', area(i_reg) fxaddpar, theader, 'CORRSCAL', 0.0000E+00 fxaddpar, theader, 'BACKFILE', bkg_phaname_src fxaddpar, theader, 'CORRFILE', 'none' fxaddpar, theader, 'RESPFILE', strtrim(response) fxaddpar, theader, 'ANCRFILE', strtrim(arf) fxaddpar, theader, 'CHANTYPE', 'PHA' fxaddpar, theader, 'DETCHANS', nchan fxaddpar, theader, 'TLMIN1', 1 fxaddpar, theader, 'TLMAX1', nchan fxaddpar, theader, 'STAT_ERR', 0 fxaddpar, theader, 'POISSERR', 'T' fxaddpar, theader, 'SYS_ERR', 0 fxaddpar, theader, 'GROUPING', 0 fxaddpar, theader, 'QUALITY', 0 mwrfits,pha_str,phaname_src,theader,/CREATE ; Write out XSPEC compatible background file ;-------------------------------------------------------------------- bkg_pha_str = replicate({CHANNEL:0.0,COUNTS:0.0},nchan) bkg_pha_str.CHANNEL = rebin_chn bkg_pha_str.COUNTS = bkg_counts day = STRMID(SYSTIME(0), 8, 2) mon = STRMID(SYSTIME(0), 4, 3) year = STRMID(SYSTIME(0), 22, 2) date = strtrim(day) +'-'+ strtrim(mon) +'-'+ strtrim(year) theader = HEADFITS('/home/bondi/chartas/acis/pileup/spec6.pha',exten=1) fxaddpar, theader, 'EXTNAME', 'SPECTRUM' fxaddpar, theader, 'HDUCLASS', 'OGIP' fxaddpar, theader, 'HDUCLAS1', 'SPECTRUM' fxaddpar, theader, 'HDUCLAS2', 'TOTAL' fxaddpar, theader, 'HDUCLAS3', 'COUNT' fxaddpar, theader, 'HDUVERS1', '1.1.0' fxaddpar, theader, 'FILENAME', bkg_phaname_src fxaddpar, theader, 'DATE', date fxaddpar, theader, 'TELESCOP', 'CHANDRA' fxaddpar, theader, 'INSTRUME', 'ACIS' fxaddpar, theader, 'FILTER', 'none' fxaddpar, theader, 'EXPOSURE', tsum fxaddpar, theader, 'AREASCAL', 1.0000E+00 fxaddpar, theader, 'BACKSCAL', bkg_area(i_reg) fxaddpar, theader, 'CORRSCAL', 0.0000E+00 fxaddpar, theader, 'BACKFILE', 'none' fxaddpar, theader, 'CORRFILE', 'none' fxaddpar, theader, 'RESPFILE', strtrim(response) fxaddpar, theader, 'ANCRFILE', strtrim(arf) fxaddpar, theader, 'CHANTYPE', 'PHA' fxaddpar, theader, 'DETCHANS', nchan fxaddpar, theader, 'TLMIN1', 1 fxaddpar, theader, 'TLMAX1', nchan fxaddpar, theader, 'STAT_ERR', 0 fxaddpar, theader, 'POISSERR', 'T' fxaddpar, theader, 'SYS_ERR', 0 fxaddpar, theader, 'GROUPING', 0 fxaddpar, theader, 'QUALITY', 0 mwrfits,bkg_pha_str,bkg_phaname_src,theader,/CREATE ; Group pha data get_lun,lun1 openw,lun1,'group.script' printf,lun1,'#!/bin/csh -f' printf,lun1,'grppha ' + strtrim(phaname_src) + ' ' + strtrim(phaname_src_gr) + ' ' + '"group min ' + strtrim(grpnum,2) +'"' + ' exit' free_lun,lun1 close,lun1 spawn,"rm " + strtrim(phaname_src_gr) spawn,"chmod u+x group.script" spawn,"group.script" endfor end