RECIPE FOR EXTRACTING XMM POINT SOURCES

Patrick Broos and Leisa Townsley

$Id: xmm_recipe.txt 3455 2009-06-19 17:28:36Z patb $


As of April 2008 CIAO and SAS are not compatible in the same shell.
The conflict seems to be in $LD_LIBRARY_PATH: CIAO wants this to be empty and 
SAS needs it to point to some SAS libraries.

Several SAS scripts expect to find Perl at /usr/local/bin/perl, however I've found
some machines where this is not the correct path.  Ask you sysadmin to make a
symbolic link at /usr/local/bin/perl if necessary.

SAS is less supported on OS-X; for example as of April 2008 the task espfilt simply 
does not exist on the Mac distribution.  So use linux when you can!!!!

In this recipe, UNIX commands are usually not indented; IDL commands are.


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
THE ONLY STEPS IN THIS RECIPE THAT ARE KNOWN TO BE SAFE TO RUN SIMULTANEOUSLY ON
MULTIPLE OBSERVATIONS ARE THE ONES INVOLVING THE ACIS EXTRACT PACAKGE (e.g.epic_extract.pro).

THE ESAS TOOLS MAKE USE OF THINGS WHICH STORE THEIR RESULTS IN A COMMON PLACE,
E.G. fkeypar stored results in its parameter file.

THE SAS TOOLS THEMSELVES MAY OR MAY NOT SUPPORT MULTIPLE INSTANCES.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


=====================================================================
SETUP

THIS RECIPE ASSUMES YOU HAVE DOWNLOADED YOUR DATA FROM THE HEASARC ARCHIVE.

Go to your project directory which should contain the XMM archive directories "ODF" and "PPS":

cd /bulk/hiawatha1/projects/HB3/data/0150960401/

AS OF 2008 THE PATH CAN NOT CONTAIN THE "+" CHARACTER BECAUSE SOME SAS TASKS WILL FAIL!!!

AS OF 2008 THE PATH CAN NOT CONTAIN THE "-" CHARACTER IF YOU'RE GOING TO USE THE ESAS PACKAGE!!!
The Perl scripts in ESAS will parse inputs like $SAS_CCF that contain a "-" as 
the start of a tool option, leading to very confusing errors.


Setup a directory in which to run ESAS, unzip the data, and configure SAS and HEASOFT:

screen -S EE
mkdir diffuse_analysis
gunzip ODF/*
setenv SAS_CCF `pwd`/diffuse_analysis/ccf.cif
setenv SAS_ODF `pwd`/ODF
loadsas
heasoft


#############################################################################
L1->L2 EVENT PROCESSING
#############################################################################

By whatever means you prever, prepare Level 2 event data that's suitable for extraction.
In the recipe below I use the mos-filter script from the XMM-ESAS package for MOS data, 
and a modified version of that tool (pn-filter-beta) for PN data.

---------------------------------------------------------------------
# mos-filter runs the SAS tasks cifbuild, odfingest, emchain and 
# espfilt (along with a few others) to filter and clean the events 
# files.  It produces a qdp plot file showing the light curve and
# indicating the accepted time intervals.  This plot should be 
# examined to get an idea of the extent of any contamination.  

cd diffuse_analysis

mos-filter     $SAS_CCF $SAS_ODF 0  |& tee -a mos-filter.log
pn-filter-beta $SAS_CCF $SAS_ODF 0  |& tee -a pn-filter.log

egrep -i "WARNING|ERROR" mos-filter.log
egrep -i "WARNING|ERROR"  pn-filter.log | egrep -v "errorcolumn|ctiCorrectionFailed"



Multiple event lists will typically be produced by mos-filter and by
pn-filter-beta because there are two MOS detectors and because some observations
have multiple "exposures".
See what cleaned event list files are available in your observation:

ls *clean*

If there is only one exposure, you'll probably see these files:

        mos1S001-clean.fits 
        mos2S002-clean.fits. 
        pnS003-clean.fits

However, you might see more than that, e.g.

        mos1S001-clean.fits
        mos2S002-clean.fits
        mos1U002-clean.fits
        mos2U002-clean.fits
        pnS003-clean.fits

Let's list the MOS1, MOS2, and PN exposures in environment variables to make the recipe more generic.

setenv MOS1_LIST `'ls' mos1*-clean.fits |grep -v oot | cut -c 4-8`
setenv MOS2_LIST `'ls' mos2*-clean.fits |grep -v oot | cut -c 4-8`
setenv PN_LIST   `'ls'   pn*-clean.fits |grep -v oot | cut -c 3-6`
setenv MOS_LIST "$MOS1_LIST $MOS2_LIST"

echo $MOS_LIST
echo $PN_LIST

---------------------------------------------------------------------
EACH OF THESE EVENT LISTS IS GOING TO PROCESSED SEPARATELY BY THE ESAS TOOLS.
THIS RECIPE ATTEMPTS TO BE GENERIC, USING $MOS1_LIST, $MOS2_LIST, $PN_LIST WHERE POSSIBLE.
YOU MAY HAVE TO MODIFY YOUR PROCEDURE IN SOME CASES.
---------------------------------------------------------------------


---------------------------------------------------------------------
Examine these plots which depict the background flare filtering that has been applied.

foreach file (*hist.qdp)
  echo $file
  qdp  $file
  end
Use /XW or /CPS or filename/ps when qpd prompts you for a device.
Type "quit" to move to the next plot.


        
---------------------------------------------------------------------
These event lists need a little more filtering to become level 2 data ready for extraction.
The PATTERN and FLAG filtering below follow the recommendations in the document
EPIC status of calibration and data analysis (XMM-SOC-CAL-TN-0018), issue 2.0 found at
http://xmm.vilspa.esa.es/docs/documents/CAL-TN-0018-2-0.pdf

Note however, that epic_extract.pro will later apply its own PATTERN and FLAG filtering.

cd ..

foreach exposure ($MOS_LIST)
  printf "\nMOS $exposure\n"
  mkdir obsmos${exposure}
  evselect table=diffuse_analysis/mos${exposure}-clean.fits:EVENTS filteredset=obsmos${exposure}/spectral.evt expression='(PATTERN<=12)&&(#XMMEA_EM)&&(PI in [200:10000])' filtertype=expression withfilteredset=yes  updateexposure=yes filterexposure=yes verbosity=0
  end
  
  
foreach exposure ($PN_LIST)
  printf "\nPN $exposure\n"
  mkdir obspn${exposure}
  evselect table=diffuse_analysis/pn${exposure}-clean.fits:EVENTS filteredset=obspn${exposure}/spectral.evt expression='(PATTERN<= 4)&&(FLAG == 0)&&(PI in [200:15000])' filtertype=expression withfilteredset=yes updateexposure=yes filterexposure=yes verbosity=0
  end

ls obs*

  
---------------------------------------------------------------------
For display purposes, create a merged MOS event list and a merged PN event list.
It's not easy to script this for an indeterminate set of event files.

'rm' MOS.evt
set expid=`echo $MOS_LIST | cut -d " " -f1`
ln -s obsmos${expid}/spectral.evt MOS.evt
foreach expid (`echo $MOS_LIST | cut -sd " " -f2-`)
  'rm'       temp.evt
  mv MOS.evt temp.evt
  echo "Merging $expid ..."
  merge  set1=temp.evt  set2=obsmos${expid}/spectral.evt  outset=MOS.evt
  end
  

'rm' PN.evt
set expid=`echo $PN_LIST | cut -d " " -f1`
ln -s obspn${expid}/spectral.evt PN.evt
foreach expid (`echo $PN_LIST | cut -sd " " -f2-`)
  'rm'      temp.evt
  mv PN.evt temp.evt
  echo "Merging $expid ..."
  merge  set1=temp.evt  set2=obspn${expid}/spectral.evt  outset=PN.evt
  end
  

---------------------------------------------------------------------
For display purposes, combine the MOS and PN data and apply an standard energy filter.

'rm' temp.evt
merge  set1=MOS.evt  set2=PN.evt outset=EPIC.evt
evselect table=EPIC.evt:EVENTS  filteredset='EPIC_300-8000.evt'  expression='(PI IN [300:8000])' filtertype=expression withfilteredset=yes updateexposure=yes filterexposure=yes verbosity=0

ds9 -bin factor 128 MOS.evt PN.evt EPIC.evt EPIC_300-8000.evt &
 



#############################################################################
SOURCE DETECTION
#############################################################################

The pipeline produces a source list, by unknown methods, possibly on data suffering from background flares.
You can browse pipeline products by opening pps/INDEX.HTM.  Some interesting pages include:

* EPIC Summary : EPIC Summary Source List
* EPIC Summary : EPIC Observation Image
* EPIC Summary : EPIC Source ds9 Regions

The list of pipeline images can be seen by opening pps/*ESKY*

For now, we're going to use the pipeline catalog in this recipe rather than dive into
source detection ourselves.



#############################################################################

#################### AE SET UP #################### 

#############################################################################

=====================================================================
GATHER DATA PRODUCTS NEEDED

All we need are the spectral.evt files in obsmos1S001/ obsmos2S002/ obspnS003/ created above.
For this XMM workflow we don't use exposure maps or PSFs, and whatever auxillary 
observation files are needed by SAS are found via $SAS_CCF and $SAS_ODF.



=====================================================================
GATHER A LIST OF THE CELESTIAL COORDINATES OF PROPOSED SOURCES


mkdir point_sources
cd    point_sources

        
We use the pipeline catalog.
gunzip -c -S FTZ  ../pps/*EPX*OBSMLI*FTZ  > pipeline_cat.fits 
 
idl |& tee -a catalog.log
  bt = mrdfits('pipeline_cat.fits', 1)

  .run ae_recipe
  sourcename = ''   ; Have ae_source_manager return the list of names it creates.
  ae_source_manager, /ADD, RA=bt.RA, DEC=bt.DEC, PROVENANCE='edetect_chain', POSITION_TYPE='emldetect', NAME=sourcename

  .run acis_extract_tools
  srclist_fn = 'temp.srclist'
  forprint, sourcename, TEXTOUT=srclist_fn, /NoCOMMENT, /SILENT
  
  ae_poke_source_property, srclist_fn, KEYWORD='ERR_DATA', VALUE=bt.RADEC_ERR, COMMENT='error circle (arcsec) around (RA,DEC)'

Before proceeding with the extraction you may wish to sort the catalog in
whatever way you need for publication. Unless you're using custom LABELs for
sources, you'll want to relabel them using the final sequence numbering scheme.

  ae_source_manager, /SORT_RA
  ae_source_manager, /SET_LABEL_AS_SEQUENCE
  
  
  
#############################################################################

######### STANDARD EXTRACTION  #################################

#############################################################################

=====================================================================
CONSTRUCT SINGLE-OBSID CATALOGS 

For each obsid and MOS detector construct a catalog that leads to non-overlapping regions.

!!! I DO NOT KNOW if these can be run in parallel on different CPUs. The SAS tools spawned
!!! may very well use temporary files that conflict.

setenv OBS MOS1 ; set prompt="$OBS %% "
setenv OBS MOS2 ; set prompt="$OBS %% "
setenv OBS PN   ; set prompt="$OBS %% "

nice idl |& tee -a construct_regions.${OBS}.log

  obsname = getenv('OBS')
  dir     = '../obs' + obsname + '/'
  epic_extract, 'all.srclist', obsname, dir+'spectral.evt', /CONSTRUCT_REGIONS

egrep -i "WARNING|ERROR|halted" construct_regions.*.log | egrep -v "DISPLAY variable|srcThetaTooLarge|No HDUNAME"

grep -h physical */$OBS/extract.reg    > ../obs${OBS}/extract.reg
grep -h physical */$OBS/background.reg > ../obs${OBS}/background.reg

ds9 -bin factor 128 ../obs${OBS}/spectral.evt -region ../obs${OBS}/extract.reg ../obs${OBS}/spectral.evt -region ../obs${OBS}/background.reg &


=====================================================================
EXTRACT SPECTRA AND MAKE RESPONSE FILES

!!! I DO NOT KNOW if these can be run in parallel on different CPUs. The SAS tools spawned
!!! may very well use temporary files that conflict.

setenv OBS MOS1 ; set prompt="$OBS %% "
setenv OBS MOS2 ; set prompt="$OBS %% "
setenv OBS PN   ; set prompt="$OBS %% "

Circa April 2008 you need to run this in a shell that does not have  SAS configured, 
but you do need the SAS environment variables set.

nice idl |& tee -a extract_spectra.${OBS}.log
 
  obsname = getenv('OBS')
  dir     = '../obs' + obsname + '/'
  epic_extract, 'all.srclist', obsname, dir+'spectral.evt', /EXTRACT_SPECTRA

   
egrep -i "WARNING|ERROR|halted" extract_spectra.*.log | egrep -v "DISPLAY variable|EnergyOutsideValidityRange|NoFilteredEvents"




=====================================================================
DECIDE WHICH SOURCES ARE INVALID & DISCARD THEM
=====================================================================

See the catalog pruning section in our ACIS recipe (recipe.txt) for a procedure to compute source validity statistics on three bands, and then prune the catalog by thresholding those statistics.

Note that the current ACIS recipe allows for only three energy bands, but Antokhin (A&A 477, 593-609 2008) recommends accepting a source if "detected" in any of four bands:
  full (0.4 - 10.0 keV)
  soft (0.4 -  1.0 keV)
  med  (1.0 -  2.5 keV)
  hard (2.5 - 10.0 keV)
We could make the code smarter if 4 bands are really desirable...




#############################################################################

#################### SECTION V: SPECTRAL FITTING #################### 

An excellent introduction to the concept of fitting X-ray spectra can be found in
Keith Arnaud's PowerPoint presentation at 
http://lheawww.gsfc.nasa.gov/~kaa/india2003/lowresspec.ppt

#############################################################################

=====================================================================
COLLATE FINAL CATALOG

  acis_extract, 'all.srclist',  COLLATED_FILENAME='tables/all.collated', REGION_FILE='all.reg', LABEL_FILENAME='label.txt'

egrep -i "WARNING|ERROR|halted" merge.log | egrep -v "no in-band|different value|not merged"

ds9 -bin factor 128 ../EPIC_300-8000.evt -region all.reg &


=====================================================================
FIT A VARIETY OF SPECTRAL MODELS TO EACH SOURCE 

Carry on with the fitting recipe in recipe.txt.





#############################################################################

#################### DIFFUSE SOURCES #################### 

#############################################################################
This recipe is derived from the example in the ESAS manual.

cd diffuse_analysis

=====================================================================
ADDITIONAL CLEANING OF EVENT DATA

EVENT LISTS
As far as I can tell the event list coming out of mos-spectra (mos1S001-clean.fits and
mos2S002-clean.fits) are the most-reduced event lists available.  
Further on-the-fly filtering is done when images are made, but I don't see any event list
derived from these.

We did not place this cleaning in the workflow earlier because we want the point source analysis
to use data from CCDs that were in their "high background mode" (high count rates at energies < 1 keV).


# At this point check for any CCD's running
# in their high background modes.  Occasionally certain CCDs act up
# and have high count rates at energies < 1 keV.  The background modeling
# software doesn't work very well on these so for the analysis of diffuse
# data they should be excluded.  Run SAS and xmmselect and create an image
# using a selection expression of: 
#      (FLAG == 0)&&(PATTERN <= 12)&&(PI in [200:1000])
# If the CCD is acting up it will be very obvious.  You can determine 
# which CCD by looking at the Users Handbook, or simply using the process 
# of elinimation.  Add &&!(CCDNR == 1), etc., until the correct CCD is
# removed from the image.  For this observation MOS2 CCD #5 is bad.  
# Also, MOS1 CCD #6 is dead since the observation was done after the 
# micrometeorite


idl
  evtfile = file_search('*clean*', COUNT=count)
  .run
  for ii=0,count-1 do begin
    fn = evtfile[ii]
    cmd = string(fn, F="(%'ciao; dmcopy ""%s[flag=0,pattern<=12,pi=200:1000]"" temp.evt clob+')")
    print, cmd
    spawn, cmd
    
    bt=mrdfits('temp.evt',1)
    hist = histogram(bt.CCDNR,MIN=1,MAX=7)
    forprint, 1+indgen(7), hist/float(median(hist)), F='(%"CCD: %d  rate: %4.2f")'
 
    spawn, 'ciao; dmcopy "temp.evt[bin detx=::200,dety=::200]" temp.img clob+'
    spawn, 'ds9 temp.img &'

    s=''
    read,fn+': ', s
  endfor
  end

The code above will display a soft image for each detector and will show a table of relative count rates
to help you judge which CCDs need to be excluded from diffuse analysis.

In this observation we need to remove CCD5 from the "clean" event list for MOS2.
(We cannot use dmcopy for this job because it seems to mess up GTI tables which makes cheese have problems.)

mv mos2S002-clean.fits temp.fits
evselect table='temp.fits:EVENTS' expression='CCDNR != 5' filteredset=mos2S002-clean.fits 

ds9 -bin factor 128 *clean.fits &


=====================================================================
SOURCE DETECTION

# Next do the source detection and create the point-source masks.  
# cheese-tight does source detection on the observation images and creates
# swiss cheese masks which can be used later.

cheese `echo $MOS_LIST|wc -w` $MOS_LIST  2 0 0.9 1.0 0.0 1 300 8000 |& tee -a cheese.log  

egrep -i --text "WARNING|ERROR" cheese.log

The images used for source detection (mos1S001-obj-im.fits and mos2S002-obj-im.fits)
are unfortunately removed by the tool.  To re-create them run:

evselect table=mos1S001-clean.fits:EVENTS imageset=mos1S001_300-8000.img withfilteredset=yes expression='(PATTERN<=12)&&(FLAG == 0)&&(PI in [300:8000])' filtertype=expression imagebinning='imageSize' imagedatatype='Int32' squarepixels=yes ignorelegallimits=yes withxranges=yes withyranges=yes xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 updateexposure=yes filterexposure=yes verbosity=0
                                                             
evselect table=mos2S002-clean.fits:EVENTS imageset=mos2S002_300-8000.img withfilteredset=yes expression='(PATTERN<=12)&&(FLAG == 0)&&(PI in [300:8000])' filtertype=expression imagebinning='imageSize' imagedatatype='Int32' squarepixels=yes ignorelegallimits=yes withxranges=yes withyranges=yes xcolumn='X' ximagesize=900 ximagemax=48400 ximagemin=3401 ycolumn='Y' yimagesize=900 yimagemax=48400 yimagemin=3401 updateexposure=yes filterexposure=yes verbosity=0


The catalog from cheese, without duplicates, can be obtained by selecting "summary band" entries:

dmcopy "../diffuse_analysis/emllist.fits[ID_BAND=0]" catalog.fits


Circa 2008 April, the tool "cheese" in ESAS does single-band detection (300-8000 eV).
This could be improved to do multi-band detection, or replaced with some other detection script.

It's not clear whether the detection tools inside cheese (eboxdetect and
emldetect) effectively combine the images that are supplied (i.e. combine the
MOS1 and MOS2 data in our case). Section 4.12.3 of the SAS manual clearly shows
passing these tools five narrow band images, and NOT passing a full-band image.
Leigh Jenkins' detection scripts do the same thing. This practice suggests that
each image is searched individually and they are somehow searched in
combination, although the "algorithm" section of the manual for eboxdetect does
not show how this would occur.




=====================================================================
EXTRACTION and BACKGROUND SPECTRA

# mos-spectra processes the filtered event files to produce background 
# spectra for the entire energy range and selected region and background 
# images for the selected region and selected band for the individual CCDs.  
# The region selection expression is in an input file and should be in 
# detector coordinates.  If the input file does not exist (as in this case
# since you want the entire FOV), nofile.reg in this case, the default is 
# to process the entire FOV.  The input energies are in eV.  

# xmm-back takes the individual ccd spectra and images from the mos-spectra
# processing and combines them into complete spectra and images for the 
# selected region.  The resultant image is in detector coordinates.  xmm-back 
# creates a QDP plot file which shows the source and model background spectra 
# for the observation.  Any discrepancies at higher energies probably indicate 
# residual soft proton contamination, unless there are really hard and bright 
# sources in the field.  The QDP files have names like: mos1S001-spec.qdp


In this recipe we omit CCD5 on MOS2 (via the 0 in the 5th flag input to mos-spectra and xmm-back below) 
to be safe, but you could leave it in if that CCD is behaving in your observation.

foreach exposure ($MOS1_LIST)
  printf "\nMOS $exposure\n"
  mos-spectra $exposure $XMM_ESAS/caldb/ 0 nofile.reg 1 500 4500 1 1 1 1 1 0 1  |& tee -a mos-spectra.log
  xmm-back    $exposure $XMM_ESAS/caldb/                500 4500 1 1 1 1 1 0 1  |& tee -a xmm-back.log
  end

foreach exposure ($MOS2_LIST)
  printf "\nMOS $exposure\n"
  mos-spectra $exposure $XMM_ESAS/caldb/ 0 nofile.reg 1 500 4500 1 1 1 1 0 1 1  |& tee -a mos-spectra.log
  xmm-back    $exposure $XMM_ESAS/caldb/                500 4500 1 1 1 1 0 1 1  |& tee -a xmm-back.log
  end
 
egrep -i --text "WARNING|ERROR" mos-spectra.log | egrep -v "NoExpoExt|exposure duration|zero normalisation|range of image"
egrep -i --text "WARNING|ERROR" xmm-back.log    | egrep -v "Hard error|Rate error"


foreach file (*spec.qdp)
  qdp $file
  end
Use /XW or /CPS or filename/ps when qpd prompts you for a device.


=====================================================================
SMOOTHED IMAGE

# rot_back uses information in a previously created count image
# in sky coordinates to rotate the detector coordinate background
# images into images in sky coordinates.

foreach exposure ($MOS1_LIST)
  printf "\nMOS $exposure\n"
  rot-im-det-sky $exposure 0 500 4500 1  |& tee -a rot-im-det-sky.log
  end
  
foreach exposure ($MOS2_LIST)
  printf "\nMOS $exposure\n"
  rot-im-det-sky $exposure 0 500 4500 1  |& tee -a rot-im-det-sky.log
  end
  
egrep -i --text "WARNING|ERROR" rot-im-det-sky.log

# comb combines the MOS1 and MOS2 images, as well as images from multiple 
# exposures.  Note that there is a hook for combining soft proton background 
# images.  Creating a soft proton image requires spectral fitting first.


comb 1 500 4500 1 `echo $MOS_LIST|wc -w` $MOS_LIST  |& tee -a comb.log

egrep -i --text "WARNING|ERROR" comb.log


# Run adapt-900 to smooth the particle-background-subtracted 
# and exposure-corrected image

adapt-900 50 0.02 0 2 500 4500 1 0  |& tee -a adapt-900.log

egrep -i --text "WARNING|ERROR" adapt-900.log


Examine the combined data, emap, and background images:

ds9 mos-obj-im-500-4500.fits mos-exp-im-500-4500.fits mos-back-im-sky-500-4500.fits &

You should be able to Frame:Match Frames:WCS and see that the three are aligned in celestial coordinates.


Rename the smoothed image.  Suffix "nsp" means "no soft proton".

mv adapt-500-4500.fits adapt-500-4500-nsp.fits


Display the smoothed image with a point source catalog overlaid.
Here we use the catalog from the point source extraction, rather than the one computed by "cheese" above.
Note that "cheese" does not mask every source it detects.
 
srcdisplay imageset=adapt-500-4500-nsp.fits boxlistset=emllist.fits sourceradius=0.01 withregionfile=true regionfile=emllist.reg 

ds9 mos-obj-im-500-4500.fits -region emllist.reg -region ../point_sources/all.reg adapt-500-4500-nsp.fits -region emllist.reg -region ../point_sources/all.reg &


=====================================================================
SPECTRAL ANALYSIS

See example script distributed with ESAS if you want to carry on with fitting diffuse emission.