SPECTRAL ANALYSIS OF PILED-UP SOURCES USING ACIS EXTRACT $Id: pileup.txt 3454 2009-06-19 17:27:40Z patb $ Patrick Broos, Konstantin Getman, Masahiro Tsujimoto Spring 2004 ===================================================================== INTRODUCTION This recipe performs a series of annular extractions of a single source for the purpose of evaluating how spectral modeling of the source is affected by pileup. This pileup analysis recipe assumes you have already run a nominal AE extraction on a source, and that you are using the directory structure and local ardlib.par file that are constructed in the standard AE recipe "recipe.txt" that is distributed with the AE package. We start from the main extraction directory, where AE was first run. The convention below is that shell commands are not indented, IDL commands are indented 2 spaces. ===================================================================== IDENTIFY THE PILED-UP SOURCE In this example, the piled-up source is called 043737.46-022928.3. Run the following 2 commands in two shells, one in which IDL will be run and one for all the shell commands in the recipe below. setenv SRCNAME "043737.46-022928.3" setenv OBS 5185 ===================================================================== SET UP IDL ENVIRONMENT idl | tee -a pileup.${OBS}.log ; Tell IDL whether pipeline event position randomization was used -- we ; almost always turn this off (=0) unless the observation is short. pipeline_randomization=0 OR pipeline_randomization=1 obs = getenv('OBS') dir = '../obs' + obs + '/' ardlib = dir + 'ardlib.par' src_name = getenv('SRCNAME') catfile = src_name + '.cat' help, obs, dir, catfile ===================================================================== BUILD A BETTER PSF USING ChaRT For pileup analysis we recommend that you use the best set of PSFs you can find. PSFs which AE builds from the PSF Library can differ significantly from ChaRT PSFs, even on-axis (see the "Caveats" section of the AE manual). Running ChaRT is not easy, but we've written a tool (below) that calculates the ChaRT parameters you need, downloads the rays, and processes them into an AE-compatible set of PSF images. idl .run acis_extract_tools ae_chart_interface, src_name, obs Then, cd to the source's observation directory and create a symbolic link that will cause the ChaRT PSFs to be used by AE: mv source.psf AE_source.psf ln -s chart_source.psf source.psf. ===================================================================== CHOOSE A SET OF ANNULI TO EXTRACT ; Guess at a set of core exclusion sizes (expressed in PSF fractions) that ; would nicely sample the range around the "best" sized excluded core. ; Making this guess is very much a recipe area needing improvement. ; It may be that an iterative approach is needed -- make a coarse pass, ; then go back and do finer sampling in the transition region. ; You might coarsely sample the PSF fraction range. inner_frac = (1+indgen(9)) / 10. OR ; You might be more clever by estimating a good exclusion radius ; (e.g. Tsujimoto http://agyo.rikkyo.ac.jp/~tsujimot/arfcorr.html), ; and then figuring out the PSF fraction that ; corresponds to that radius (e.g. by applying a circular region to the ; PSF image). None of this is automated at present ... ; Choose an outer PSF fraction. Consider carefully how much contamination from neighboring ; sources you are willing to tolerate. outer_frac = 0.99 all_frac = [outer_frac,inner_frac] excluded_pcnt = [ 0,inner_frac] * 100 num_frac = n_elements(all_frac) all_name = "annulus" + string(excluded_pcnt*10, F='(%"%3.3d")') + ":" + string(replicate(outer_frac*1000,num_frac), F='(%"%3.3d")') outer_frac_name = all_name[0] forprint, all_name ; These names are annulus???:???, where the 3 spaces gives room for the PSF fraction ; corresponding to the annulus inner radius (core excluded) and outer radius, to one decimal place. save, FILENAME=src_name+'/pileup.sav' ===================================================================== BUILD POLYGONS FOR A VARIETY OF PSF FRACTIONS ; Build catalog with multiple instances of the source, each with a ; different PSF fraction. ; We set RA=0, DEC=0 so that AE will used the stored coordinates. ra = replicate(0, num_frac) dec = replicate(0, num_frac) sourcename = replicate(src_name, num_frac) psf_energ = replicate(1.4967, num_frac) ; Write sources to catalog. fmt = '(A,1x,F10.6,1x,F10.6,1x,F5.3,1x,F7.5)' forprint, TEXTOUT=catfile, sourcename, ra, dec, all_frac, psf_energ, F=fmt, /NoCOMMENT ; Make extraction directories & regions acis_extract, EXTRACTION_NAME=all_name, catfile, obs, dir+'spectral.evt', /CONSTRUCT_REGIONS, EMAP_FILENAME=dir+'obs.emap', ASPECT_FN=dir+'obs.asol', PIPELINE_RANDOMIZATION=pipeline_randomization ===================================================================== MAKE ANNULAR POLYGON SOURCE REGIONS We want to edit the extract.reg files for all annulus>0 extractions so they contain "polygon annuli", e.g. polygon(...) <- the "outer_frac" (e.g. 99%) polygon from annulus000:???/extract.reg -polygon(...) <- the "inner_frac" polygon * The hard, but clear way is to manually edit the extract.reg files: foreach file ($SRCNAME/$OBS/annulus*/extract.reg) open $file end * The fast way is to use some IDL code: .run ; Read ds9 header line and outer polygon from annulus000:???. fl_reg_00 = './'+src_name+'/'+obs+'/'+outer_frac_name+'/extract.reg' readcol, fl_reg_00, reg_00, F='A', delimiter='|' for i=1, n_elements(all_name)-1 do begin ; Read inner polygon and prepend minus sign. fl_reg_curr='./'+src_name+'/'+obs+'/'+all_name[i]+'/extract.reg' print, fl_reg_curr readcol, fl_reg_curr, reg_curr, F='A', delimiter='|', skipline=1 reg_curr[0]='-'+reg_curr[0] ; Rename single-polygon region file, and write polygon annulus. temp=[reg_00[0],reg_00[1],reg_curr] fl_reg_old='./'+src_name+'/'+obs+'/'+all_name[i]+'/extract.reg_old' file_move, fl_reg_curr, fl_reg_old, /OVERWRITE forprint, textout=fl_reg_curr, temp, F='A', /silent, /nocomment endfor end You can check that all the annular region files have two polygons in them via: grep -c polygon $SRCNAME/*/annulus*/extract.reg ===================================================================== EXTRACT SPECTRA and TIMING * Look at regions, and make sure polygons don't cross! (In your non-IDL shell...) grep -v circle $SRCNAME/$OBS/annulus000*/extract.reg > temp.reg foreach reg ($SRCNAME/$OBS/annulus*/extract.reg) grep '-' $reg >> temp.reg end ds9 -log $SRCNAME/$OBS/neighborhood.evt -regionfile temp.reg & * Extract source spectrum. acis_extract, EXTRACTION_NAME=all_name, catfile, obs, dir+'spectral.evt', /EXTRACT_SPECTRA, EMAP_FILENAME=dir+'obs.emap', ASPECT_FN=dir+'obs.asol', ASPHIST_DIR=dir+'asphist', ARDLIB_FILENAME=ardlib, PBKFILE=dir+'obs.pbkfile', MSKFILE=dir+'obs.mskfile' acis_extract, EXTRACTION_NAME=all_name, catfile, obs, dir+'spectral.evt', /TIMING acis_extract, EXTRACTION_NAME=all_name, catfile, obs, /EXTRACT_SPECTRA, /PLOT * In the above plot stage, many of the plots are not interesting, but a few are: The "PSF Fraction" vs. "source #" and "Mean ARF" vs. "source #" plots should correspond to the annuli you defined. The Source Median Energy plot may vary with the size of the excluded core for two reasons. First, the Chandra PSF varies with energy. Second, this Median Energy statistic is contaminated by background and the relative background level varies with the size of the excluded core. In the plot titled "Area of Source Region" the areas "via polyfill" reflect the inner polygon area, not the annular area. ===================================================================== DEFINE A BACKGROUND SPECTRUM We should already have a background.pi file from the regular extraction of the source previously run. AE will find and use this background. NOTE HOWEVER, that backgrounds built by the tool ae_better_backgrounds are constructed for a specific extraction region, i.e. the nominal region with no excluded core. If you have a significant background component from neighboring cataloged sources, then the optimal background region will be different for each annular extraction region. ===================================================================== MERGE THE OBSERVATIONS. PERFORM SPECTRAL FITTING * Final products. acis_extract, EXTRACTION_NAME=all_name, MERGE_NAME=all_name, catfile, /MERGE_OBSERVATIONS * Fitting. acis_extract, MERGE_NAME=all_name, catfile, /FIT_SPECTRA, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm' exit plot_spectra.pl $SRCNAME/annulus*/spectral_models/nogrp_tbabs_vapec/ldata.ps ===================================================================== PLOT FIT RESULTS VS PSF FRACTION idl src_name = getenv('SRCNAME') restore, src_name+'/pileup.sav', /VERBOSE ; Use COLLATED_FILENAME stage to collect fit results. reportfile = 'annulus.collated' acis_extract, MERGE_NAME=all_name, catfile, COLLATED_FILENAME=reportfile bt = mrdfits(reportfile,1) ; Get some photometry values for band 0 (assumed to be 0.5-8 keV). ; Confirm we're using the "full" band: band_full = 0 print, 'Using the energy band ', bt[0].ENERG_LO[band_full], bt[0].ENERG_HI[band_full] NET_CNTS = bt.NET_CNTS[0] SRC_CNTS = bt.SRC_CNTS[0] MEAN_ARF = bt.MEAN_ARF[0] PSF_FRACTION= bt.PSF_FRAC inner_radius = bt.SRC_RAD inner_radius[0] = 0 ; Move to source directory. cd, src_name ; Make some plots. xtit1 = 'size of discarded core (% of PSF)' xtit2 = 'approximate radius of discarded core (sky pixels)' ;; Plot narrow band and summed fluxes (from AE) vs excluded core size. function_1d, idf, excluded_pcnt, bt.FLUX2[9], LINE=6, PSYM=1, DATASET='0.5-1 keV' function_1d, idf, excluded_pcnt, bt.FLUX2[10], LINE=6, PSYM=1, DATASET='1-2 keV' function_1d, idf, excluded_pcnt, bt.FLUX2[11], LINE=6, PSYM=1, DATASET='2-4 keV' function_1d, idf, excluded_pcnt, bt.FLUX2[12], LINE=6, PSYM=1, DATASET='4-6 keV' function_1d, idf, excluded_pcnt, bt.FLUX2[13], LINE=6, PSYM=1, DATASET='6-8 keV' function_1d, idf, excluded_pcnt, total(bt.FLUX2[9:13],1), LINE=6, PSYM=6, DATASET='SUM of bands, 0.5-8keV', XTIT=xtit1, YTIT='FLUX2 from AE', TIT=src_name forprint, bt[0].ENERG_LO[9:13], bt[0].ENERG_HI[9:13], bt[0].FLUX2[9:13] ; Plot some XSPEC parameters. ;CHANGE THESE TO BE APPROPRIATE FOR YOUR MODEL: param_names = ['NH' ,'KT' ,'KT2' ,'NORM' ,'NORM2' ,'F0P5_8' ,'F0P5_2' ,'F2_8' ] uerr_names = ['NH_ERRU','KT_ERRU','KT2_ERRU','NORMERRU','NOR2ERRU','F0P5_8U','F0P5_2U','F2_8U'];upper end of errors lerr_names = ['NH_ERRL','KT_ERRL','KT2_ERRU','NORMERRL','NOR2ERRL','F0P5_8L','F0P5_2L','F2_8L'];lower end of errors .run for ii=0,n_elements(param_names)-1 do begin if tag_exist( bt, param_names[ii], INDEX=param_col ) then begin if tag_exist( bt, uerr_names [ii], INDEX=uerr_col ) AND $ tag_exist( bt, lerr_names [ii], INDEX=lerr_col ) then avg_error = (bt.(uerr_col)+bt.(lerr_col))/2. $ else avg_error = 0 id1=0L & id2=0L function_1d, id1, excluded_pcnt, bt.(param_col), ERROR=avg_error, LINE=6, PSYM=1, XTIT=xtit1, YTIT=param_names[ii], TIT=src_name ; function_1d, id2, inner_radius, bt.(param_col), ERROR=avg_error, LINE=6, PSYM=1, XTIT=xtit2, YTIT=param_names[ii], TIT=src_name endif endfor end ;; The following two plots are not going to be exactly flat because PSF_FRACTION is measured at 1.5 keV whereas NET__CNTS and MEAN_ARF are across the full band. id1=0L & id2=0L function_1d, id1, excluded_pcnt, NET_CNTS/PSF_FRACTION, LINE=6, PSYM=1, XTIT=xtit1, YTIT='NET_CNTS/PSF_FRACTION', TIT=src_name function_1d, id2, inner_radius, NET_CNTS/PSF_FRACTION, LINE=6, PSYM=1, XTIT=xtit2, YTIT='NET_CNTS/PSF_FRACTION', TIT=src_name id1=0L & id2=0L function_1d, id1, excluded_pcnt, MEAN_ARF/PSF_FRACTION, LINE=6, PSYM=1, XTIT=xtit1, YTIT='Mean ARF/PSF_FRACTION', TIT=src_name function_1d, id2, inner_radius, MEAN_ARF/PSF_FRACTION, LINE=6, PSYM=1, XTIT=xtit2, YTIT='Mean ARF/PSF_FRACTION', TIT=src_name function_1d, id3, excluded_pcnt, PSF_FRACTION, LINE=6, PSYM=1, XTIT=xtit1, YTIT='PSF_FRACTION', TIT=src_name !!!! To Display errors you must select Edit->Show Errors !!!! If the y-axis label is off the screen (e.g. in flux plots) press the "X-edit" button (lower right) and change Left Margin from 8 to ~12.