; name : xpsf.pro ; author : george chartas ; date : 5 July 2000 ; purpose : ; arf files generated by mkarf ; are modified to account for the dependence of ; photon scattering with energy. This effect is ; important whenever extraction regions are used ; in spectral analysis. ; One application of this tool is performing ; spectral fits to piled-up spectra. ; A spectrum extracted from an annulus centered ; on the piled-up source will be less affected by pile-up. ; XPSF simulates the effective area (EA) of CHANDRA ; for a specified location and extraction region ; and source morphology using XSPEC and MARX. ; ; The ratio f(E) = EA(ROI) / EA(2pi) is the ; arf correction function ; ; XPSF creates a new arf file with a corrected effective ; area column xtrad_in = par_get("xpsf.par","xtrad_in",status1) xtrad_out = par_get("xpsf.par","xtrad_out",status1) psfsize = par_get("xpsf.par","psfsize",status1) OutputDir = (par_get("xpsf.par","OutputDir",status1)) DetectorType = (par_get("xpsf.par","DetectorType",status1)) ccdid = fix((par_get("xpsf.par","ccdid",status1))) DetIdeal = (par_get("xpsf.par","DetIdeal",status1)) DetOffsetX = par_get("xpsf.par","DetOffsetX",status1) DetOffsetY = par_get("xpsf.par","DetOffsetY",status1) DetOffsetZ = par_get("xpsf.par","DetOffsetZ",status1) SourceOffsetZ = par_get("xpsf.par","SourceOffsetZ",status1) SourceOffsetY = par_get("xpsf.par","SourceOffsetY",status1) NumRays = (par_get("xpsf.par","NumRays",status1)) dNumRays = (par_get("xpsf.par","dNumRays",status1)) SourceFlux = fix(par_get("xpsf.par","SourceFlux",status1)) SpectrumType = (par_get("xpsf.par","SpectrumType",status1)) SpectrumFile = (par_get("xpsf.par","SpectrumFile",status1)) SourceType = (par_get("xpsf.par","SourceType",status1)) SGaussSigma = (par_get("xpsf.par","S-GaussSigma",status1)) SBetaCoreRadius = (par_get("xpsf.par","S-BetaCoreRadius",status1)) SBetaBeta = (par_get("xpsf.par","S-BetaBeta",status1)) xspecmodel = (par_get("xpsf.par","xspecmodel",status1)) xspecparam = (par_get("xpsf.par","xspecparam",status1)) xspec_model_file = (par_get("xpsf.par","xspec_model_file",status1)) GratingType = (par_get("xpsf.par","GratingType",status1)) arf = (par_get("xpsf.par","arf",status1)) new_arf = (par_get("xpsf.par","new_arf",status1)) Dithermodel = "NONE" ; Simulate psf 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, " dNumRays=",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")' 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 printf,lun,model_pars(i) printf,lun,'dummyrsp 0.05 10 2048' 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 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 = strtrim(NumRays) OutputDir = strtrim(OutputDir) GratingType = strtrim(GratingType) DetectorType = strtrim(DetectorType) DetIdeal = strtrim(DetIdeal) SpectrumType = strtrim(SpectrumType) SpectrumFile = strtrim(SpectrumFile) SourceType = strtrim(SourceType) DitherModel = strtrim(DitherModel) cmd = string(NumRays,dNumRays,OutputDir,GratingType,DetectorType, $ DetIdeal,SourceFlux,SpectrumType,SpectrumFile,SourceType,SourceOffsetZ, $ SourceOffsetY,SGaussSigma,SBetaCoreRadius,SBetaBeta,DitherModel, FORMAT=f ) spawn, 'echo "' + cmd + '" >> marxlog.txt ' spawn, "source /bulk/pkg/asc/ascds.csh ;" + cmd,result ; create psf ypsf = read_marx_file(strtrim(Outputdir) + '/xpixel.dat') zpsf = read_marx_file(strtrim(Outputdir) + '/ypixel.dat') ene = read_marx_file(strtrim(Outputdir) + '/energy.dat') pha = read_marx_file(strtrim(Outputdir) + '/pha.dat') ; find centroid of psf ; ypsf_m = mean(ypsf) ; zpsf_m = mean(zpsf) if (DetectorType eq 'ACIS-S' and ccdid eq 7) then begin ypsf_m1 = 263.5 - (SourceOffsetY/(0.492/60.)) zpsf_m1 = 507.4 - (SourceOffsetZ/(0.492/60.)) endif if (DetectorType eq 'ACIS-I' and ccdid eq 3) then begin ypsf_m1 = 959.4 - (SourceOffsetZ/(0.492/60.)) zpsf_m1 = 951.2 + (SourceOffsetY/(0.492/60.)) endif if (DetectorType eq 'ACIS-I' and ccdid eq 1) then begin ypsf_m1 = - (1024.- 959.4 + 21.7) - (SourceOffsetZ/(0.492/60.)) zpsf_m1 = 951.2 + (SourceOffsetY/(0.492/60.)) endif if (DetectorType eq 'ACIS-I' and ccdid eq 0) then begin ypsf_m1 = 2048. - (-(SourceOffsetZ/(0.492/60.)) - (- 959.4 + 21.7)) zpsf_m1 = 2048. - (951.2 - 20.5 + (SourceOffsetY/(0.492/60.))) endif if (DetectorType eq 'ACIS-I' and ccdid eq 2) then begin ypsf_m1 = 1024. - (959.4 - (SourceOffsetZ/(0.492/60.))) zpsf_m1 = 2048. - 951.2 + 20.5 - (SourceOffsetY/(0.492/60.)) endif ; plot,ypsf,zpsf,psym=3 ; stop sub = 0.5 ; 0.125 ; 0.25 psf = make_image(ypsf,zpsf,xbinsize=sub*1.01626,ybinsize=sub*1.01626,xrange=[ypsf_m1-psfsize,ypsf_m1+psfsize],yrange=[zpsf_m1-psfsize,zpsf_m1+psfsize]) ; bin size of psf is sub*1.01626*0.492 arcsec psf_max = max(psf,ind) psf_index = ind ymax= ind mod n_elements(psf(*,0)) zmax= ind/n_elements(psf(0,*)) ; determine the total number of events in the entire psf ntot = n_elements(ypsf) ; determine the total number of events within an annulus of inner radius xtrad_in ; and outer radius xtrad_out centered on the psf rad_evt = sqrt((ypsf - ypsf_m1)^2. + (zpsf - zpsf_m1)^2.) index = where(rad_evt le xtrad_out and rad_evt ge xtrad_in) nxtr = n_elements(index) psf_cor = (1.d0*nxtr)/(1.d0*ntot) ; determine the spectrum of all events in the simulated psf; spec1 binsize = 0.1 ; binsize in keV spec1 = histogram(ene,min=0.1,max=10,bin=binsize) ; determine the spectrum of events within a circle of radius xtrad ; centered on the psf; spec2 spec2 = histogram(ene(index),min=0.1,max=10,bin=binsize) energy = findgen(n_elements(spec1))*0.1 erase !noeras=-1 ymax = max(spec1) ymin = min(spec1) !p.charsize = 1.5 set_viewport,0.15,0.9,0.1,0.4 !xtitle='Energy (keV)' !ytitle='Counts per bin' plot_io,energy,spec1,psym=10,yrange = [ymin,ymax] plot_io,energy,spec2/psf_cor,psym=10,yrange = [ymin,ymax] set_viewport,0.15,0.9,0.5,0.9 !xtitle='' !ytitle='Arf Correction Factor' arf_cor = (spec2*1.)/(spec1*1.) plot,energy,arf_cor,psym=10 ; calculate the modified effective area values at the same energies ; as the user supplied arf orig_arf = mrdfits(arf,1) arf_energ_lo = orig_arf.energ_lo arf_energ_hi = orig_arf.energ_hi arf_specresp = orig_arf.specresp arf_energ_mid = (arf_energ_lo + arf_energ_hi)/2. arf_cor_sm = smooth(arf_cor,3) arf_cor_int = interpol(arf_cor_sm,energy,arf_energ_mid) new_arf_specresp = arf_cor_int*arf_specresp nelem = n_elements(arf_specresp) ; correct .arf file ; read in original arf file header orig_arf_header0 = HEADFITS(arf,exten=0) orig_arf_header1 = HEADFITS(arf,exten=1) h0 = orig_arf_header0 h1 = orig_arf_header1 writefits,new_arf,0,h0 structure = replicate({ENERG_LO:0.0,ENERG_HI:0.0,SPECRESP:0.0},nelem) structure.ENERG_LO = arf_energ_lo structure.ENERG_HI = arf_energ_hi structure.SPECRESP = new_arf_specresp mwrfits,structure,new_arf,h1 end