; name : xdeconv.pro ; author : george chartas ; date : 20 March 2000 ; ; purpose: this idl program extracts images ; from for one or more Chandra observations. ; In the case of two of more observations the images ; are realigned and co-added to produce the total stacked image. ; A maximum-likelihood deconvolution technique, ; using the appropriate simulated PSF, ; is applied to each individual observation. ; In the case of two of more observations the resulting ; deconvolved images are realigned ; and co-added to produce the total deconvolved image. ; ; For each observation, a PSF appropriate for ; the focus , aim-point position and source spectrum is created ; employing the simulation tool MARX v2.2 (Wise et al., 1997). ; Simulated PSF's are binned to a sub-pixel scale of 0.25*0.5". ;----------------------------------------------------------------- @readevents.pro file = par_get("xdeconv.par","file",status1) nevt = fix(par_get("xdeconv.par","nevt",status1)) xce = (par_get("xdeconv.par","xce",status1)) xcw = (par_get("xdeconv.par","xcw",status1)) ycs = (par_get("xdeconv.par","ycs",status1)) ycn = (par_get("xdeconv.par","ycn",status1)) psfsize = (par_get("xdeconv.par","psfsize",status1)) OutputDir = (par_get("xdeconv.par","OutputDir",status1)) GratingType = (par_get("xdeconv.par","GratingType",status1)) DetectorType = (par_get("xdeconv.par","DetectorType",status1)) DetIdeal = (par_get("xdeconv.par","DetIdeal",status1)) SourceFlux = fix(par_get("xdeconv.par","SourceFlux",status1)) SpectrumType = (par_get("xdeconv.par","SpectrumType",status1)) SpectrumFile = (par_get("xdeconv.par","SpectrumFile",status1)) SourceType = (par_get("xdeconv.par","SourceType",status1)) SGaussSigma = (par_get("xdeconv.par","S-GaussSigma",status1)) SBetaCoreRadius = (par_get("xdeconv.par","S-BetaCoreRadius",status1)) SBetaBeta = (par_get("xdeconv.par","S-BetaBeta",status1)) xspecmodel = (par_get("xdeconv.par","xspecmodel",status1)) xspecparam = (par_get("xdeconv.par","xspecparam",status1)) xspecrange = (par_get("xdeconv.par","xspecrange",status1)) xspec_model_file = (par_get("xdeconv.par","xspec_model_file",status1)) niter = (par_get("xdeconv.par","niter",status1)) DitherModel = "NONE" evtfile='' struc = {TYPES,A:1D,B:1D,C:1D,D:1D,E:1D,F:1D,G:'string'} dn2evfile = '/bulk/pkg/xray/idl_lib/pro/acis.dn2ev' split_threshold = 14 tsum = 0 ; open filelist containing evt files get_lun, lun0 openr, lun0, file point_lun, lun0, 0 nl = numlines(file) source_circle_xc = dblarr(nl+1) source_circle_yc = dblarr(nl+1) source_circle_xc(*) = 0 source_circle_yc(*) = 0 ; ====================== ; Loop around event files for ievt = 0,nevt-1 do begin readf, lun0,struc tm = struc.A xc = struc.B yc = struc.C SourceOffsetZ = struc.D SourceOffsetY = struc.E roll = struc.F evtfile = strtrim(struc.G,2) time = tm tsum = tsum + time evtfile = 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 source_circle_xc(ievt) = xc source_circle_yc(ievt) = yc ; region 0 rect xllc = xc - xce xurc = xc + xcw yllc = yc - ycs yurc = yc + ycn rect = [xllc,yllc,xurc,yurc] specfile = 'tmp.pha' src_file="src.fits" 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=rect, $ ONTIME=ontime, USE_QUAL=0, /DN ,ANNULUS=0 run = GetProperty( 'X', DDS_DATA=x ) run = GetProperty( 'Y', DDS_DATA=y ) ; Determine psf's using MARX ; first remove current marx parameter file and cp original one spawn, '/usr/bin/rm marxlog.txt marx.par' spawn, 'cp /bulk/pkg/asc/marx/marx.par .' f='("marx NumRays=",I0, " OutputDir=", A, " GratingType=", A, " DetectorType=", A, " DetIdeal=", A, " SourceFlux=",F0, " SpectrumType=", A, " SpectrumFile=", A, " SourceType=", A, " SourceOffsetZ=",F0, " SourceOffsetY=",F0, " S-GaussSigma=",F0, " S-BetaCoreRadius=",F0, " S-BetaBeta=",F0, " DitherModel=", A, " >> marxlog.txt")' ; CREATE model_pars = str_sep(xspecparam,' ') npars = n_elements(model_pars) model_pars = model_pars*1.d0 file2 = strtrim(xspec_model_file) get_lun, lun openw, lun,file2 printf,lun,strtrim(xspecmodel) for i = 0, npars-1 do begin printf,lun,model_pars(i) endfor printf,lun,'dummyrsp 0.05 10 2048' printf,lun,strtrim(xspecrange) printf,lun,'plot model' printf,lun,'setplot com wdata model' printf,lun,'plot model' printf,lun,'exit' printf,lun,'/' free_lun,lun close,lun spawn ,"rm model.qdp" ; spawn ,"setenv LD_LIBRARY_PATH /usr/local/openwin3/lib:/usr/local/X11R6.1/lib:/opt/SUNWspro/lib ; setenv PGPLOT_TYPE /null ; /bulk/pkg/xanadu.v10/SunOS_5.5_sparc/bin/xspec " + strtrim(xspec_model_file),result spawn ,"setenv PGPLOT_TYPE /null ; /bulk/pkg/lheasoft/SunOS_5.6_sparc/bin/xspec - " + strtrim(xspec_model_file),result spawn ," /bulk/pkg/asc/marx/bin/xspec2marx model.qdp > " + strtrim(SpectrumFile) NumRays = 10000 ; OutputDir = strtrim(OutputDir) GratingType = strtrim(GratingType) DetectorType = strtrim(DetectorType) DetIdeal = strtrim(DetIdeal) SpectrumType = strtrim(SpectrumType) SpectrumFile = strtrim(SpectrumFile) SourceType = strtrim(SourceType) cmd = string(NumRays,OutputDir,GratingType,DetectorType, $ DetIdeal,SourceFlux,SpectrumType,SpectrumFile,SourceType,SourceOffsetZ, $ SourceOffsetY,SGaussSigma,SBetaCoreRadius,SBetaBeta,DitherModel, FORMAT=f ) spawn, 'echo "' + cmd + '" >> marxlog.txt ' ; spawn, cmd spawn, "source /bulk/pkg/asc/ascds.csh ;" + cmd,result xc0 = source_circle_xc(0) yc0 = source_circle_yc(0) if (ievt eq 0) then begin xpos = *x ypos = *y ; create image for current obsid 1000x1000 each pixel is 0.25*0.492*1.01626 arcsec = 0.25*0.5 arcsec im = make_image(xpos,ypos,xbinsize=.25*1.01626,ybinsize=.25*1.01626,xrange=[xc0-xce,xc0+xcw],yrange=[yc0-ycs,yc0+ycn]) ; create psf for current obsid ypsf = read_marx_file(strtrim(Outputdir) + '/xpixel.dat') zpsf = read_marx_file(strtrim(Outputdir) + '/ypixel.dat') ; find centroid of psf ypsf_m = mean(ypsf) zpsf_m = mean(zpsf) psf = make_image(ypsf,zpsf,xbinsize=.25*1.01626,ybinsize=.25*1.01626,xrange=[ypsf_m-psfsize,ypsf_m+psfsize],yrange=[zpsf_m-psfsize,zpsf_m+psfsize]) psf_max = max(psf,ind) psf_index = ind ymax= ind mod n_elements(psf(*,0)) zmax= ind/n_elements(psf(0,*)) psf = psf(ymax-psfsize:ymax+psfsize,zmax-psfsize:zmax+psfsize) surface,psf ; rotate and reflect psf psf = reverse(psf,2) psf = rot(psf,roll,1,CUBIC=-0.5) ; normalize psf norm = total(psf) psf = double(psf/norm) ; Deconvolve image with maximum likelihood deconv = im ; Niter = 20 for i=1,Niter do Max_Likelihood, im, psf, deconv, FT_PSF=psf_ft psf_deconv = psf ; Deconvolve PSF with maximum likelihood for i=1,Niter do Max_Likelihood, psf, psf, psf_deconv, FT_PSF=psf_ft d1 = n_elements(im(*,0)) d2 = n_elements(im(0,*)) projx = dblarr(d1) projy = dblarr(d2) imx = dblarr(d1) imy = dblarr(d2) for i =0,d1-1 do projx(i) = total(deconv(i,*)) for i =0,d2-1 do projy(i) = total(deconv(*,i)) for i =0,d1-1 do imx(i) = total(im(i,*)) for i =0,d2-1 do imy(i) = total(im(*,i)) plot,projx oplot,imx,linestyle=2 ss_im = im ss_deconv = deconv ss_psf_deconv = psf_deconv endif if (ievt gt 0) then begin xtmp = *x - ( source_circle_xc(ievt) - source_circle_xc(0) ) ytmp = *y - ( source_circle_yc(ievt) - source_circle_yc(0) ) xpos = [xpos,xtmp] ypos = [ypos,ytmp] ; create image for current obsid im = make_image(xtmp,ytmp,xbinsize=.25*1.01626,ybinsize=.25*1.01626,xrange=[xc0-xce,xc0+xcw],yrange=[yc0-ycs,yc0+ycn]) ss_im = ss_im + im ; create psf for current obsid ypsf = read_marx_file(strtrim(Outputdir) + '/xpixel.dat') zpsf = read_marx_file(strtrim(Outputdir) + '/ypixel.dat') ; find centroid of psf ypsf_m = mean(ypsf) zpsf_m = mean(zpsf) psf = make_image(ypsf,zpsf,xbinsize=.25*1.01626,ybinsize=.25*1.01626,xrange=[ypsf_m-psfsize,ypsf_m+psfsize],yrange=[zpsf_m-psfsize,zpsf_m+psfsize]) psf_max = max(psf,ind) psf_index = ind ymax= ind mod n_elements(psf(*,0)) zmax= ind/n_elements(psf(0,*)) psf = psf(ymax-psfsize:ymax+psfsize,zmax-psfsize:zmax+psfsize) psf = reverse(psf,2) psf = rot(psf,roll,1,CUBIC=-0.5) surface,psf norm = total(psf) psf = double(psf/norm) ; deconvolve image with maximum likelihood deconv = im for i=1,Niter do Max_Likelihood, im, psf, deconv, FT_PSF=psf_ft ; deconvolve psf with maximum likelihood psf_deconv = psf for i=1,Niter do Max_Likelihood, psf, psf, psf_deconv, FT_PSF=psf_ft ; find centroid of core of deconv image and shift deconv image to match first image ss_deconv = ss_deconv + deconv ss_psf_deconv = ss_psf_deconv + psf_deconv d1 = n_elements(im(*,0)) d2 = n_elements(im(0,*)) projx = dblarr(d1) projy = dblarr(d2) imx = dblarr(d1) imy = dblarr(d2) for i =0,d1-1 do projx(i) = total(deconv(i,*)) for i =0,d2-1 do projy(i) = total(deconv(*,i)) for i =0,d1-1 do imx(i) = total(im(i,*)) for i =0,d2-1 do imy(i) = total(im(*,i)) plot,projx oplot,imx endif endfor free_lun,lun0 close,lun0 function_2d,id,ss_deconv end