;+ ; NAME: clean.pro ; function clean,dirty_map,mod_pat,params,count_rate,TAPER=weights ; ; PURPOSE: ; To extract a "clean" map from a "dirty" map ; ; METHOD: ; Iteratively decomposes the image into single bright pixels, ; subtracting out the side lobes for each point. ; ; Clean does not know about collimators, pitches or angles. It takes ; the modulation patterns as they are given. It cannot tell which ; mod_pats are from which collimator, so weighting (aka tapering) ; must be done via a vector which "knows" which angles and ; collimators are which. ; ; INPUTS: ; Dirty_map = back-projection map created from binned or unbinned ; photon list--must be 1-dimensional of size N^2 ; mod_pat_str = structure containing ; mod_pat, the same array used to create dirty_map, ; = N^2 x M matrix of modulation patterns, where ; N is size of map ; M is number of time bins or photons ; ncoll = vector giving the collimators used in mod_pat ; nbins = vector of sizes of each array within mod_pat ; For example, if mod_pat includes collimators in the ; order [3,8,5,2], that is vector ncoll, and if there ; are 1000 bins for collimators 3 and 2, and 100 for 5 and 8, ; nbins=[1000,100,100,1000] ; ; params = structure created by clean_params.pro containing: ; ; niter = number of iterations to perform ; frac = amount of PSF to subtract in each iteration ; chi_sq_crit = max allowed chi_sq ; ; OPTIONAL INPUTS: ; count_rate = the count-rate profile of size M used for making ; the dirty map. The model count rate will be put in its place. ; If it is not defined, then no model counts will be computed. ; TAPER=weights, where weights is a vector of size M giving the ; factors by which each mod_pat is weighted for cleaning and ; modeling. ; ; OUTPUTS; ; A structure containing: ; clean_map ; clean_components ; clean_amplitudes ; clean_counts ; params ; chi-squared (if count_rate is defined.) ;- ; ALGORITHM: ; 1. Copy dirty_map to map, set N=0 ; 2. Find the brightest point (x,y) in map, if negative, quit. ; 3. Look at component list to see if (x,y) was done before ; a. If it was, extract the PSF from cache ; b. If not, then use backprojection to create the PSF, add to cache ; 4. Subtract PSF*frac from map, set N=N+1, save (x,y), increm clean_map ; 5. If count_rate exists, compute modulation profile and chi_sq ; 6. If chi_sq > chi_sq_crit or N < Niter, go to 2; ; 7. Add up all the clean components, replace count_rate with the ; model count_rate and return. ; ; ; function countrate_profile,map,mod_pat ; Computes modulation profile (count rate) for a map, using the mod_pat ; INPUTS: ; map must have been reformed to 1-D (like 1st index of mod_pat). ; mod_pat = N^2 x M matrix of modulation patterns profile = map # mod_pat return,profile end function back_proj, xy, mod_pat, mean_sq ; ; Computes the back-projection map for a point source ; Used to make a PSF for CLEAN ; Mean_sq is a vector of mean squares, one for each mod_pat map ; Nominally, this is 1/12 for all, but in practice there are variations ; of a few percent even for slit=slat. ; a = reform(mod_pat(xy,*)) ; this is a 1-d profile norm=mean_sq*n_elements(a) b_proj = mod_pat # (a/norm) return,b_proj end function chi_square,model_prof,obs_prof ; RETURNS A VECTOR OF CHI-SQUARED COMPONENTS chisq_vector=(model_prof-obs_prof)^2/model_prof return,chisq_vector end ;CLEAN PROGRAM function clean,dirty_map,mod_pat_str,params,f_obs=count_rate,TAPER=weights niter=params.(0) frac=params.(1) chi_sq_crit=params.(2) mod_pat=mod_pat_str.(0) he=hessi_params() s2s=he.(3) s2p=s2s/(1.+s2s) ; slit/pitch ratio mean_sq=mod_pat_str.(3) ; 1. Copy dirty_map to map, initialize stuff NXY=n_elements(dirty_map) NX=sqrt(NXY) map=reform(dirty_map,NXY) clean_map=map*0 iter=0 Ncache=0 ncomp=0 xy_cache=[-1] psf_cache=map chi_sq_tot=1.e20 map_max=1. clean_flux=0 obs_flux=1.e20 label=' ITER NCACHE XY(0) MAP_MAX TOTAL(CHI_SQ) TOTAL(CLEAN_COUNTS) if keyword_set(count_rate) then begin obs_flux=total(count_rate) print,'Observed flux=',obs_flux endif while ((chi_sq_tot GT chi_sq_crit) AND (iter LT niter) AND (map_max GT 0) $ AND (clean_flux LT obs_flux)) do begin ; 2. Find the brightest point (x,y) in map, if negative, then exit loop max_pt: map_max=max(map) if (map_max GT 0) then begin xy=where(map EQ map_max) xy=xy(0) ; must make it a scalar, selecting only first xy ; 3. Look at xy_cache to see if (x,y) was done before index=where((xy_cache EQ xy),n_before) ; a. If it was, extract the PSF from cache if (n_before GT 0) then begin psf=psf_cache(*,index) endif else begin ; b. If not, then use backprojection to create the PSF, add to cache psf=back_proj(xy, mod_pat, mean_sq) psf=psf/max(psf) psf_cache=[[psf_cache],[psf]] ; same # of columns as xy_cache xy_cache=[xy_cache,xy] ; This and previous must match up ;tvscl,reform(psf,NX,NX),0 Ncache=Ncache+1 endelse ; 4. Subtract PSF*frac*map_max from map, increment iteration count and map ;print,xy,where(psf EQ max(psf)) ;tvscl,reform(map,NX,NX),1 map=map-psf*(frac*map_max) clean_map(xy)=clean_map(xy)+frac*map_max iter=iter+1 ncomp=ncomp+1 ; 5. If count_rate exists, compute modulation profile and chi_sq if (Keyword_set(count_rate)) then begin plot,count_rate,psym=10,xsty=1 clean_counts = countrate_profile(clean_map,mod_pat) oplot,clean_counts,psym=2 clean_flux=total(clean_counts) chi_sq = chi_square(clean_counts,count_rate) chi_sq_tot=total(chi_sq) endif ; 6. If chi_sq > chi_sq_crit or N < Niter OR clean_flux < obs_flux, go to 2 endif else begin print,"Negative Maximum. Exiting Loop." endelse if ((iter mod 30) EQ 1) then print,label print,iter,Ncache,xy(0),map_max,total(chi_sq),total(clean_counts) endwhile ; 7. Return the result if (n_elements(count_rate) GT 0) then begin count_rate=(clean_counts) print,'Total chi square=',total(chi_sq) print,'Reduced chi square=',total(chi_sq)/n_elements(chi_sq) endif ca=fltarr(ncache) ; compress the clean amplitude table: for i=0,ncache-1 do $ ca(i)=total(clean_map(xy_cache(i+1))) cc=xy_cache(1:*) output=create_struct('clean_map',clean_map,'cc', $ cc,'ca',ca, 'chi_sq',chi_sq,'clean_counts',clean_counts, $ 'params',params) return, output end