SPECTRAL ANALYSIS OF PILED-UP SOURCES USING ACIS EXTRACT
$Id: pileup.txt 4488 2013-10-02 15:31:29Z psb6 $
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 $SRCNAME/pileup.${OBS}.log
  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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
As of Oct 2013, the ae_chart_interface tool has been removed because ChaRT+MARXv5 cannot simulate PSFs for data processed with the EDSER sub-pixel positioning algorithm.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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 ae         
  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*1, F='(%"%2.2d")') + ":" +  string(replicate(outer_frac*100,num_frac), F='(%"%2.2d")')
  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 the nearest percent.  

  save, FILENAME=src_name+'/pileup.sav'
  

=====================================================================
BUILD POLYGONS FOR EACH OF THE PSF FRACTIONS WE ARE USING

  ; 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'


  
=====================================================================
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

You should think carefully about what background region is appropriate for you annuli.  The existing AE extraction will have defined a background region, depicted by background.emap and background.evt, and extracted a background spectrum, background.pi.  By default, all the annular extrations will use this existing background.

If your annuli suffer significant contamination from nearby sources, then you have a very difficult situation.  Each annuli will have different levels of contamination, and thus each requires a different background region chosen (by the ae_better_backgrounds tool) to correctly subtract that contamination.  This recipe does not yet have instructions for how to run the ae_better_backgrounds tool for each of those annuli.

If your piled source is well isolated from significant neighbors, then a single background region is appropriate for all annuli.  However, you should review the background region that AE produced, by looking at background.emap and background.evt in ds9, to see if the piled source is adequately masked.  The ae_better_backgrounds tool uses the PSF model to understand what portions of the field contain significant light from the source we are extracting, and to exclude those regions from the background.  However, the PSF model has a finite spatial footprint, and thus the background region for a bright source can sometimes contain a few percent of the light from that source.  In such cases you may prefer to hand-draw a background region, background.reg, and then run the background stage of AE.  Or you may choose to assume the background is zero by removing all the background files.


=====================================================================
PRODUCE THE DATA PRODUCTS REQUIRED FOR SPECTRAL FITTING AND FIT EACH ANNULAR SPECTRUM

* 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 = ['NH1'     ,'KT1'     ,'KT2'     ,'NORM1'   ,'NORM2'   ,'FH8' ,'FH2' ,'F28' ]
  uerr_names  = ['NH1_ERRU','KT1_ERRU','KT2_ERRU','NOR1ERRU','NOR2ERRU','FH8U','FH2U','F28U'];upper end of errors
  lerr_names  = ['NH1_ERRL','KT1_ERRL','KT2_ERRU','NOR1ERRL','NOR2ERRL','FH8L','FH2L','F28L'];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.




