RECIPE FOR ACIS EXTRACT USAGE FOR A FIELD WITH MULTIPLE MISALIGNED OBSERVATIONS $Revision: 4207 $ $Date: 2012-04-05 15:01:45 -0400 (Thu, 05 Apr 2012) $ Patrick Broos patb@astro.psu.edu with assistance from Leisa Townsley townsley@astro.psu.edu ############################################################################# PLEASE SEND ME FEEDBACK ON THIS RECIPE IF YOU FIND MISTAKES OR THINK OF IMPROVEMENTS!!! I SPEND MORE TIME WRITING CODE AND DOCUMENTATION THAN ACTUALLY EXTRACTING ACIS DATA. THIS RECIPE WILL IMPROVE MUCH MORE RAPIDLY IF EVERYONE USING IT REPORTS PROBLEMS. Patrick Broos ############################################################################# ############################################################################# ####################### NOTES ABOUT THIS RECIPE ############################# ############################################################################# The "master" list of sources is the ASCII file 'all.srclist' maintained by ae_source_manager; the coordinates of those sources are stored in the source.stats FITS files in each source directory. The sources observed by each ObsID are listed in the 5-column "catalog" named obsXXXX/obs.cat which is constructed by the tool ae_make_catalog; these single-ObsID catalogs are needed because we're adjusting PSF fractions for crowding which can differ between ObsIDs (due to different off-axis angles). The "tee" syntax used in this recipe to start IDL causes CIAO and IDL output to be printed both to the screen and to the log file specified. Note the use of "|&" to redirect stdout and stderr. This "tee" mechanism is HIGHLY RECOMMENDED for all AE calls since it forms a very important record of your work. Most CIAO and HEASOFT commands run by AE are echoed to this log file, giving the observer at least a partial record of what was done to their data. More importantly, error messages produced by CIAO/HEASOFT tools and by AE itself may be in the log file. It's not obvious how to find those messages however. In this recipe these logs are searched using "egrep" for the words WARNING or ERROR. Various undocumented AE options are used in this recipe to speed up processing of large catalogs by skipping various computations when they are not needed. The recipe refers to a number of other data analysis recipes that we have developed for use in our own work, generally noted as "(not released to the public)". In this recipe, UNIX commands are indented two spaces; IDL commands are indented four spaces. ############################################################################# #################### SECTION I: SET UP #################### ############################################################################# ============================================================================= Record in your notes the version of this procedure that you are executing, which you can find at the top of this file. ============================================================================= ============================================================================= CONFIGURE YOUR COMPUTING ENVIRONMENT ============================================================================= Once and forever, edit your .cshrc to add the following lines so that $PROMPT, if available, will set your unix prompt to show the ObsID number being processed in each 'screen' window: if ($?PROMPT == 1) then set prompt="$PROMPT" endif ===================================================================== GATHER DATA PRODUCTS NEEDED BY AE --------------------------------------------------------------------- EVENT DATA The strategy we use in OUR OWN work is to perform a partial preliminary AE run for the purposes of (1) evaluating the validity of the proposed sources, via the PROB_NO_SOURCE statistic computed by AE, and (2) improving the positions of those sources, via AE's position estimates. The event list produced by our own L1->L2 processing (referenced by the symlink "spectral.evt") has undergone two different styles of "cleaning": 1. Events satisfying cleaning criteria that are appropriate for ALL point sources have been removed. 2. Events satisfying cleaning criteria that are appropriate for only FAINT point sources have been flagged by setting a bit in their STATUS word. The two criteria currently in this category are: * Events attributed to "flaring pixels", identified by the aggressive acis_detect_afterglow algorithm * Events identified by the Clean55 algorithm for data taken in Very Faint mode (5x5 pixel event islands), which are attributed to instrumental background. AE will apply a STATUS=0 filter to all point source extractions that it judges would benefit from the aggressive cleaning criteria (#2 above), namely those with a low count rate. No STATUS filter is applied by AE for high-rate sources because such cleaning erroneously removes events from the source. --------------------------------------------------------------------- EXPOSURE DATA For each ObsID, AE requires an exposure map that corresponds to the event data you are extracting. For example, if your event list includes only data from the I-array CCDs, then your exposure map must include only those same CCDs. If you wish to hide a portion of the field from AE then you should apply a mask to the exposure map, and then remove from the event list events where the exposure map is zero. AE requires a directory containing aspect histogram files for each CCD used in the observation, named using the convention ccd0.asphist, ccd1.asphist, ... See the "Getting Started" and "EXTRACT_SPECTRA Stage" sections of the AE manual. Aspect histograms are required by the CIAO tool mkarf, and are a by-product of constructing an exposure map. Below, a symbolic link is made to the "asphist" directory left over from using the AE tool ae_make_emap to construct an exposure map; if you build the exposure map another way then simply make a directory named "asphist" and copy the aspect histograms there, using the specified naming convention. --------------------------------------------------------------------- CALIBRATION INFORMATION The ARF-generating commands in CIAO require an "ardlib.par" parameter file that has been configured to contain observation specific information, as discussed in the CIAO threads. Thus, AE requires that you provide such an ardlib.par file for each observation you are extracting. You should have these observation-specific arlib.par files on-hand already since they were required for exposure map construction. See the CIAO threads, and the "Getting Started" and "EXTRACT_SPECTRA Stage" sections of the AE manual for information on configuring ardlib.par. --------------------------------------------------------------------- STANDARD PATHS USED IN RECIPE This AE recipe requires specific standard paths to the input data products (event list, exposure map, aspect file, aspect histogram directory, ardlib.par, PBK file, and MSK file) for each observation. Below we call a UNIX shell script supplied with AE that constructs the standard directory structure required by this recipe, using symbolic links ("ln -s") with standard names to point to your actual data products (which can have any names). The script presumes that each ObsID's L2 files are found in a directory with a name of the form obsid_* , and presumes that the L2 files in that directory follow the naming conventions we adopt in our own L1-->L2 reduction process (not the naming conventions of the CXC's L2 files). You can edit the script as required to match the naming conventions you use. Our own projects usually use observations taken with the I-array aimpoint, with 6 CCDs enabled. We often choose to extract only the I-array CCDs, and our L1-->L2 reduction process generates event and emap products for the I-array. The AE setup script setup_AE_for_iarray_data.csh will configure such an extraction. Move to whatever directory you wish to contain the AE extractions, and pass the desired AE setup script a shell "glob pattern" that finds the directories containing the L2 products for your target. cd /data/extract/ /setup_AE_for_iarray_data.csh "../pointing_*/obsid_*" ===================================================================== GATHER A LIST OF THE CELESTIAL COORDINATES OF PROPOSED SOURCES Create a directory named "point_sources.noindex" in which AE will be executed. (The ".noindex" suffix in the name is useful on Mac machines to tell Spotlight not to waste time indexing the contents of the directory.) This should be a SIBLING to the obsXXXX directories created above, i.e. the directory tree will look like this: /data/extract/ target.evt (a symlink to an event list suitable for visual review of regions) obs1874/ obs1875/ obs1876/ obs1877/ obs3750/ point_sources.noindex/ ... cd /data/extract/ ln -s ../iarray.target.evt target.evt mkdir point_sources.noindex cd point_sources.noindex mkdir tables ln -s Pb_full_band.collated tables/Pb_template Throughout the recipe we'll be using the tool ae_source_manager to add, remove, and reposition sources. This tool maintains a list of all the sources in the file "all.srclist" so that it matches the set of source extraction directories. Gather celestial coordinates (in decimal degrees) for the most complete set of proposed sources you can find, by whatever means you choose (e.g. wavdetect, human examination of the data, reconstructions, catalogs from other observatories, etc.). Read them into IDL as DOUBLE PRECISION vectors using appropriate tools. Using the AE tool ae_source_manager, create source names from the coordinates, store source information in AE source directories, and maintain a list of the sources in "all.srclist". The optional PROVENANCE and POSITION_TYPE inputs to ae_source_manager are used to record how the source was detected, and what method has been used to estimate the position. PROVENANCE values are stored as FITS keywords (via fxaddpar.pro), and thus are restricted only by the FITS standard (e.g. strings have to be less than 66 characters, I think). The values are for use by the observer, not by AE itself, so use whatever naming conventions are convenient. The optional LABEL input to ae_source_manager defines the source labels that will be shown in ds9 region files produced by AE. This label will also be shown in many messages and tables produced by AE. This label does not change as the source is repositioned, whereas the source NAME generated by ae_source_manager will be updated to reflect the source's celestial coordinates. Below are four examples. 1. A FITS catalog produced by the wavdetect tool could provide the source coordinates. idl |& tee -a get_sourcelist.log bt = mrdfits('wavdetect_cat.fits', 1) .run ae ae_source_manager, /ADD, RA=bt.RA, DEC=bt.DEC, PROVENANCE='wavdetect', POSITION_TYPE='wavdetect' exit 2. A ds9 region file (ASCII) written in the "X Y" format using celestial coordinates in decimal degrees could provide the source coordinates. You must read the coordinates in double precision! idl |& tee -a get_sourcelist.log readcol, 'my_catalog.reg', ra, dec, COMMENT=';',F='D,D' .run ae ae_source_manager, /ADD, RA=ra, DEC=dec, PROVENANCE='my_catalog', POSITION_TYPE='user' exit 3. A ds9 region file (ASCII) containing ellipses that represent wavdetect sources and circles that represent sources added by the user could provide the source coordinates, and the PROVENANCE property (to distinguish wavdetect and user-defined sources). idl |& tee -a get_sourcelist.log .run ae ; Read the lines from the region file. reg_filename = 'catalog.cel.reg' num_lines = numlines(reg_filename) lines = strarr(num_lines) openr, regunit1, reg_filename, /GET_LUN readf, regunit1, lines free_lun, regunit1 ; Parse "ellipse" regions for RA & DEC values. ; This regular expression should work for positive or negative DEC, ; space or comma separation between parameters, ; and with or without () in the region specification. result = stregex(lines,'ellipse[^0-9]*([0-9]+\.[0-9]+)[^-0-9]*(-*[0-9]+\.[0-9]+)',/SUB,/EXT) ra = double(reform(result[1,*])) dec = double(reform(result[2,*])) ; Add the ellipse sources to the AE catalog. index = where(ra NE 0 AND dec NE 0) ae_source_manager, /ADD, RA=ra[index], DEC=dec[index], PROVENANCE='wavdetect', POSITION_TYPE='wavdetect' ; Parse "circle" regions for RA & DEC values. result = stregex(lines,'circle[^0-9]*([0-9]+\.[0-9]+)[^-0-9]*(-*[0-9]+\.[0-9]+)',/SUB,/EXT) ra = double(reform(result[1,*])) dec = double(reform(result[2,*])) ; Add the circle sources to the AE catalog. index = where(ra NE 0 AND dec NE 0) ae_source_manager, /ADD, RA=ra[index], DEC=dec[index], PROVENANCE='user', POSITION_TYPE='user' exit 4. An IDL structure built by custom code could provide the source coordinates, provenance information, and source labels. idl |& tee -a get_sourcelist.log restore, /V, 'proposed.target.sav' .run ae ae_source_manager, /ADD, RA=target_cat.ra, DEC=target_cat.dec, LABEL=target_cat.LABEL, PROVENANCE=target_cat.provenance, POSITION_TYPE=position_type exit ===================================================================== CREATE A LIST OF ObsIDS IN THE target, AND A SHELL FOR EACH ObsID cd /data/extract/point_sources.noindex setenv OBS_LIST "" foreach dir (../obs*) if (! -d $dir) continue set obs=`basename $dir | sed -e "s/obs//"` setenv OBS_LIST "$OBS_LIST $obs" end echo $OBS_LIST Make sure all your ObsID's are printed to the screen in the above "echo". Extraction runs for each ObsID in your target can proceed in parallel. It is convenient to create a shell for each ObsID using virtual windows in the Unix "screen" utility. Each shell will store its ObsID number in the environment variable $OBS and for convenience will show the ObsID in its shell prompt. Create a new 'screen' session named after your target (***replace below with your target name***), and create a 'screen' window for each ObsID: setenv TARGET screen -S AE_${TARGET} -d -m -t merge sleep 3 foreach obs ($OBS_LIST) screen -S AE_${TARGET} -X setenv OBS ${obs} screen -S AE_${TARGET} -X setenv PROMPT "${obs} %% " screen -S AE_${TARGET} -X screen -t obs${obs} sleep 1.5 end screen -S AE_${TARGET} -r The "screen" utility accepts commands that begin by holding down CONTROL and pressing the "a" key, which I'll denote by ^a. These commands will help you navigate among the "screen" windows you have just created: ^a" select from list of windows using arrow keys ^an next window ^ap previous window ^aw list the existing windows ^a^c create a new window and switch to it ^aA name the current window The best feature of "screen" is that processes running in "screen" windows continue running if the interactive screen session itself dies. For example, if you have AE processing running on a remote machine under "screen" and connection to that machine from your local machine is lost, then the AE processing will continue. You can easily regain your connection to the running processes by ssh'ing to the remote machine again, then "attaching" to the existing screen windows via: screen -S AE_${TARGET} -r You can find all existing screen sessions via: screen -ls Sometimes, it is convenient to intentionally "detach" from a screen session using this "screen" command: ^ad detach from "screen" session For example, you might detach from your "screen" session when you go home so that from home you can log into your desk machine, attach to the "screen" session, and check on the progress of your processing. Even if you forget to detach before you go home, you can force a detach from home using the -d -r options---see the man page for "screen". When your data processing is completely finished, you can destroy each "screen" window by exiting its shell in the normal way, e.g. via "exit" or "^d". --------------------------------------------------------------------- A note about running on multiple computers: Most tasks in AE will generally run fastest on the computer that is hosting the disk on which the data are stored. However, if you have more ObsIDs than the number of "cores" on that computer, and if your data are visible to other computers (via a networked filesystem), then you will want to spread your ObsID-based processing across multiple computers. You can simply create the screen session shown above on each computer, and choose which ObsIDs will be processed by each computer. ############################################################################# ############# (OPTIONAL) SECTION II: IMPROVING THE CATALOG ################## ############################################################################# The initial catalog (Section I) may be a liberal set of source candidates obtained from liberal Wavdetect runs, human examination of the data, analysis of image reconstructions, etc. In this section, we extract the candidates, judge their significance, and discard those found to be not significant. We also reposition sources, since their original position estimates are often of poor quality. Our 9-step strategy on improving the quality of the catalog is listed below. Obviously, an observer could decide to skip some or all of these steps. #1. Run the initial AE stages in order to construct PSF images and non-overlapping extraction regions. #2. Review the membership of the catalog, i.e. look for missed sources, spurious sources, or grossly mis-positioned sources. At a minimum the review would consist of visual examination of the extraction regions plotted on the X-ray data (in ds9). If desired, the observer might perform image reconstruction on selected crowded regions in order to identify new sources. #3. Remove "unobserved" sources. #4. Prune source candidates that are excessively crowded in all their ObsIDs. #4.5 Decide whether the astrometric alignment among your obsids is acceptable. #5. Search for appropriately-sized background regions. Steps 6-9 are repeated until the catalog is stable. #6. Prune insignificant source candidates using the PROB_NO_SOURCE (called Pb here) statistic computed on an optimal subset of the extractions available for each source. #7. Scan for sources dominated by afterglow events. #8. Improve photometry for sources suffering from photon pileup, so that we can more effectively prune insignificant nearby candidates. #9. Run the AE stages needed to get position estimates. Adjust the source positions using the most appropriate estimate. ===================================================================== SETUP --------------------------------------------------------------------- The processing in this major section requires RMF files only to map energies to PI channels; the actual response is not used. Since RMF files take a long time to calculate, we are going to supply AE with a "generic" RMF (for an arbitrary location on the detector) built by hand. For example, run the command below in one of the ObsID "screen" windows configured above: mkacisrmf infile=CALDB outfile=generic.rmf wmap=none energy="0.3:9.886:0.005" channel=1:1024:1 chantype=PI ccd_id=3 chipx=500 chipy=500 gain=CALDB asolfile=../obs${OBS}/obs.asol obsfile=../obs${OBS}/spectral.evt logfile=none clob+ --------------------------------------------------------------------- Assign a LABEL to each source. If your source detection procedure has already assigned useful source labels that have been passed into ae_source_manager, then skip this step. Otherwise, sort and label the sources in some convenient way. Until the final catalog is determined, you'll be repeatedly visually reviewing the sources in ds9 and probably making notes about individual sources. It may often be necessary to locate a particular source label in ds9. This can be difficult if there are 1000 sources, but will be easier if the sources are sorted in some way. For a convention sort by RA, run idl .run ae ae_source_manager, /SORT_RA ae_source_manager, /SET_LABEL_AS_SEQUENCE Alternatively, the IDL code below will sort and label the catalog into groups of 200 sources lying in annuli around the center of the field. Within each group the sources are sorted by their azimuthal position. With this sorting you can locate a given source, say #520, by first finding the annulus that has sources in the 400-599 range, then scanning azimuthally to find source 520. To guide the eye, ae_source_manager writes the annular regions to the file bullseye.reg which you can load into ds9 whenever needed. idl .run ae ae_source_manager, /SORT_BULLSEYE ae_source_manager, /SET_LABEL_AS_SEQUENCE ===================================================================== CONSTRUCT INITIAL SINGLE-ObsID AE CATALOGS and EXTRACTION APERTURES (Step #1) For EACH ObsID construct an AE "catalog" that produces non-overlapping extraction regions. Choose PIPELINE_RANDOMIZATION=1 (true) or =0 (false) to match your data reduction. CXC-produced event lists have in the past had position randomization turned on, whereas one of the main reasons that we reprocess the data is to turn this randomization off, since it blurs the imaging. We only leave it on for very short (<5 ksec) observations. Since we assume that you are following our data reduction recipe, it is set to 0 in the following example. Below, we save time by using the option /REGION_ONLY to restrict PSF generation to one energy, since we're not yet worried about aperture correction. The /REUSE_NEIGHBORHOOD option specifies that "neighborhood" event files around each source are to be re-used, not re-created, when they exist from a previous pass in order to reduce processing time. Source spectra are extracted from these neighborhoods, NOT directly from the EVTFILE specified below and specified later in ae_standard_extraction. Thus, if you need to modify your event data after ae_make_catalog has been run, then YOU MUST OMIT THE /REUSE_NEIGHBORHOOD option in order for the modifications to be used!! Note that /REUSE_NEIGHBORHOOD is NOT used in Section III where we perform the final extraction processing. The obvious dangers presented by this design are the price we pay for efficient iteration of this catalog improvement loop. The processing of multiple ObsIDs can be run in parallel on different CPUs, using your screen windows defined above. The CIAO/HEASOFT tools spawned by each AE session use different parameter directories and should not interfere with each other. The only files at risk of simultaneous access by multiple AE sessions are the "source.stats" files in the source directories; the chance of problems here is very low. nice idl |& tee -a ae_make_catalog_step1.${OBS}.log .run ae .run ae_make_catalog, getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', /REGION_ONLY, /REUSE_NEIGHBORHOOD, PIPELINE_RANDOMIZATION=0 save, source_not_observed, FILE='../obs'+getenv("OBS")+'/source_not_observed.sav' exit end Due to the ds9 session that IDL spawned, after exiting IDL you'll need to type CTRL^C to get the "tee" process to exit. Using any of your screen windows, review warning or error messages in all log files, ignoring uninteresting ones: egrep -i "WARNING|ERROR|halted" ae_make_catalog.*.log | egrep -v "DISPLAY variable|arithmetic error|no in-band data|No HDUNAME|has no rows" ===================================================================== REVIEW THE CATALOG (Step #2) VISUAL REVIEW IS HIGHLY RECOMMENDED ON EVERY ITERATION THROUGH THIS RECIPE. You can find problems with your catalog, problems with your data, problems with the AE tools, problems with the PSF models, and problems with CIAO. Remember however that you are going to be calculating crowding metrics and source significance statistics below and using them to prune the catalog. Thus, in this visual review you do not need to take on the burden of deciding subjectively which crowded or weak sources should be eliminated. Also, we will shortly remove candidate sources that are not observed in any extraction you have started above (e.g. sources detected on the S-array when you are extracting only I-array data). At a minimum you should visually examine the extraction regions in the ds9 sessions that were created for each ObsID by ae_make_catalog. * Did you miss any obvious sources? * Are the extraction regions reasonably well aligned with the data? --------------------------------------------------------------------- OPTIONAL Single-source Image Reconstruction If desired, you can perform image reconstruction on some of the source neighborhoods in order to identify new sources. Start by listing the source neighborhoods you want to extract into the ASCII file neighborhoods_to_recon.srclist. Then for EACH ObsID, extract the events in a square "neighborhood" around each source. The optional NEIGHBORHOOD_SIZE parameter specifies the width and height of the source neighborhood, in arcseconds. If omitted, neighborhoods are small on-axis and large off-axis. nice idl |& tee -a extract_events.${OBS}.log .run obsname = getenv('OBS') dir = '../obs' + obsname + '/' acis_extract, 'neighborhoods_to_recon.srclist', obsname, dir+'spectral.evt', /EXTRACT_EVENTS, EMAP_FILENAME=dir+'obs.emap', NEIGHBORHOOD_SIZE=70 exit end egrep -i "WARNING|ERROR|halted" extract_events.*.log Then merge the neighborhoods from multiple ObsIDs and perform the reconstructions. If desired, the optional parameter MAXLIKELIHOOD_ITERATIONS will specify the number of reconstruction iterations you wish to perform and the points in the algorithm at which you'd like to save the reconstructed image. In the code below we save images after 10, 20, and 100 iterations. nice idl |& tee -a check_positions.log .run acis_extract, 'neighborhoods_to_recon.srclist', MERGE_NAME='position', /MERGE_OBSERVATIONS, /MERGE_FOR_POSITION, OVERLAP_LIMIT=0.10, ENERGY_RANGE=[0.5,7], /SKIP_SPECTRA, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' save_list = [10,20,100] acis_extract, 'neighborhoods_to_recon.srclist', /CHECK_POSITIONS, MAXLIKELIHOOD_ITERATIONS=save_list, /SKIP_CORRELATION exit end egrep -i "WARNING|ERROR|halted" check_positions.log NOTE: Be aware that the FITS HDU numbering scheme used by ds9 differs from that used by CIAO. For example, suppose you save 3 reconstructed images at various points in the algorithm. Then the file neighborhood.img will have 4 images in it (the data plus three reconstructions). You would look at these 4 FITS HDUs in ds9 with this 0-based syntax: ds9 -log "neighborhood.img[0]" "neighborhood.img[1]" "neighborhood.img[2]" "neighborhood.img[3]" However, to reference one of these images in CIAO you would use a 1-based syntax; for example to copy the last reconstruction you would use: dmcopy "neighborhood.img[4]" recon.img ===================================================================== REMOVE "UNOBSERVED" SOURCES (Step #3) Let's remove and archive any sources that are not observed in any extraction you have started above (e.g. sources detected on the S-array when you are extracting only I-array data). Run this tool only ONCE for the target, NOT for each ObsID! nice idl |& tee -a unobserved_step3.log .run acis_extract, 'all.srclist', /MERGE_OBSERVATIONS, /SKIP_PSF, /SKIP_NEIGHBORHOOD, /SKIP_SPECTRA, /SKIP_TIMING acis_extract, 'all.srclist', COLLATED_FILENAME='temp.collated', LABEL_FILENAME='label.txt' bt = mrdfits('temp.collated', 1) ind = where((bt.NUM_OBS EQ 0), count) if (count GT 0) then begin print, count, F="(%'\nThe %d sources listed in unobserved.srclist were not observed; if desired, use ae_source_manager to remove and archive them.')" forprint, TEXTOUT='unobserved.srclist', SUBSET=ind, bt.CATALOG_NAME, bt.LABEL, F="(%'%s (%s)')", /NoCOMMENT endif end ;If you decide to remove the sources in unobserved.srclist use these commands: .run ae readcol, 'unobserved.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, /REMOVE, NAME=sourcename, TRASH_DIR='unobserved_sources/' exit cntrl-c ===================================================================== PRUNE EXCESSIVELY CROWDED SOURCE CANDIDATES (Step #4) Although we are a long way from spectral analysis, we will extract and re-extract backgrounds throughout this "validation" section (Section II) because two statistics that depend on the background (NET_CNTS and Pb) are used in the pruning decisions in this section. See the sub-section titled "PERFORM STANDARD EXTRACTION ON EACH ObsID" in Section III for a discussion of the BETTER_BACKGROUNDS option used below. !!! NOTE THAT THE AE CALLS IN THIS SECTION ARE NOT IDENTICAL TO THAT FOUND IN Section III !!! --------------------------------------------------------------------- Build decrowded extraction apertures and extract src spectra. Build background regions and extract background spectra. Using your screen windows, for each ObsID perform a standard extraction of the events. This takes many hours (~128 srcs/hr). NOTE: in your first pass through the code below, you can skip the call to ae_make_catalog because you just did that in Step #1 above (IF you didn't alter the catalog by removing unobserved sources in Step #3 above). For all subsequent passes, though, you do need to include this step, because the catalog will have changed due to pruning. nice idl |& tee -a ae_validation_step4.${OBS}.log .run ae .run ae_make_catalog , getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', /REUSE_NEIGHBORHOOD, /REGION_ONLY, PIPELINE_RANDOMIZATION=0 ae_standard_extraction, getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', /REUSE_NEIGHBORHOOD, TIMING=0, GENERIC_RMF_FN='generic.rmf', /BETTER_BACKGROUNDS, NUMBER_OF_PASSES=1 exit end ****** REVIEW THE ARDLIB CONTENTS which AE prints for you to make sure that it's using appropriate files. ****** Due to the ds9 session that IDL spawned, after exiting IDL you'll need to type CTRL^C to get the "tee" process to exit. Using any one of your screen windows (after all the ObsID's have finished), review warning or error messages in all log files, ignoring uninteresting ones: egrep -i "WARNING|ERROR|halted" ae_validation_step4.*.log | egrep -v "DISPLAY variable|no in-band data|No HDUNAME|Program caused arithmetic error|error=|dmimgpick bug|ARF was computed to be zero|has no rows|spans multiple|not observed|spectra will be reused|sources were not extracted|sources not in this observation|ran out of candidate|one emap pixel|Background spectrum|BACKGROUND_IMBALANCE_THRESHOLD" |more For the FIRST time that ae_standard_extraction is run, every source will produce this warning message: WARNING! Background spectrum for ObsID xxxx is missing. The following message is simply warning that a source is severely crowded in a particular ObsID. ae_make_catalog: WARNING! these 35 pairs of sources remain severely crowded (OVERLAP > 0.1). The recipe will soon prune sources that are severely crowded in _all_ the ObsIDs that observed them. For a message like this: ERROR: Source 022606.83+615343.5 (4_542) is not in field of view, but extraction 022606.83+615343.5/446/ exists! you should delete the offending extraction subdirectory (ObsID 446 in this example). This comes about when a source moves a little, just onto or off of a given ObsID. --------------------------------------------------------------------- Adjust the background scaling target for each source. (This tool is adjusting the keywords BKSCL_LO, BKSCL_GL, and BKSCL_HI in source.stats.) Run this tool only ONCE for the target, NOT for each ObsID! Run this tool even if you only have a single ObsID. idl |& tee -a ae_adjust_backscal_range_step4.log .run ae .run ae_adjust_backscal_range, OVERLAP_LIMIT=0.10, MIN_NUM_CTS=100, RERUN_SRCLIST_FILENAME='rerun.srclist' exit end egrep -i "WARNING|ERROR|halted" ae_adjust_backscal_range_step4.log | egrep -v "best extraction" |more If adjustments are recommended, the tool will produce rerun.srclist containing the sources that require new background extractions (most sources will require new backgrounds at the beginning). We will deal with these in Step #5. --------------------------------------------------------------------- Merge in the full energy band to find out which sources suffer overlap in all their ObsIDs. Run this tool only ONCE for the target. Run this tool even if you only have a single ObsID. This takes several minutes. nice idl |& tee -a merge_full_step4.log .run energy_range = [0.5,7] acis_extract, 'all.srclist', MERGE_NAME='Pb_full_band', /MERGE_OBSERVATIONS, MERGE_FOR_PB=0, MIN_NUM_CTS=3, OVERLAP_LIMIT=0.10, ENERGY_RANGE=energy_range, EBAND_LO=energy_range[0], EBAND_HI=energy_range[1], /SKIP_PSF, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'all.srclist', MERGE_NAME='Pb_full_band', COLLATED_FILENAME='tables/Pb_full_band.collated', MATCH_EXISTING='tables/Pb_template', LABEL_FILENAME='label.txt' exit end egrep -i "WARNING|ERROR|halted" merge*.log | egrep -v "DISPLAY variable|no in-band|different value|not merged|accepting only" --------------------------------------------------------------------- Identify sources whose apertures overlap in all ObsIDs. Run this tool only ONCE for the target. idl |& tee -a pruning_overlap.log .r ae .run ae_analyze_pb, /ANALYZE, /IGNORE_PB, GENERIC_RMF_FN='generic.rmf' ae_analyze_pb, /MAKE_REGIONS, /IGNORE_PB exit end wc -l prune.srclist The ae_analyze_pb tool above will list the sources nominated to be pruned in the file prune.srclist, and will spawn a ds9 session with sources color-coded and tagged by their pruning status. IF ONLY A FEW SOURCES WERE NOMINATED TO BE PRUNED, THEN SKIP AHEAD TO THE NEXT MAJOR SECTION ("DECIDE WHETHER THE ASTROMETRIC ALIGNMENT AMONG YOUR ObsIDs IS ACCEPTABLE, Step #4.5"). If you're pruning a lot (e.g. the first time through to remove overlaps), keep going here. --------------------------------------------------------------------- Review the proposed pruning. The CROWDED (orange) and REVIVED (tan) groups contain sources that fail your chosen OVERLAP_LIMIT requirement. Crowded sources generally come in pairs, but it's generally undesirable to prune both sources in such pairs since if one of them was allowed to carry on then it may survive. Whenever both a source and its nearest neighbor are slated for pruning the ae_analyze_pb tool changes the classification of the stronger of the pair to "REVIVED" (tan) and removes it from the proposed pruning list. You may wish to visually confirm that this algorithm made the choices you desire. During your visual review you can reduce clutter by selecting and deleting the ds9 group "cat" (the plus symbols marking catalog positions). You can hide/show all the source labels by toggling the menu item Region:Show Regions Text. NOTE that, when you have multiple observations, AE must choose *which* extraction polygon to display for each source. In an effort to draw attention to sources that are using crowded data AE chooses to display the polygon from the MOST crowded of the observations. --------------------------------------------------------------------- BEFORE WE MODIFY THE SOURCE LIST, create a "pass" directory to archive files related to the AE runs we have completed using the current catalog. We archive a copy of "all.srclist" _before_ we modify it below. setenv PASS pass1 mkdir $PASS cp all.srclist $PASS --------------------------------------------------------------------- Prune invalid sources from the master list "all.srclist". If any of your sources came up "unobserved" and you don't want to keep them around, add them to prune.srclist and they will be removed. idl |& tee -a pruning_overlap.log .run ae .run readcol, 'prune.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, /REMOVE, NAME=sourcename, TRASH_DIR='$PASS/invalid_sources/' exit end Now archive your extraction notes, plus all logs, region files, collated results, etc. In the example below, validation_notes.txt is a textfile containing notes about your data analysis, so change the name as appropriate. cp validation_notes.txt $PASS cp label.txt $PASS mv `'ls' -1 *.srclist | grep -v all.srclist` $PASS mv *.log *.reg *.sav $PASS mv *.collated $PASS mv tables/*.collated $PASS chmod -R a-w $PASS 'ls' -ltr RETURN TO THE BEGINNING OF THIS SECTION (Step #4) to re-extract the modified catalog. ===================================================================== DECIDE WHETHER THE ASTROMETRIC ALIGNMENT AMONG YOUR ObsIDs IS ACCEPTABLE (Step #4.5) At this point, we have an opportunity to very easily estimate the astrometric offsets among your ObsIDs, which directly influences the tightness of the multi-ObsID PSF realized by your on-axis sources. Run the following in a shell in which $OBS_LIST is defined. idl |& tee -a ae_interObsID_astrometry.log .run ae .run ae_interObsID_astrometry, strtrim(strsplit(getenv('OBS_LIST'), /EXTRACT), 2) exit end The report printed shows estimated shifts between every pair of ObsIDs in $OBS_LIST. Each source that is on-axis, not crowded, and reliably detected in both ObsIDs provides a mean data position, with errors, in both ObsIDs and thus provides an estimate of the astrometric offset between the two ObsIDs, with errors. This tool reports an optimally weighted average of those offset estimates (see eq. 5-6 in Bevington). Note that AE's single-extraction position uncertainties decrease as the number of counts extracted grows. Thus, in the tool ae_interObsID_astrometry, offset estimates involving a pair of high-count extractions will be weighted heavily. We can also (optionally) estimate the astrometric offset between each ObsID and a reference catalog that you may have available, such as 2MASS. That requires the machinery of the match_xy.pro tool (distributed with the TARA package); the ae_interObsID_astrometry call would take the following form: idl |& tee -a ae_interObsID_astrometry.log .run match_xy .run ae .run event2wcs_astr = get_astrometry_from_eventlist('../../tangentplane_reference.evt') twomass_cat = build_FITS_cat('../../counterparts/Twomass/2mass_highqual.fits', event2wcs_astr, RA_EXPRESSION='tb.RAJ2000', DEC_EXPRESSION='tb.DEJ2000', RA_ERROR_EXPRESSION='(tb.errmaj+tb.errmin)/2', DEC_ERROR_EXPRESSION='(tb.errmaj+tb.errmin)/2', NAME='TWOMASS') ae_interObsID_astrometry, strtrim(strsplit(getenv('OBS_LIST'), /EXTRACT), 2), ASTROMETRY=event2wcs_astr, REF_CATALOG=twomass_cat exit end NOW IS A GOOD TIME TO CONSIDER WHETHER YOUR OBSID'S ARE HAVE "GOOD ENOUGH" ASTROMETRY---a classic management decision about spending limited resources (to mess with the event data some more) for an un-quantified benefit. Naturally, this astrometry analysis may be more accurate later in this recipe, when spurious candidate sources have been pruned and extraction apertures have been well centered on the sources. We put this quick check here because if you do find an astrometric offset that you cannot live with, then knowing that as soon as possible will avoid wasted work. ===================================================================== SEARCH FOR APPROPRIATELY-SIZED BACKGROUND REGIONS (Step #5) As discussed in the AE manual, when a source has multiple extractions there are good reasons to constrain the relative sizes of the background regions. We are forced to iteratively adjust the allowed range of background region sizes (BACKSCAL). This painful iteration is begun here, BEFORE pruning the catalog with a Pb threshold, because we believe that Pb cannot be accurately calculated when the background is not accurately estimated. An excessively large, not "local", bkg region could overestimate the bkg of a source and produce an erroneously large Pb value. Conversely, an excessively small bkg region would not contain very many bkg counts, forcing Pb to be larger than it would be with a larger bkg sample. The most recent execution of the tool ae_adjust_backscal_range in the previous section will have listed any sources requiring new background extractions in the file rerun.srclist. If that file is empty, then skip ahead to the next section ("PRUNE INSIGNIFICANT SOURCE CANDIDATES"). wc -l rerun.srclist --------------------------------------------------------------------- Re-extract backgrounds for the sources in rerun.srclist Using your screen sessions, run this for each ObsID. nice idl |& tee -a ae_validation_step5.${OBS}.log .run ae .run ae_standard_extraction, SRCLIST_FILENAME='rerun.srclist', getenv("OBS"), EVTFILE='spectral.evt', EXTRACT_SPECTRA=0, TIMING=0, GENERIC_RMF_FN='generic.rmf', /BETTER_BACKGROUNDS, /REUSE_MODELS end Do not exit the IDL session yet; as you iterate through this section you can repeat this call by typing ".go", instead of doing another cut-and-paste. egrep -i "WARNING|ERROR|halted" ae_validation_step5.*.log | egrep -v "DISPLAY variable|no in-band data|No HDUNAME|Program caused arithmetic error|error=|dmimgpick bug|ARF was computed to be zero|has no rows|spans multiple|not observed|spectra will be reused|sources were not extracted|sources not in this observation|ran out of candidate|one emap pixel|previous session|subset of full catalog|BACKGROUND_IMBALANCE_THRESHOLD" |more --------------------------------------------------------------------- Check the background scaling target for each source and adjust if necessary Run this tool only ONCE for the target, NOT for each ObsID! Run this tool even if you only have a single ObsID. idl |& tee -a ae_adjust_backscal_range_step5.log .run ae .run ae_adjust_backscal_range, OVERLAP_LIMIT=0.10, MIN_NUM_CTS=100, RERUN_SRCLIST_FILENAME='rerun.srclist' end Do not exit the IDL session yet; as you iterate through this section you can repeat this call by typing ".go", instead of doing another cut-and-paste. egrep -i "WARNING|ERROR|halted" ae_adjust_backscal_range_step5.log | egrep -v "arithmetic error|excessive OVERLAP" REPEAT THIS SUBSECTION AS LONG AS SIGNIFICANT NUMBERS OF SOURCES REQUIRE ADJUSTMENTS. Go back to the beginning of Step #5 and adjust backgrounds again. I usually move on when I'm down to <50 sources needing adjustment. If the number of sources requiring adjustment stops decreasing, you may find that each source's scaling range is oscillating between two sets of values, and thus not converging. You can detect this problem by examining how the sequence of scaling ranges for individual sources have evolved over time, by running this at the shell: foreach src ( `cat rerun.srclist | cut -f1 -d ' '` ) grep "$src" ae_adjust_backscal_range_step5.log printf '\n\n' end When all sources remaining in rerun.srclist are oscillating then stop trying to adjust them and move on. ===================================================================== LOOK FOR MISSING DATA PRODUCTS Verify that every extraction now has an up-to-date background computed: idl .run readcol, 'all.srclist', sourcename, FORMAT='A', COMMENT=';' src_file_list = file_search(sourcename, 'source.pi') bkg_file_list = file_dirname(src_file_list, /MARK) + 'background_pixels.reg' bkg_file_exists = file_test(bkg_file_list) ind = where(bkg_file_exists EQ 0, count) if (count GT 0) then begin print, 'These better background files are missing.' forprint, TEXTOUT=2, bkg_file_list, SUBSET=ind endif ind = where(bkg_file_exists, count) src_file_list = src_file_list[ind] bkg_file_list = bkg_file_list[ind] print, count, ' bkg regions found' src_file_time = (file_info(src_file_list)).MTIME bkg_file_time = (file_info(bkg_file_list)).MTIME ind = where(src_file_time GT bkg_file_time, count) if (count GT 0) then begin print, 'These better backgrounds are obsolete.' forprint, TEXTOUT=2, bkg_file_list, SUBSET=ind endif print, systime(0,min(bkg_file_time,imin)), bkg_file_list[imin], F="(%'The oldest background was made on %s: %s')" print, systime(0,max(bkg_file_time,imax)), bkg_file_list[imax], F="(%'The newest background was made on %s: %s')" exit end ===================================================================== PRUNE INSIGNIFICANT SOURCE CANDIDATES (Step #6) --------------------------------------------------------------------- MERGE extractions from multiple ObsIDs and calculate Pb in three energy bands. Assuming the target involves more than one ObsID, we're going to merge using AE's /MERGE_FOR_PB screening algorithm which seeks to choose, for each source, the subset of the extractions which minimize Pb in the specified energy band. However, in the code that uses the Pb values to prune the catalog we are going to impose a minimum requirement (3 counts) on SRC_CNTS, to avoid very weak sources even if they have good Pb values. This brings the potential for a nasty trap. A source may have a subset of extractions which produce an acceptable Pb and enough SRC_CNTS, but the _minimum_ Pb occurs using a smaller set of extractions which fails the SRC_CNTS requirement. In short, the Pb pruning must be allowed only after a source has accumulated enough counts to pass our minimum requirement. We communicate that requirement to the MERGE stage via the MIN_NUM_CTS parameter shown below. We will do this MERGE and compute Pb for three energy bands: full (0.5 - 7 keV) soft (0.5 - 2 keV) hard (2.0 - 7 keV) Since we're using the MERGE_NAME option to save each band's data products separately, you can run the three MERGE stages below in parallel if desired. The MERGE stage also prunes extractions whose OVERLAP property exceeds the parameter OVERLAP_LIMIT. That parameter lets you control how much you're willing to tolerate overlapping extraction regions (counts shared among multiple sources). Run the following MERGE stages in parallel, using any 3 of your screen windows. If you're working with just 1 ObsID, you'll have to open a new screen window (ctrl-a ctrl-c to create the window then ctrl-a A to name it) so you'll have 3 total. In the first screen window: nice idl |& tee -a merge_soft.log .run energy_range = [0.5,2] acis_extract, 'all.srclist', MERGE_NAME='Pb_soft_band', /MERGE_OBSERVATIONS, /MERGE_FOR_PB, MIN_NUM_CTS=3, OVERLAP_LIMIT=0.10, ENERGY_RANGE=energy_range, EBAND_LO=energy_range[0], EBAND_HI=energy_range[1], /SKIP_PSF, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'all.srclist', MERGE_NAME='Pb_soft_band', COLLATED_FILENAME='tables/Pb_soft_band.collated', MATCH_EXISTING='tables/Pb_template' exit end In the second screen window: nice idl |& tee -a merge_hard.log .run energy_range = [2.0,7] acis_extract, 'all.srclist', MERGE_NAME='Pb_hard_band', /MERGE_OBSERVATIONS, /MERGE_FOR_PB, MIN_NUM_CTS=3, OVERLAP_LIMIT=0.10, ENERGY_RANGE=energy_range, EBAND_LO=energy_range[0], EBAND_HI=energy_range[1], /SKIP_PSF, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'all.srclist', MERGE_NAME='Pb_hard_band', COLLATED_FILENAME='tables/Pb_hard_band.collated', MATCH_EXISTING='tables/Pb_template' exit end In the third screen window: nice idl |& tee -a merge_full.log .run ae .run energy_range = [0.5,7] acis_extract, 'all.srclist', MERGE_NAME='Pb_full_band', /MERGE_OBSERVATIONS, /MERGE_FOR_PB, MIN_NUM_CTS=3, OVERLAP_LIMIT=0.10, ENERGY_RANGE=energy_range, EBAND_LO=energy_range[0], EBAND_HI=energy_range[1], /SKIP_PSF, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'all.srclist', MERGE_NAME='Pb_full_band', COLLATED_FILENAME='tables/Pb_full_band.collated', MATCH_EXISTING='tables/Pb_template', REGION_FILE='all.reg', LABEL_FILENAME='label.txt' ; Generate a region file that marks where bright PSF hooks would appear in image reconstructions. ; We doubt that sources fainter than 4 counts could survive P_B pruning when situated in the wing of a bright source. ; We omit generating PSF hook regions for extraction more than 5' off-axis, because beyond that the PSF is so big that it swallows up the hook anyway. ae_make_psf_hook_regions, COLLATED_FILENAME='tables/Pb_full_band.collated', HOOK_CNTS_THRESHOLD=4, THETA_THRESHOLD=5 ; Build image reconstructions for the sources that lie close to sources that are generating bright hooks. acis_extract, 'near_psf_hook.srclist', MERGE_NAME='Pb_full_band', /MERGE_OBSERVATIONS, /MERGE_FOR_PB, MIN_NUM_CTS=3, OVERLAP_LIMIT=0.10, ENERGY_RANGE=energy_range, EBAND_LO=energy_range[0], EBAND_HI=energy_range[1], /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'near_psf_hook.srclist', MERGE_NAME='Pb_full_band', /CHECK_POSITIONS exit end In a single shell after all of the above have finished: egrep -i "WARNING|ERROR|halted" merge*.log | egrep -v "DISPLAY variable|no in-band|No HDUNAME|different value|not merged|accepting only" These are ok: merge_hard.log:WARNING: all 4 extractions have excessive OVERLAP; accepting only those with OVERLAP comparable to the best extraction (i.e. OVERLAP<=0.1). The full_band collation above creates the file label.txt that records the mapping from AE source label to sexagesimal source name. ------------------------- ***If your computer dies in the middle of one of the above MERGE or COLLATE stages you can tail the log to see what source was being processed when the computer died, then copy all.srclist and edit the copy to remove everything before the source in question. Then just re-call the stage using your modified sourcelist. If a file got corrupted when the computer died, you'll know right away because the stage will fail on the first source. ------------------------- --------------------------------------------------------------------- Use the Pb statistic to propose a set of sources to be pruned. The code below stratifies the catalog into ranges of Pb, and applies a Pb threshold to identify insignificant sources that should be pruned. Tallies of the sources in each strata are shown in a summary table. You may specify your own choice of Pb boundaries and the Pb threshold that via the optional parameters BOUNDARIES and INVALID_THRESHOLD. idl |& tee -a pruning.log .r ae .run ae_analyze_pb, /ANALYZE, GENERIC_RMF_FN='generic.rmf' ae_analyze_pb, /MAKE_REGIONS exit end ********** If you have converged (no sources left to remove, i.e. no red, magenta, or orange sources in Pbslice.reg or listed in the Pb table produced above), skip down to the end of Step #6 now and read about "convergence". ********** --------------------------------------------------------------------- Review the proposed pruning. Although the automated machinery usually does a good job of estimating backgrounds and identifying insignificant source candidates, an obviously-valid source will sometimes make it into the proposed pruning list (prune.srclist). This occurs when the background region for the source is inappropriately positioned, or contains very few counts. Thus, YOU REALLY SHOULD VISUALLY REVIEW EVERY SOURCE THAT IS SLATED TO BE PRUNED! The ds9 session that comes up at the end of ae_analyze_pb uses color to stratify the catalog by the Pb ranges printed by ae_analyze_pb. You can reduce visual clutter by selecting and deleting the ds9 group "cat" (the plus symbols marking catalog positions). You can hide/show the source labels by toggling the menu item Region:Show Regions Text. NOTE that, when you have multiple observations, AE must choose *which* extraction polygon to display for each source. In an effort to draw attention to sources that are using crowded data AE chooses to display the polygon from the MOST crowded of the observations. Alternatively, there are three ways you could examine *individual* sources that are slated to be pruned (i.e. in prune.srclist): 1. Display Pbslice.reg on the project event list; pan to each source. acis_extract, COLLATED_FILENAME='tables/Pb_full_band.collated', SRCLIST_FILENAME='prune.srclist', /SHOW_REGIONS, DISPLAY='../target.evt', REGION_FILE='Pbslice.reg' 2. Display each extraction of a source in a separate frame. acis_extract, COLLATED_FILENAME='tables/Pb_full_band.collated', SRCLIST_FILENAME='prune.srclist', /SHOW_REGIONS An extra ds9 frame shows the merged neighborhood.evt data, i.e. the subset of the extractions chosen to optimize Pb (in the MERGE named 'Pb_full_band'). (NOTE that in SHOW_REGIONS the background region defined by ae_better_backgrounds is shown by blue pluses that mark the set of emap pixels that define the background region. The blue pluses won't line up with pixel centers in ds9 if the bin size or phase used by ds9 differs from that used in the emap, but you can get a rough idea of the background region.) 3. Display reconstructed images remaining from your most recent repositioning work. acis_extract, COLLATED_FILENAME='tables/Pb_full_band.collated', SRCLIST_FILENAME='prune.srclist', /SHOW_REGIONS, MERGE_NAME='position' Beware, however, that the polygons and catalog positions (+) shown will represent the source positions BEFORE your most recent repositioning work. -------------------------------------------------- Occasionally prune sources whose position and photometry are consistent with the "PSF hook" When your Pb pruning has nearly converged, AFTER you've repositioned sources (Step #9 below), visually review a list of sources (near_psf_hook.srclist) that lie close to sources that are generating bright hooks. In the review, make notes about these potential hook artifacts, including your decisions on whether they should be kept as real or deleted as hook artifacts. It is important to consider the brightness of the source -- in crowded fields, real sources might fall in the PSF hook regions (annular "pandas") of bright nearby sources. We distinguish real sources from hook artifacts by the arbitrary decision that, to be real, a source should have >~2x the counts expected from the PSF hook. idl acis_extract, COLLATED_FILENAME='tables/Pb_full_band.collated', SRCLIST_FILENAME='near_psf_hook.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' exit Annular "panda" regions depict the zones where the CXC expects power from the PSF hook to appear **in image reconstructions**. In the observed data, the hook is of course blurred beyond the boundaries of the panda region by the Chandra PSF. Thus, the best way to visualize the hook is in an image reconstruction, which the tool will display (in the last ds9 frame), when available. Each bright source that has a set of panda regions will also have a text region that reports the number of hook counts expected in each panda region. Beware that these estimates are not corrected for pile-up, and that CXC has only a rough idea of the PSF fraction expected to lie in the hook. Unfortunately, there's no way to indentify which hook count estimate goes with which region .... The color used for a set of panda regions encodes the number of counts expected in the hook across all ObsIDs in the Pb_full_band merge, as follows: hook counts color ----------------------- <5 'green' 5:10 'cyan' 10:20) 'yellow' 20:40) 'red' >40 'magenta' If you decide that a source's position and photometry is consistent with a PSF hook, and you decide to remove the source, then add its telephone number to prune.srclist. Beware that pile-up in the bright central source can make the hook brighter than you might think. Make a list of sources that are likely spurious (really due to the PSF hook of the nearby bright source) and make an executive decision about deleting them. -------------------------------------------------- Occasionally flag sources dominated by afterglows When your Pb pruning has nearly converged, AFTER you've repositioned sources (Step #9 below), scan for sources that are dominated by ACIS afterglow events that slipped through the L1->L2 processing, using the procedure in Step #7 below. Again, don't do this until very near the end of Pb iterations, but when it needs to be done, this is where it will be done. Below we will remove sources with Pb_revised > 0.01 using the file agr_prune.srclist. --------------------------------------------------------------------- Archive this pass BEFORE WE MODIFY THE SOURCE LIST, create a "pass" directory to archive files related to the AE runs we have completed using the current catalog. We archive a copy of "all.srclist" _before_ we modify it below. **Change "pass2" below to be appropriate for whatever pass you're on!** setenv PASS pass2 mkdir $PASS cp all.srclist $PASS --------------------------------------------------------------------- Prune invalid sources from the master list "all.srclist". When you're happy that prune.srclist contains the sources you want to discard, then use the source_manager tool to remove them from all.srclist and to archive the invalid source directories. Also remove afterglow-dominated sources if you ran that step in this pass. idl |& tee -a pruning.log .run ae .run ;Sources with high Pb or other problems: readcol, 'prune.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, /REMOVE, NAME=sourcename, TRASH_DIR='$PASS/invalid_sources/' exit end You can move a source to a better "by-eye" position using the following. This is where you need the source name from prune.srclist. Be sure to keep the D for double precision at the end of the decimal coordinates. ***So change the source NAME, RA, and DEC to be appropriate for your source below!*** idl |& tee -a pruning.log .run ae .run ;Sources that need to be moved: ae_source_manager, /MOVE, NAME='022649.59+621536.0', RA=36.70676D, DEC=+62.26009D, POSITION_TYPE='eye' exit end Remove afterglow-dominated sources, if necessary: idl |& tee -a pruning.log .run ae .run ;Sources with afterglows: readcol, 'agr_prune.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, /REMOVE, NAME=sourcename, TRASH_DIR='$PASS/invalid_sources/' exit end Now archive your extraction notes, plus all logs, region files, collated results, etc. In the example below, validation_notes.txt is a textfile containing notes about your data analysis, so change the name as appropriate. cp validation_notes.txt $PASS cp ../status_filtering_patch_notes.txt $PASS cp label.txt $PASS mv `'ls' -1 *.srclist | grep -v all.srclist` $PASS mv *.log *.reg *.sav $PASS cp $PASS/psf_hook.reg . mv *.collated $PASS mv tables/*.collated $PASS chmod -R a-w $PASS 'ls' -ltr NOW, WE WILL RE-EXTRACT THE MODIFIED CATALOG. THIS STARTS A NEW "PASS." --------------------------------------------------------------------- Build decrowded extraction apertures and extract src spectra. This should all sound familiar -- you've done it before... Using your screen windows, for EACH ObsID: nice idl |& tee -a ae_validation.${OBS}.log .run ae .run ae_make_catalog , getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', /REUSE_NEIGHBORHOOD, /REGION_ONLY, PIPELINE_RANDOMIZATION=0 save, source_not_observed, FILE='../obs'+getenv("OBS")+'/source_not_observed.sav' ; Extract source aperture (no backgrounds yet). ae_standard_extraction, getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', /REUSE_NEIGHBORHOOD, TIMING=0, GENERIC_RMF_FN='generic.rmf', EXTRACT_BACKGROUNDS=0 exit end ****** REVIEW THE ARDLIB CONTENTS which AE prints for you to make sure that it's using appropriate files. ****** In a single shell after all of the above have finished: egrep -i "WARNING|ERROR|halted" ae_validation.*.log | egrep -v "DISPLAY variable|no in-band data|No HDUNAME|Program caused arithmetic error|error=|dmimgpick bug|ARF was computed to be zero|has no rows|spans multiple|not observed|spectra will be reused|sources were not extracted|sources not in this observation|ran out of candidate|one emap pixel" |more --------------------------------------------------------------------- Once your sourcelist is mostly cleaned of insignificant sources, start correcting for photon pile-up. If the last thing you did before starting this pass was to go through the repositioning step (Step #9) for the first time for this target, then you must have cleaned away most of the insignificant sources after several Pb pruning iterations. It is now time to get serious about finalizing your sourcelist and source properties. That means you need to deal with piled-up sources. You should scan for extractions suffering from photon pile-up and apply an algorithm to correct the photometry of those extractions (using the instructions in Step #8). That single-ObsID corrected photometry will replace AE's photometry in the step below, where spatial models of all the sources are constructed to guide the construction of background regions. Accuracy in the single-ObsID photometry allows accuracy in the spatial models, which allows accuracy in the background estimates. Insert HERE the pile-up procedure shown in Step #8 (only execute it through pile-up correction -- you don't have a lightcurve to correct yet). As explained in Step #8, an existing pile-up correction may have to be REPEATED (at this point in the flow) on subsequent passes. So if you have been through repositioning with this target, all subsequent Pb pruning iterations will involve executing Step #8 at this point, to improve pile-up correction upon each iteration. So leave the screen sessions that you started in Step #8 running for now, in case you need them on your next Pb iteration. Obviously you get to skip all this if your first execution of Step #8 reveals that none of your sources are piled up. --------------------------------------------------------------------- Build background regions and extract background spectra. For EACH ObsID: nice idl |& tee -a ae_validation_bkgds.${OBS}.log .run ae .run ; Restore the vector that records which sources are not observed (to save execution time). restore, '../obs'+getenv("OBS")+'/source_not_observed.sav' ; Define and extract backgrounds. ae_standard_extraction, getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', TIMING=0, /REUSE_NEIGHBORHOOD, GENERIC_RMF_FN='generic.rmf', EXTRACT_SPECTRA=0, /BETTER_BACKGROUNDS, NUMBER_OF_PASSES=1 exit end In a single shell after all of the above have finished: egrep -i "WARNING|ERROR|halted" ae_validation_bkgds.*.log | egrep -v "DISPLAY variable|no in-band data|No HDUNAME|Program caused arithmetic error|error=|dmimgpick bug|ARF was computed to be zero|has no rows|spans multiple|not observed|is missing|spectra will be reused|sources were not extracted|sources not in this observation|ran out of candidate|one emap pixel|BACKGROUND_IMBALANCE_THRESHOLD" |more You can verify that your pile-up corrections were used by the BETTER_BACKGROUNDS algorithm with the following grep: egrep -i recon_spectrum ae_validation_bkgds.*.log --------------------------------------------------------------------- Adjust the background scaling range for each source. Run this tool only ONCE for the target, NOT for each ObsID! Run this tool even if you only have a single ObsID. idl |& tee -a ae_adjust_backscal_range.log .run ae .run ae_adjust_backscal_range, OVERLAP_LIMIT=0.10, MIN_NUM_CTS=100, RERUN_SRCLIST_FILENAME='rerun.srclist' end Do not exit the IDL session yet; as you iterate through this section you can repeat this call by typing ".go", instead of doing another cut-and-paste. If adjustments are recommended, the tool will produce rerun.srclist containing the sources that require new background extractions. Re-extract these backgrounds in each ObsID. idl |& tee -a ae_validation_bkgds.${OBS}.log .run ae .run ae_standard_extraction, SRCLIST_FILENAME='rerun.srclist', getenv("OBS"), EVTFILE='spectral.evt', EXTRACT_SPECTRA=0, TIMING=0, GENERIC_RMF_FN='generic.rmf', /BETTER_BACKGROUNDS, /REUSE_MODELS end Do not exit the IDL session yet; as you iterate through this section you can repeat this call by typing ".go", instead of doing another cut-and-paste. Try to repeat the cycle of ae_adjust_backscal_range and ae_standard_extraction runs above until no more sources require adjustment. If the number of sources requiring adjustment stops decreasing, you may find that each source's scaling range is oscillating between two sets of values, and thus not converging. You can detect this problem by examining how the sequence of scaling ranges for individual sources have evolved over time, by running this at the shell: foreach src ( `cat rerun.srclist | cut -f1 -d ' '` ) grep "$src" ae_adjust_backscal_range.log printf '\n\n' end When all sources remaining in rerun.srclist are oscillating then stop trying to adjust them and move on. --------------------------------------------------------------------- The catalog pruning process should be interleaved with one or two recalculations of source positions, by inserting HERE the procedure shown in Step #9. Note that these two operations---catalog pruning and source repositioning---interact. For example: * As crowded pairs of sources are pruned away, the dominant member will often be left at an inappropriate position that should be recalculated. * Pruning a source may allow a neighbor's aperture to grow, changing the set of event used to calculate the centroid position. * Repositioning a source may change its overlap with a neighbor, its background, or the number of counts in its aperture. Any of these changes can force it to become not significant. Thus, pruning and repositioning operations should be interleaved. Since the repositioning step requires lots of resources (both computing and human) we recommend performing it only twice, once when you're down to ~20 sources in prune.srclist, then again when you're back down to ~10 sources in prune.srclist. --------------------------------------------------------------------- ITERATE: Go back to the beginning of Step #6 and repeat the cycle of merging, Pb calculation, pruning, and re-extraction until the catalog is stable. Interleave repositioning (Step #9) and checks for flare-dominated sources (Step #7). This is important because local background estimates for the surviving sources may change when the catalog changes---the act of removing a source from the catalog is a declaration that those counts are no longer considered to be from sources but are considered to be background. When a source's position changes, AE thinks of it as a whole new source, so its extraction region and background get re-extracted. The "new" source may not pass Pb pruning or may affect its neighbors differently. As you iterate through this section, the need for background scaling adjustments will fade away for stable sources. However, sources whose apertures are modified as the catalog is changed will tend to come back into play, requiring more adjustments. CONVERGENCE: How do you know when you can escape the cycle of Pb pruning and achieve nirvana? When the Pb table gives no red, magenta, orange, or tan sources, you're almost there. At that point, take the repositioning road (skip to Step #9 below). With those final positions in hand, start a new Pb pass, looping around through merging and Pb calculation to make sure that the final move didn't cause one of your sources to fail Pb. If you get no red, magenta, orange, or tan sources in the Pb table again, you are well and truly converged and should GO TO STEP #10 to complete this procedure. Yippee! ===================================================================== SCAN FOR SOURCES DOMINATED BY AFTERGLOW EVENTS (Step #7) This is performed occasionally as part of Pb pruning (Step #6). idl ; Identify source extractions suspected to contain afterglow events; sort by Pb calculated with suspect events removed. .r ae .run ae_afterglow_report, 'all.srclist', MERGE_NAME='Pb_full_band', BAND_NUM=0, /SORT_BY_PB exit end idl |& tee agr.log ; Re-run tool on just the suspect sources to produce a cleaner report. .r ae .run ae_afterglow_report, 'agr.srclist', MERGE_NAME='Pb_full_band', BAND_NUM=0, suspect_name, Pb_revised, ag_count, ag_fraction save, suspect_name, Pb_revised, ag_count, ag_fraction, FILE='agr.sav' exit end egrep -i "WARNING|ERROR|halted" agr.log wc -l agr_prune.srclist Sources severely contaminated with afterglow events are recorded in the file agr_prune.srclist, which is later used to prune the catalog . YOU SHOULD REVIEW (in agr.log) EACH SOURCE LISTED IN agr_prune.srclist. Rarely, a bright source will generate enough spurious afterglow identifications to cause it to mistakenly appear in agr_prune.srclist. If you need to, edit agr_prune.srclist to override its decisions about which sources to prune because of afterglows. If agr_prune.srclist doesn't exist, that means there is nothing to prune. The ae_afterglow_report tool prints a table showing each pair of extracted events whose timestamps differ by less than the specified separation, specified in units of CCD frames (max_frame_separation). For each pair the timestamp separation (dfrm), CHIP position offset (dx,dy), event energies (energy), and CCD frame counter (frame MOD 1000) are shown. A classic afterglow sequence would look something like the example below. In this case there is a sequence of 7 events in nearly-consecutive frames (dfrm=1, 2, or 3) at the same location (dx=dy=0) with decreasing event energy. The recalculated Pb is >2%, so we would remove this source. 022611.36+620743.5/Pb_full_band/ (2_1485): 6 suspected afterglow events (46.2% of 13 SRC_CNTS, 230.9% of background); P_b if suspects were removed: 0.02144 dfrm dx dy energy frame MOD 1000 1 0 0 4306-> 2283 649 3 0 0 2283-> 996 650 1 0 0 996-> 852 653 2 0 0 852-> 659 654 2 0 0 659-> 555 656 2 0 0 555-> 461 658 Bright or flaring sources will of course falsely trigger this simple algorithm since real pairs of X-rays often have small separation in time and space. For such a source there will typically be a large number of suspicious pairs, e.g. 115 suspects as shown in the example below. 182031.97-161030.4 We suspect 115 of 456 in-band events (25.2%) are afterglows. More than 20 suspect events; skipped printing! ===================================================================== IMPROVE PHOTOMETRY FROM PILED EXTRACTIONS (Step #8) This step, part of the iterative Pb pruning (Step #6), should be first performed shortly after your first execution of the repositioning step (Step #9). This step may need to be repeated in subsequent pruning passes, as described below. Events are lost from a source suffering significant photon pileup, biasing its photometry downward. Accurate photometry for piled extractions affects our assessment of the validity of other source candidates. When the wing OR READOUT STREAK of a piled extraction contaminates the aperture of a source candidate, AE attempts to model the bright extraction (including its readout streak) and then build a background estimate for the contaminated aperture that accounts for this contamination. When such a model is based on piled photometry, the background estimate will be biased downward, making the contaminated source appear more significant that it really is. Here we try to improve the photometry of piled extractions by constructing a simulated event list that approximates the events ACIS would have detected if pile-up effects were absent, and then performing photometry on that simulated data. Pile-up is detected and corrected on a single-ObsID basis, since pile-up of the same source can vary significantly among observations (and even within a single observation, but we're ignoring that). We start by identifying extractions that may suffer significant pile-up by thresholding the extraction property RATE_3x3 (an estimate of the OBSERVED count rate in a cell of size 3x3 CCD pixels centered on the source position). Using your screen windows, for EACH ObsID: idl |& tee -a pileup_screening.${OBS}.log .run ae .run ae_pileup_screening, getenv('OBS') exit end This tool prints a table of likely-piled extractions, and prepares the files needed for pile-up reconstruction. If any of the rows (extractions) in the table printed above end with "STATUS=0", then you should stop and contact Patrick Broos!!! There is, of course, no magic threshold on RATE_3x3 where pile-up "begins". We have observed that the pile-up tool has produced an event rate correction of 1.11 for an extraction where RATE_3x3 was 0.075 ct/frame.## Thus, we recommend that you consider pile-up correction for extractions with RATE_3x3 > 0.05 (ct/frame). Your least complex path forward is to simply execute the pile-up correction procedure, pileup_reconstruction.txt, on every extraction identified by the screening above. However, pileup_reconstruction.txt explains that it would be reasonable to skip pile-up work on those extractions that you judge are unlikely to have much effect on the Pb calculation of other source candidates. (Pile-up correction to improve the flux, spectrum, lightcurve, etc. of the bright source itself will be performed later in the photometry recipe.) As long as an extraction does merit pile-up correction, you will need to decide on each pruning/repositioning pass whether the pile-up correction should be REPEATED, using the following guidelines: * If the piled source is repositioned, then the correction should be repeated. * If a crowded neighbor of the piled source is repositioned or pruned, then the aperture and background estimation for the piled extraction may change. We want to propagate those changes to the photometry estimate for the piled source (which is key to finding backgrounds for the close neighbors that we might need to prune). Currently, the only way to recompute that photometry is to repeat the pile-up correction itself. Don't worry---it should run faster the second time. REMEMBER THAT THE PILEUP SCREENING CALL ABOVE is preparing INPUT files for the procedure pileup_reconstruction.txt; repeat the pileup_screening call above in every "pass" in which you are going to execute pileup_reconstruction.txt. When all your pile-up work is finished, RETURN TO THE PLACE IN Step #6 where you jumped to this step. FYI, each pile-up reconstruction you perform will save corrected photometry to a file /EPOCH_XXXX/model.photometry, which will be used in future runs of the background algorithm instead of /EPOCH_XXXX/source.photometry. Thus, in subsequent pruning iterations, neighbors of the piled source will enjoy more accurate estimates of the background they suffer. ## The extraction mentioned above was in NGC 3603: 111509.34-611602.0 (p1_3368 ) in ObsID 0633 (THETA=0.3; PSF_FRAC=0.66): 291 (ct) 0.075 (ct/frame) ===================================================================== IMPROVE SOURCE POSITIONS (Step #9) This is performed occasionally as part of Pb pruning (Step #6), right after a_a_b_r. This exercise is long and detailed and should be done with a clear head. Expect to spend at least half a day on it. --------------------------------------------------------------------- Compute "mean data," "correlation," and "maximum likelihood" position estimates The task of choosing the best position estimate for each source remains a research topic. Our current policy is: 1. For uncrowded on-axis sources, we prefer Mean Data Positions (DATA). 2. For uncrowded off-axis sources, we prefer PSF Correlation Positions (CORR). 3. We prefer Maximum Likelihood Image Reconstruction Positions (ML) when the source is sitting on a sloped background caused by the wings of a bright neighbor, because CORR and DATA positions will try to crawl up the slope. However, it is difficult devise an algorithm that can accurately identify such sources. For now we are identifying DATA and CORR estimates that are likely to be influenced by neighbors by testing for a PSF_FRAC (exposure-weighted average PSF fraction) that is significantly reduced from our nominal PSF fraction (assumed to be 0.90 here). In a target with multiple observations one must consider (for each source) which set of observations are going to be used (merged together) in the process of estimating source position. The MERGE stage has two options that are relevant to optimizing the DATA position estimate for each source: 1. The /MERGE_FOR_POSITION option chooses the set of source extractions that minimizes the formal position error, which is a calculation involving the "spread" of the merged extracted counts (influenced by both the PSFs of the merged observations, and the size of the apertures) normalized by the square root of the number of extracted counts (i.e. the common sigma/sqrt(N) "standard error of the mean" calculation). Minimizing this position error calculation thus trades off PSF size against the number of counts extracted. 2. The OVERLAP_LIMIT option prunes severely crowded extractions before they are even considered by the position error minimization above. We take the view that extractions with severe overlap (i.e. those failing an OVERLAP_LIMIT criterion) are unlikely to produce accurate photometry, since a large fraction of extracted counts are shared with a neighbor's extraction. If the photometry is suspect, then it seems reasonable to also suspect the results of the /MERGE_FOR_POSITION algorithm, which uses photometry information. Thus, we conclude that whenever /MERGE_FOR_POSITION is enabled one should also impose a reasonable OVERLAP_LIMIT criterion. Recall that we have three types of position estimates available (DATA, CORR, ML); what sort of merge is best for each type? Clearly, the /MERGE_FOR_POSITION option is appropriate for DATA positions, since it directly minimizes the formal error on DATA position. The CORR and ML position algorithms are very different from the mean DATA algorithm, in that they make use of the entire neighborhood (not just the extracted counts). Nevertheless, it seem reasonable to assume that the observational characteristics favored by /MERGE_FOR_POSITION---tight PSFs and large signal---are also beneficial to the correlation and image reconstruction techniques. Thus, we conclude that all three position estimates will be computed on a data set merged with both the /MERGE_FOR_POSITION option and the OVERLAP_LIMIT criterion. Using that merge, choose which type of position estimate is preferred for each source and produce the three sublists need_data.srclist, need_corr.srclist, and need_recon.srclist which we will use in separate AE sessions. nice idl |& tee -a check_positions_subdivide.log .r ae .run acis_extract, 'all.srclist', MERGE_NAME='position', /MERGE_OBSERVATIONS, /MERGE_FOR_POSITION, OVERLAP_LIMIT=0.10, ENERGY_RANGE=[0.5,7], /SKIP_SPECTRA, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'all.srclist', MERGE_NAME='position', COLLATED_FILENAME='tables/position.collated', LABEL_FILENAME='label.txt' ae_improve_positions, /CHOOSE_PREFERRED_METHOD exit end egrep -i "WARNING|ERROR|halted" check_positions_subdivide.log | egrep -v "DISPLAY variable|no in-band|No HDUNAME|different value|not merged|accepting only" To save time, we compute data positions only for the sources for which we're going to prefer the data position (rate = 0.028 minutes per source on Hiawatha) and we compute correlation positions only for the sources for which we're going to prefer the CORR position (rate = 0.040 minutes per source on Hiawatha). Sometimes ML positions fail, so for those sources we compute CORR positions as well. ML positions used to take MUCH longer (rate = 0.73 minutes per source on Hiawatha) but Pat has recently sped up the code by reducing the size of the neighborhood being reconstructed and reducing the size of the PSF image (may need to beware of this, look for ringing or other artifacts in the recon). Run the following stages in parallel, using any 3 of your screen windows. nice idl |& tee -a check_positions_data.log .run ; Compute DATA positions. acis_extract, 'need_data.srclist', MERGE_NAME='position', /CHECK_POSITIONS, /SKIP_RECONSTRUCTION, /SKIP_CORRELATION exit end nice idl |& tee -a check_positions_corr.log .run ; Compute CORR positions. acis_extract, 'need_corr.srclist', MERGE_NAME='position', /CHECK_POSITIONS, /SKIP_RECONSTRUCTION exit end nice idl |& tee -a check_positions_ml.log .run ; Compute ML (and CORR) positions. acis_extract, 'need_recon.srclist', MERGE_NAME='position', /CHECK_POSITIONS exit end ~~~~~~~~~OPTIONAL~~~~~~~~~ If you're in a hurry you can take the trouble to segment the source lists and divide the work among multiple CPUs. ÊThe code below divides each of the three sourcelists into 8 segments. Ê idl .r ae number_of_cpus = 8 ae_split_srclist, number_of_cpus, SRCLIST_FILENAME='need_recon.srclist', 'need_recon' ae_split_srclist, number_of_cpus, SRCLIST_FILENAME='need_corr.srclist' , 'need_corr' ae_split_srclist, number_of_cpus, SRCLIST_FILENAME='need_data.srclist' , 'need_data' exit Then execute the shell commands suggested by ae_split_srclist to create a separate screen window for each segment writing. In each of those screen windows, paste the following commands to run an IDL session (with its own log file) that processes the segment identified by the shell variable $SEGMENT: nice idl |& tee -a check_positions${SEGMENT}.log Ê .run segment = getenv('SEGMENT') acis_extract, Ê'need_data'+segment+'.srclist', MERGE_NAME='position', /CHECK_POSITIONS, /SKIP_RECONSTRUCTION, /SKIP_CORRELATION acis_extract, Ê'need_corr'+segment+'.srclist', MERGE_NAME='position', /CHECK_POSITIONS, /SKIP_RECONSTRUCTION acis_extract, 'need_recon'+segment+'.srclist', MERGE_NAME='position', /CHECK_POSITIONS exit end ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Êegrep -i "WARNING|ERROR|halted" check_positions*.log | egrep -v "DISPLAY variable|arithmetic error|No HDUNAME" You may see this warning, which occurs when the ML algorithm could not calculate a position: ÊWARNING! ÊBrightest recon pixel near source is not a local maximum; ML position estimate skipped. ------------------------------------------- OPTIONAL In any pass, you can easily and quickly check the astrometry of your event data. ------------------------------------------- Run the following in a shell in which $OBS_LIST is defined. idl |& tee -a ae_interObsID_astrometry.log .run match_xy .run ae .run event2wcs_astr = get_astrometry_from_eventlist('../../tangentplane_reference.evt') twomass_cat = build_FITS_cat('../../counterparts/Twomass/2mass_highqual.fits', event2wcs_astr, RA_EXPRESSION='tb.RAJ2000', DEC_EXPRESSION='tb.DEJ2000', RA_ERROR_EXPRESSION='(tb.errmaj+tb.errmin)/2', DEC_ERROR_EXPRESSION='(tb.errmaj+tb.errmin)/2', NAME='TWOMASS') ae_interObsID_astrometry, strtrim(strsplit(getenv('OBS_LIST'), /EXTRACT), 2), ASTROMETRY=event2wcs_astr, REF_CATALOG=twomass_cat exit end --------------------------------------------------------------------- Decide which source position estimate will actually be *used* for each source, based on the preferences established above and on the availability of the estimates. *** Note that the ae_improve_positions tool will plot the distributions of the proposed changes in RA and DEC, for each of the three types of position estimates, and print the plots to ae_improve_positions_ra.ps and ae_improve_positions_dec.ps. *** nice idl |& tee -a check_positions_choose.log .r ae .run acis_extract, 'all.srclist', MERGE_NAME='position' , COLLATED_FILENAME='tables/position.collated' acis_extract, 'all.srclist', MERGE_NAME='Pb_full_band', COLLATED_FILENAME='tables/Pb_full_band.collated', MATCH_EXISTING='tables/Pb_template' ae_improve_positions, /CHOOSE_POSITION end ; Review the histograms showing position changes to look for something systematic. exit The ae_improve_positions tool above defines a "suspicion metric" to prioritize our visual review (since we'll probably get tired of looking at sources), and creates the files use_catalog.srclist, no_estimate.srclist, use_recon.srclist, use_corr.srclist, and use_data.srclist that order our visual reviews by the suspicion metric. 1. The basic suspicion metric is the distance the source will move, normalized by the aperture size. 2. Sources with a reasonably close neighbor are at some risk that their position CORR and DATA estimates are corrupted, and are thus assigned a higher suspicion than isolated ones. 3. Among the recon sources, we may wish the observer to review the existence of sources with large Pb, so they are assigned a higher suspicion than more significant sources. --------------------------------------------------------------------- Visually review these five groups of sources The following subsections show how to review the position estimates for each group of sources. The following subsections show how to review the position estimates for each group of sources. See recipe.txt for detailed explanations about each subsection. Bring your sourcelists into a text editor, e.g. open use*.srclist open no_estimate.srclist I annotate the relevant sourcelists directly, following each source name and label with a semicolon, then notes about what position that source should use or other pertinent information. I review all the sourcelists making notes about where sources should move (I append a semicolon to every source I review, even if I don't make notes about it, so I'll know which ones I looked at). After all the reviewing is done, I go back and do all the text editing to put sources into different sourcelists if I want to use a different position. I never remove a source from a given sourcelist, I just comment it out once I've added it (to the end) of a different sourcelist. Each subsection below shows two ways the SHOW_REGIONS stage can be used to conduct the visual review. The first displays each extraction of a source in ds9, and the second displays the region file review.reg (from ae_improve_positions) on top of the project-level event list. In these ds9 sessions, region shape encodes position type + = catalog position diamond = mean data position circle = correlation position box = ML recon position and region color encodes the source group COLOR : SRCLIST FILE : CLASS ---------------------------------------------------------------------- yellow: use_eye.srclist : retaining positions previously defined by the observer (FITS keyword POSNTYPE=='eye') red : no_estimate.srclist : MISSING the preferred position estimate green : use_recon.srclist : moving to RECON position blue : use_corr.srclist : moving to CORRELATION position cyan : use_data.srclist : moving to MEAN DATA position tan : use_catalog.srclist : position is stable - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Visually review the BY-EYE POSITIONS (+) you determined before - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - (None will exist if this is your first pass through repositioning). You should look at all of these to make sure that your decision from the last time through is still valid. Ideally, it's good to move these to one of the calculated positions (RECON, DATA, or CORR) if at all possible. idl acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_eye.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' AND/OR (in a separate IDL session): idl acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_eye.srclist', /SHOW_REGIONS, DISPLAY='../target.evt', REGION_FILE='review.reg' If you want to adjust the position again, append a semicolon to the line listing the source and record your new position (in decimal degrees); you will later use that information to manually reposition the source. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Visually review the sources that are missing their preferred position estimates - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - idl acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='no_estimate.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' AND/OR (in a separate IDL session): idl acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='no_estimate.srclist', /SHOW_REGIONS, DISPLAY='../target.evt', REGION_FILE='review.reg' Catalog positions (+) are shown for all sources. You should investigate why these sources have no position estimate available. * ML positions (boxes) will be missing for sources that report this error message: "WARNING! Brightest recon pixel near source is not a local maximum; ML position estimate skipped." * CORR positions (circles) will be missing for sources that report this error message: "'WARNING! Correlation position cannot be computed; peak lies on search boundary.'" * DATA positions (diamonds) should never be missing. For crowded sources missing their ML positions, this is often because they are trying to climb up the wings of their bright neighbor. Sometimes we choose a secondary recon peak for the source position or even delete the source altogether if there is no recon peak to associate with it. This is accomplished with ae_soure_manager, /MOVE or ae_source_manager, /REMOVE -- this will be shown below. Another good idea is to see how the source in question fares in the recon of a nearby source (e.g. the bright one it's trying to climb). YOU SHOULD MAKE A DECISION ON WHAT POSITION TO ADOPT FOR EACH SOURCE (RECON, DATA, CORR, current catalog, or a by-eye estimate) and record that decision by annotating no_estimate.srclist. DS9 note: you can position the cursor where you want then type "c" and ds9 will spawn a new window showing you that position. Set WCS coords to "degrees" and you'll have what you need to assign a "by-eye" position. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Visually review the RECON POSITIONS (box) - - - - - - - - - - - - - - - - - - - - - acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_recon.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' AND/OR (in a separate IDL session): acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_recon.srclist', /SHOW_REGIONS, DISPLAY='../target.evt', REGION_FILE='review.reg' The proposed repositioning you are reviewing is from the cross (catalog) to the box (reconstruction). WHEN THE PIXELS IN THE RECONSTRUCTED IMAGE ARE TOO SMALL, YOU WILL SEE MULTIPLE PEAKS THAT ARE TOO CLOSE TO BE RESOLVABLE SOURCES. I SUGGEST USING A LITTLE SMOOTHING IN DS9 TO PRODUCE A LOCAL MAXIMUM, AND THEN DEFINING A "BY-EYE" POSITION AT THAT PEAK. The file use_recon.srclist consists of two blocks of sources, each sorted by the magnitude of the proposed shift. The first block (suspicion > 10) are low-significance sources (large PROB_NO_SOURCE). When you get tired of looking at the first block, advance to the first source in the second block by entering its label; then use the "r" (rejoin the sequence) command to get the SHOW tool to advance down the list from there. You will see that some weak sources in this list seem to be positioned far from any recon peak in the _current_ reconstruction. Presumably the "recon_detect" reconstructions that produced these candidate sources did have peaks at these positions. We believe this is simply an illustration of the fact that low-level features in reconstructions will vary with changes in the bin phase (or bin size), or with changes in the PSF. If you find a source that should NOT use the RECON position, then append a semicolon followed by your preferred position type (DATA, CORR, current catalog) or a new by-eye set of coordinates (in decimal degrees). - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Visually review the CORR POSITIONS (circle) - - - - - - - - - - - - - - - - - - - - - - acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_corr.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' AND/OR (in a separate IDL session): acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_corr.srclist', /SHOW_REGIONS, DISPLAY='../target.evt', REGION_FILE='review.reg' The proposed repositioning you are reviewing is from the cross (catalog) to the circle (correlation). The file use_corr.srclist consists of two blocks of sources, each sorted by the magnitude of the proposed shift. The first block (DISTANCE_REG2REG < SRC_RAD) are sources that are at some risk that their position estimates are corrupted by a reasonably close neighbor. When you get tired of looking at the first block, advance to the first source in the second block by entering its label; then use the "r" (rejoin the sequence) command to get the SHOW tool to advance down the list from there. You will find that for some weak sources an extraction region at the proposed correlation position looks like it would encompass *fewer* counts than the current extraction region. Often in such cases the CORR position seems to be moving toward a "hole" in the background. We have not investigated this phenomenon, but we suspect that it is an artifact caused by the finite extent of the PSF image interacting with the "hole" in the background. Suppose the correlation shifts a little towards the "hole" such that a few background counts on the opposite edge of the PSF footprint stop participating in the correlation and no new counts are added to the calculation. The loss of those counts will have only a tiny effect on the numerator of the correlation, because the PSF was tiny at the edge of the footprint. However the denominator of the correlation appears to include a term that depends on the *total number* of counts under the PSF footprint; thus the denominator will decrease and the correlation may increase. In short, the correlation calculation rewards shifts that avoid counts where the PSF is very small, and thus a "hole" in the background tends to attract the correlation position. You can override these proposed positions, or not ... it's hard to say what's most correct. If you find a source that should NOT use the CORR position, then append a semicolon followed by your preferred position type (RECON, DATA, current catalog) or a new by-eye set of coordinates (in decimal degrees). - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Visually review the DATA POSITIONS (diamond) - - - - - - - - - - - - - - - - - - - - - - - acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_data.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' AND/OR (in a separate IDL session): acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='use_data.srclist', /SHOW_REGIONS, DISPLAY='../target.evt', REGION_FILE='review.reg' The proposed repositioning you are reviewing is from the cross (catalog) to the diamond (mean data). The file use_data.srclist consists of two blocks of sources, each sorted by the magnitude of the proposed shift. The first block (DISTANCE_REG2REG < SRC_RAD) are sources that are at some risk that their position estimates are corrupted by a reasonably close neighbor. When you get tired of looking at the first block, advance to the first source in the second block by entering its label; then use the "r" (rejoin the sequence) command to get the SHOW tool to advance down the list from there. If you find a source that should NOT use the DATA position, then append a semicolon followed by your preferred position type (RECON, CORR, current catalog) or a new by-eye set of coordinates (in decimal degrees). We sometimes find that sources extracted with ~90% extraction regions (thus about to be assigned a mean data position) are in fact pretty crowded with nearby sources. You can do a reconstruction of those sources if you want. To do recon's on extra sources: Make a small sourcelist, extra_recons.srclist, containing the sources that need reconstruction. nice idl |& tee -a check_positions.log .run acis_extract, 'extra_recons.srclist', MERGE_NAME='position', /CHECK_POSITIONS acis_extract, COLLATED_FILENAME='tables/position.collated', SRCLIST_FILENAME='extra_recons.srclist', /SHOW_REGIONS, /OMIT_BKG_REGIONS, REGION_FILENAME='psf_hook.reg' exit end Then copy these sources to the appropriate use_*.srclist. Generally they should move to their recon positions, so they should be appended to use_recon.srclist. Be sure to comment them out of use_data.srclist if you want them to use a different position. - - - - - - - - - - - - Modify your sourcelist to reflect the decisions you have made. - - - - - - - - - - - - During your reviews you should have recorded your decisions to change AE's recommended position update by adding notes to rows in use_catalog.srclist, no_estimate.srclist, use_recon.srclist, use_corr.srclist, and use_data.srclist. Now you must implement your decisions by "moving" sources to the correct srclist file as needed. We implement such a "move" by COPYING the source to the correct file and then COMMENTING out (with ';') the row where the source originally appears. Any sources with new by-eye source positions should be listed as new rows at the end of use_eye.srclist. This is an organizational aid that will help when you're building the IDL calls to move sources below. - - - - - - - - - - - - Verify your sourcelists - - - - - - - - - - - - When you think that the source lists use_eye.srclist use_recon.srclist use_data.srclist use_corr.srclist use_catalog.srclist are correct, verify that together they contain the correct number of uncommented sources. egrep -v '^;|^[:space:]*$' all.srclist | wc -l egrep -v '^;|^[:space:]*$' use_catalog.srclist use_recon.srclist use_corr.srclist use_data.srclist use_eye.srclist | wc -l The second of these calls may produce fewer sources than the first; the difference should be the sources that you're about to remove (below). --------------------------- Archive the current catalog BEFORE WE MODIFY THE SOURCE LIST, create a "pass" directory to archive files related to the AE runs we have completed using the current catalog. We archive a copy of "all.srclist" _before_ we modify it below. Since RA and Dec are stored in source.stats, we archive those files at this stage as well, since we are about to move everyone's positions. setenv PASS pass???? mkdir $PASS cp all.srclist $PASS 'rm' list foreach src (`egrep -v '^;|^[:space:]*$' all.srclist | cut -d " " -f1`) echo $src/source.stats >> list end tar -czvf $PASS/source.stats.ztar --files-from list --------------------------------- Apply the position review changes First, add any new sources, then remove sources that you want to delete. Then move the sources in use_recon.srclist, use_corr.srclist, and use_data.srclist. (The POSNTYPE property of these sources will be changed by ae_source_manager below.) Finally, apply the "by-eye" positions that you came up with for the new sources that you just added to use_eye.srclist and any sources that were in use_eye.srclist from a previous pass and that need to be moved again. Use the convention of POSITION_TYPE='eye' so that the ae_improve_positions tool will not override your decision in the future. BE VERY CAREFUL TO USE DOUBLE-PRECISION FOR THE SOURCE COORDINATES (the "D" appended to the numbers in the example below)! Following are some examples, with templates. idl |& tee -a update_positions.log .run ae ;Add sources: ***Be sure to invent new source labels!*** ;For example: ae_source_manager, /ADD, RA=274.59639D, DEC=-16.817775D, LABEL='p1_added1', POSITION_TYPE='eye' ae_source_manager, /ADD, RA=D, DEC=+D, LABEL='p1_a1', POSITION_TYPE='eye' ;Delete sources: ;For example: ae_source_manager, /REMOVE, NAME='022743.30+621218.1', TRASH_DIR='$PASS/invalid_sources/' ae_source_manager, /REMOVE, NAME='', TRASH_DIR='$PASS/invalid_sources/' ;Apply any (DOUBLE PRECISION) by-eye positions you may have generated, e.g.: ;For example: ae_source_manager, /MOVE, NAME='022649.59+621536.0', RA=36.70676D, DEC=+62.26009D, POSITION_TYPE='eye' ae_source_manager, /MOVE, NAME='', RA=D, DEC=+D, POSITION_TYPE='eye' ;Reposition the uncommented sources that need to move to an AE-estimated position: .run readcol, 'use_recon.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, MERGE_NAME='position', /UPDATE_POSITIONS_RECON, NAME=sourcename readcol, 'use_corr.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, MERGE_NAME='position', /UPDATE_POSITIONS_CORR, NAME=sourcename readcol, 'use_data.srclist', sourcename, FORMAT='A', COMMENT=';' ae_source_manager, MERGE_NAME='position', /UPDATE_POSITIONS_DATA, NAME=sourcename exit end egrep -i "WARNING|ERROR|halted" update_positions.log OK errors are (e.g.): ERROR! Could not read 022523.73+615953.5/position/source.stats if a source got pruned before these /UPDATE_POSITIONS runs were done. Now archive your notes, plus all logs, region files, collated results, e.g.: cp validation_notes.txt $PASS cp ../status_filtering_patch_notes.txt $PASS cp label.txt $PASS mv `'ls' -1 *.srclist | grep -v all.srclist` $PASS mv *.log *.reg *.sav *.ps $PASS cp $PASS/psf_hook.reg . mv *.collated $PASS mv tables/*.collated $PASS chmod -R a-w $PASS 'ls' -ltr Although you jumped to this step (#9) from near the end of Step #6, you cannot return to that point because a re-extraction is necessary after moving sources above. Thus, you should return to the middle of Step #6 where the following comment is found: << NOW, WE WILL RE-EXTRACT THE MODIFIED CATALOG. THIS STARTS A NEW "PASS." >> =========================================================================================== PREPARE CATALOG DATA PRODUCTS TO RELEASE (Step #10) --------------------------------------------------------------------- Archive your final work from Step #6 **Change "pass2" below to be appropriate for whatever pass you're on!** setenv PASS pass2 mkdir $PASS cp all.srclist $PASS cp ../status_filtering_patch_notes.txt $PASS cp validation_notes.txt $PASS cp label.txt $PASS mv `'ls' -1 *.srclist | grep -v all.srclist` $PASS mv *.log *.reg *.sav $PASS cp $PASS/psf_hook.reg . mv *.collated $PASS mv tables/*.collated $PASS chmod -R a-w $PASS 'ls' -ltr --------------------------------------------------------------------- SORT THE CATALOG BY RA idl |& tee -a sort_catalog.log .run ae .run ae_source_manager, /SORT_RA exit end --------------------------------------------------------------------- Re-generate final data products using the sorted catalog We assume here that you completed Step #6 with a Pb calculation, i.e. that the three Pb merges correspond to the final extraction apertures. Although you may have performed a repositioning step near the end of Step #6, the position **uncertainty** estimates are stale, because you *moved* some sources after doing a MERGE_FOR_POSITION. Thus, we must repeat that MERGE_FOR_POSITION now and collate the results. nice idl |& tee -a sorted_collations.log .run ae .run acis_extract, 'all.srclist', MERGE_NAME='position', /MERGE_OBSERVATIONS, /MERGE_FOR_POSITION, OVERLAP_LIMIT=0.10, ENERGY_RANGE=[0.5,7], /SKIP_SPECTRA, /SKIP_TIMING, GENERIC_RMF_FN='generic.rmf' acis_extract, 'all.srclist', MERGE_NAME='position' , COLLATED_FILENAME='tables/position.collated', REGION_FILE='all.reg', LABEL_FILENAME='label.txt' acis_extract, 'all.srclist', MERGE_NAME='Pb_full_band', COLLATED_FILENAME='tables/Pb_full_band.collated', MATCH_EXISTING='tables/Pb_template' acis_extract, 'all.srclist', MERGE_NAME='Pb_soft_band', COLLATED_FILENAME='tables/Pb_soft_band.collated', MATCH_EXISTING='tables/Pb_template' acis_extract, 'all.srclist', MERGE_NAME='Pb_hard_band', COLLATED_FILENAME='tables/Pb_hard_band.collated', MATCH_EXISTING='tables/Pb_template' ; Rebuild Pbslice.sav for the RA sort. ae_analyze_pb, /ANALYZE, GENERIC_RMF_FN='generic.rmf' exit end grep "{cat}" all.reg > xray_positions.reg grep polygon all.reg > polygons.reg egrep -i "WARNING|ERROR|halted" sorted_collations.log | egrep -v "No HDUNAME|different value" --------------------------------------------------------------------- Now construct a FITS table that contains whatever important catalog information that you wish to release to collaborators. idl |& tee -a catalog_release.log .run restore, 'Pbslice.sav' spawn, 'dmcopy "tables/position.collated[cols CATALOG_NAME,LABEL,RA,DEC,ERX_DATA,ERY_DATA,ERR_DATA,POSNTYPE]" temp.collated clob+' bt_position = mrdfits('temp.collated',1, theader) failure = (n_elements(bt_position) NE num_sources) || (total(CATALOG_NAME NE strtrim(bt_position.CATALOG_NAME,2)) GT 0) if failure then begin print, 'ERROR! Pb and position collations were done on different source lists!' retall endif ; Build a single table suitable for releasing the catalog. cat = replicate(create_struct(bt_position[0], 'Pb_min',0., 'Pb_full_band',0., 'Pb_soft_band',0., 'Pb_hard_band',0.), num_sources) struct_assign, bt_position, cat cat.Pb_min = Pb_min cat.Pb_full_band = reform(Pb[0,*]) cat.Pb_soft_band = reform(Pb[1,*]) cat.Pb_hard_band = reform(Pb[2,*]) cat_fn = 'catalog.fits' mwrfits, cat, cat_fn, theader, /CREATE print, 'Wrote '+cat_fn exit end Check that catalog.fits is sensible. In the dmstat call below, make sure that the numbers are reasonable (e.g. RA and Dec min and max span your field). dmstat catalog.fits | more The columns in catalog.fits are: *source name (built from sexagesimal J2000 coordinates) *source label (a unique source identifier built from the pointing number and a running sequence number) *source position in J2000 decimal RA and Dec *radial position error in arcsec * origin of position (from image reconstruction, mean data within the extraction aperture, or PSF correlation) * 4 estimates of P_B, the probability that the source is a background fluctuation. The most significant sources (those that will go in our "primary" sourcelist) have P_B < 0.3%; less significant sources (those that will go in our "tentative" sourcelist) have 0.3% < P_B < 1.0%. For columns that are not numbers, the dmstat call above will return "Column type not supported." That's ok. The only columns that should have any null values are PB_SOFT_BAND and PB_HARD_BAND. --------------------------------------------------------------------- The files you may wish to formally release include: Pbslice.reg all.reg xray_positions.reg polygons.reg final Pbslice table ../target.evt catalog.fits Now, finally, close out this pass, archiving important files. Later passes will symlink to some of these products, so use a standard name for this final catalog pass: setenv PASS pass_final_catalog mkdir $PASS cp all.srclist $PASS cp validation_notes.txt $PASS cp ../status_filtering_patch_notes.txt $PASS cp label.txt $PASS mv `'ls' -1 *.srclist | grep -v all.srclist` $PASS mv *.log *.reg *.sav $PASS cp $PASS/psf_hook.reg . mv *.collated $PASS mv tables/*.collated $PASS mv catalog.fits $PASS chmod -R a-w $PASS 'ls' -ltr For each source still active in all.srclist, write-protect the "named merges" that we have created so far in the source directories so they won't be mistakenly overwritten later: foreach src (`egrep -v '^;|^[:space:]*$' all.srclist | cut -d " " -f1`) chmod -R a-w ${src}/Pb_full_band ${src}/Pb_soft_band ${src}/Pb_hard_band ${src}/position end --------------------------------------------------------------------- Make useful region files accessible from various convenient places in the directory tree ln -s pass_final_catalog/*.reg . pushd .. ln -s point_sources.noindex/all.reg point_sources.noindex/xray_positions.reg point_sources.noindex/polygons.reg . cd .. ln -s extract/point_sources.noindex/all.reg extract/point_sources.noindex/xray_positions.reg extract/point_sources.noindex/polygons.reg . cd counterparts ln -s ../extract/point_sources.noindex/all.reg ../extract/point_sources.noindex/xray_positions.reg ../extract/point_sources.noindex/polygons.reg . popd ============================================================================== CREATE Postscript IMAGES FOR EACH POINTING, perhaps to be used for publication (Step #11) THIS IS UNTESTED, AND MAY TURN OUT TO BE NOT VERY USEFUL!!! In each pointing directory, idl if ~file_test('polygons.reg') then message, 'You must create a symlink to polygons.reg!' .r ae ;; Create a randomly named scratch directory that will not conflict with another instance of AE. ;; Some tools that use these temp filenames (e.g addarf and addrmf) are executed from a directory ;; other than cwd, and thus require that the path to tempdir is absolute, not relative. .run repeat begin session_name = string(random()*1E4, F='(I4.4)') tempdir = 'AE' + session_name +'.noindex/' tempdir = filepath(tempdir, /TMP) endrep until (NOT file_test(tempdir)) file_mkdir, tempdir print, 'Using temporary directory: ', tempdir run_command, /INIT, PARAM_DIR=tempdir ae_send_to_ds9, my_ds9, NAME='single_pointing_'+session_name, OPTION_STRING=' -geometry 1500x1500 -linear -invert' run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s height 1200")') run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s width 1200")') for jj=0,1 do begin if (jj EQ 0) then begin scene = 'iarray.img.smooth.norm' run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s scale limits 0 3E-8")') endif else begin scene = 'central_1.img' run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s scale limits 0 5")') endelse filename = file_search('pointing_*/full_band/'+scene, COUNT=num_pointings) fdecomp, filename, disk, item_path, item_name, item_qual png_filename = item_path+item_name+'.png' for ii=0,num_pointings-1 do begin ;if file_test(png_filename[ii]) then begin ; print, png_filename[ii]+' exists ...' ; continue ;endif ae_send_to_ds9, my_ds9, filename[ii], 'polygons.reg' run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s zoom to 1")') run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s regions select all")') run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s regions color red")') run_command, /QUIET, string(my_ds9, F='(%"xpaset -p %s regions select none")') run_command, /QUIET, string(my_ds9, png_filename[ii], F='(%"xpaset -p %s saveimage png %s")') print, png_filename[ii] endfor ; ii endfor ; jj end ############################################################################# ######### SECTION III: FINAL EXTRACTION ################################# ############################################################################# !!!! COMMANDS IN THIS SECTION HAVE VITAL DIFFERENCES FROM THOSE IN SECTION II !!!! Cut and paste from this section, NOT from your notes from Section II. !!!! !!!! ========================================== PERFORM STANDARD EXTRACTION ON EACH ObsID --------------------------------------------------------------------- Construct PSFs at 5 energies, plus RMFs, for each source. This takes many hours. Using your screen windows, for EACH observation, execute: nice idl |& tee -a ae_standard_extraction.${OBS}.log .run ae .run ae_make_catalog , getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', PIPELINE_RANDOMIZATION=0 save, source_not_observed, FILE='../obs'+getenv("OBS")+'/source_not_observed.sav' ; Extract source aperture (no backgrounds yet). ae_standard_extraction, getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', EXTRACT_BACKGROUNDS=0 exit end In a single shell after all of the above have finished: egrep -i "WARNING|ERROR|halted" ae_standard_extraction.*.log | egrep -v "DISPLAY variable|no in-band data|No HDUNAME|Program caused arithmetic error|error=|dmimgpick bug|ARF was computed to be zero|has no rows|spans multiple|not observed|sources not in this observation|ran out of candidate|one emap pixel|GRATTYPE|smoothed light curve" |more --------------------------------------------------------------------- Complete your work on single-ObsID pile-up reconstructions Identify extractions that may suffer significant pile-up. Using your AE screen windows, for EACH ObsID: idl |& tee -a pileup_screening.${OBS}.log .run ae .run ae_pileup_screening, getenv('OBS') exit end This tool prints a table of likely-piled extractions, and prepares the files needed for pile-up reconstruction. If any of the rows (extractions) in the table printed above end with "STATUS=0", then you should stop and contact Patrick Broos!!! Here, we want to execute the pile-up modeling procedure for EVERY EXTRACTION EXPECTED TO BE PILED (RATE_3x3 > 0.05 ct/frame), so that you will have a pile-up corrected spectrum on-hand ready for spectral fitting. If you still have a pile-up screen session from the validation recipe, then make sure you still have a screen window for every piled extraction---you may have closed some windows for sources that did not require pile-up modeling in the validation recipe. (You can close an entire screen session with the command "^a^\".) For complicated multi-ObsID cases, I save notes on these tests in check_for_pileup.txt, where I listed all sources showing >=0.05 cts/frame in each ObsID. **** NOTE that these pile-up reconstructions (which take place in EPOCH_XXXX/ directories) are NOT propagated into the photometry merge found later in this section. In other words, we are collating (and probably publishing) piled photometry for the bright sources. That's not as bad as it sounds, since broad-band photometry is a very poor way to describe these bright sources; instead their spectra should be carefully fit with models that account for pile-up. Complex pile-up reconstruction for scientific analyses is outside the scope of this recipe. In some cases, you may want to combine similar ObsIDs and reconstruct that merged spectrum. In other cases, you may want to time-slice one ObsID and reconstruct the spectrum from each slice. **** There is, of course, no magic threshold on RATE_3x3 where pile-up "begins". We have observed that the pile-up tool has produced an event rate correction of 1.11 for an extraction where RATE_3x3 was 0.075 ct/frame.## Thus, we recommend that you consider pile-up correction (using the procedure pileup_reconstruction.txt) for extractions with RATE_3x3 > 0.05 (ct/frame). ## The extraction mentioned above was in NGC 3603: 111509.34-611602.0 (p1_3368 ) in ObsID 0633 (THETA=0.3; PSF_FRAC=0.66): 291 (ct) 0.075 (ct/frame) --------------------------------------------------------------------- Build background regions and extract background spectra. Using your screen windows, for EACH observation, execute: nice idl |& tee -a ae_standard_extraction_bkgds.${OBS}.log .run ae .run ; Restore the vector that records which sources are not observed (to save execution time). restore, '../obs'+getenv("OBS")+'/source_not_observed.sav' ; Define and extract backgrounds. ae_standard_extraction, getenv("OBS"), SOURCE_NOT_OBSERVED=source_not_observed, EVTFILE='spectral.evt', TIMING=0, EXTRACT_SPECTRA=0, /BETTER_BACKGROUNDS exit end In a single shell after all of the above have finished: egrep -i "WARNING|ERROR|halted" ae_standard_extraction_bkgds.*.log | egrep -v "DISPLAY variable|no in-band data|No HDUNAME|Program caused arithmetic error|error=|dmimgpick bug|ARF was computed to be zero|has no rows|spans multiple|not observed|sources not in this observation|ran out of candidate|one emap pixel|GRATTYPE|smoothed light curve|BACKGROUND_IMBALANCE_THRESHOLD" |more --------------------------------------------------------------------- Adjust the background scaling range for each source. Run this tool only ONCE for the target, NOT for each ObsID! Run this tool even if you only have a single ObsID. idl |& tee -a ae_adjust_backscal_range.log .run ae .run ae_adjust_backscal_range, OVERLAP_LIMIT=0.10, MIN_NUM_CTS=100, RERUN_SRCLIST_FILENAME='rerun.srclist' end Do not exit the IDL session yet; as you iterate through this section you can repeat this call by typing ".go", instead of doing another cut-and-paste. If adjustments are recommended, the tool will produce rerun.srclist containing the sources that require new background extractions. Re-extract these backgrounds in each ObsID. idl |& tee -a ae_validation_bkgds.${OBS}.log .run ae .run ae_standard_extraction, SRCLIST_FILENAME='rerun.srclist', getenv("OBS"), EVTFILE='spectral.evt', EXTRACT_SPECTRA=0, TIMING=0, GENERIC_RMF_FN='generic.rmf', /BETTER_BACKGROUNDS, /REUSE_MODELS end Do not exit the IDL session yet; as you iterate through this section you can repeat this call by typing ".go", instead of doing another cut-and-paste. Try to repeat the cycle of ae_adjust_backscal_range and ae_standard_extraction runs above until no more sources require adjustment. If the number of sources requiring adjustment stops decreasing, you may find that each source's scaling range is oscillating between two sets of values, and thus not converging. You can detect this problem by examining how the sequences of scaling ranges for individual sources have evolved over time, by running this at the shell: foreach src ( `cat rerun.srclist | cut -f1 -d ' '` ) grep "$src" ae_adjust_backscal_range.log printf '\n\n' end When all sources remaining in rerun.srclist are oscillating then stop trying to adjust them and move on. --------------------------------------------------------------------- Verify that every extraction now has an up-to-date background computed: idl .run readcol, 'all.srclist', sourcename, FORMAT='A', COMMENT=';' src_file_list = file_search(sourcename, 'source.pi') bkg_file_list = file_dirname(src_file_list, /MARK) + 'background_pixels.reg' bkg_file_exists = file_test(bkg_file_list) ind = where(bkg_file_exists EQ 0, count) if (count GT 0) then begin print, 'These better background files are missing.' forprint, TEXTOUT=2, bkg_file_list, SUBSET=ind endif ind = where(bkg_file_exists, count) src_file_list = src_file_list[ind] bkg_file_list = bkg_file_list[ind] print, count, ' bkg regions found' src_file_time = (file_info(src_file_list)).MTIME bkg_file_time = (file_info(bkg_file_list)).MTIME ind = where(src_file_time GT bkg_file_time, count) if (count GT 0) then begin print, 'These better backgrounds are obsolete.' forprint, TEXTOUT=2, bkg_file_list, SUBSET=ind endif print, systime(0,min(bkg_file_time,imin)), bkg_file_list[imin], F="(%'The oldest background was made on %s: %s')" print, systime(0,max(bkg_file_time,imax)), bkg_file_list[imax], F="(%'The newest background was made on %s: %s')" exit end ===================================================================== REVIEW SINGLE-ObsID EXTRACTION RESULTS ======================================= For EACH observation, execute the following /PLOT stages (do these one ObsID at a time to keep all the plots straight). Each plot is automatically saved as a Postscript file; the filename is shown in the IDL plot window. You are free to rescale them (or switch from linear to log on either or both axes) to make them more informative, and then re-save (File->Print). You may want to record information about them in photometry_notes.txt. The number of points in all these plots should be the number of sources in the catalog; this allows you to middle-click on a point to recover the index of that source, in case you need to investigate outliers. When no data are available for a source (e.g. when the source was not extracted in the ObsID you're plotting), its point appears at a "null" location, typically (0,0). idl obsname = getenv('OBS') collated_filename = '../obs' + obsname + '/all.collated' acis_extract, 'all.srclist', obsname, /SINGLE_OBS, VERBOSE=0 , COLLATED_FILENAME=collated_filename acis_extract, '', obsname, /CONSTRUCT_REGIONS, /PLOT, COLLATED_FILENAME=collated_filename These plots show characteristics of the PSF images and extraction apertures. * The 'PSF pixel size' vs. 'Off-axis Angle' plot (saved as psf_pixsize_vs_theta.ps) reminds us that AE chooses larger pixels for its PSF images as we move off-axis, where lower resolution is needed to reasonably sample the growing PSF. * The 'CROPFRAC' vs. 'Off-axis Angle' plot (saved as cropfrac_vs_theta.ps.) reminds us that the footprint of AE's PSF images are "small" (in order to conserve disk space and CPU time), and thus that a few percent of the PSF power falls outside of the PSF image. We believe that higher-than-typical values are a flaw arising from our use of MARX to generate PSFs; when a source dithers across a chip gap its core sees a lot of dead time that its far wings do not see, increasing the fraction of MARX events lying outside the PSF image. * The 'PSF Fraction' vs. 'Off-axis Angle' plot (saved as psffrac_vs_theta.ps) helps you see the level of crowding in your catalog. Reduced apertures are seen as outliers below the main locus. * The 'Offset between catalog position and PSF centroid' plot (saved as psfoffset_vs_theta.ps) reminds us that the centroid of the Chandra PSF is a poor estimator for the source position off-axis. In the plot, we believe that outliers above the main locus arise from PSF images that are truncated by the edge of the field. * The 'Area of Source Aperture' vs. 'Off-axis Angle' plot (saved as aperturearea_vs_theta.ps) reminds us that the areas of (and thus the backgrounds in) off-axis extraction regions are ~100 times larger than on-axis regions. Outliers below the main locus are simply reduced apertures. The plots described above are largely educational. However, experienced AE users are encouraged to continue scanning these plots because an atypical instance of plots like these may be our first indication of something wrong with a data set, with CIAO, or with the AE software. acis_extract, '', obsname, /EXTRACT_SPECTRA, /PLOT, COLLATED_FILENAME=collated_filename These plots show characteristics of the source extraction. * The 'Mean ARF' vs. 'Off-axis Angle' plot (saved as meanarf_vs_theta.ps) depicts a main locus governed by vignetting, and outliers below the main locus corresponding to reduced apertures. * In the 'Fraction of Source Counts on Primary CCD' plot (saved as multiple_ccds.ps), outliers < 1 represent the rare source whose dither pattern touches multiple CCDs. If you care, you can note which sources are affected. I often write down the zero-based indices for strongly-affected sources (those with only 50% of their counts on the primary CCD). * The 'PSF Fraction' plot (saved as psffrac.ps) reminds us that the Chandra PSF is energy-dependent, i.e. that the extraction aperture corresponds to a very different PSF fraction at low (red) and high (blue) energies. * For a target with little overlap among ObsIDs, the 'Catalog/Data Offset' vs. 'Off-axis Angle' plot (saved as dataoffset_vs_theta.ps) reminds us that the position estimates made by AE (e.g. RECON and CORR positions) differ significantly from simple mean data estimates at large off-axis angles. For a target with overlapping ObsIDs, this plot is not very useful because the "data" position is calculated using only the single ObsID you are examining but the "catalog" position often comes from multi-ObsID merged data. acis_extract, '', obsname, /EXTRACT_BACKGROUNDS, /PLOT, COLLATED_FILENAME=collated_filename These plots show characteristics of the background extraction. * The 'Background Surface Brightness' vs. 'Off-axis Angle' plot (saved as bkgd_sb_vs_theta.ps) serves as a sanity check on the scaling of the background extractions. Any correlation with off-axis angle should have an astrophysical explanation. Outliers above the main locus should be sources suffering significant background from the wings of a neighboring source. * The 'Background Surface Brightness' vs. 'Flux in Source Aperture' plot (saved as bkgd_sb_vs_flux.ps) reminds us that the background algorithms do not do a perfect job of excluding power from the point source itself---in this plot you may see a small correlation between the background surface brightness and the flux of the corresponding source. * The 'Scaling of Bkg Extraction' histogram (saved as backscal_hist.ps) shows the relative size of the bkg region and src aperture (BACKSCAL). You may want to display this as a log-log plot. BACKSCAL will often range from <10 to >1000. * The 'Counts in Background Region' histogram (saved as bkgd_counts_hist.ps) may show a peak at 5 counts, because in the struggle among ObsIDs to influence the *single* bkgd scaling range that all ObsIDs must respect, each ObsID will vote to increase the scaling upper-limit if it cannot find 5 or more counts. It will likely show another peak around 100 counts, the target for merged photometry. You may need to change binsize. exit Archive the saved plots. mv *.ps ../obs${OBS}/ ===================================================================== CUSTOM EXTRACTIONS An unimportant weak neighbor can sometimes force the extraction aperture for a bright and important star to be reduced, discarding a large fraction of its signal. In such cases, you may choose to enlarge the extraction aperture(s) for the bright star, either by drawing an asymmetric aperture that avoids the weak neighbor or by letting the aperture encompass the neighbor and ignoring its contamination of the extraction. If the aperture is modified and the "BETTER_BACKGROUNDS" option was used, then you must rebuild the background region for that source. Consult with Pat on the precise procedure to use. ############################################################################# ################### SECTION IV: PHOTOMETRY ################################# ############################################################################# ===================================================================== MERGE PRODUCTS ACROSS ALL ObsIDS AND PERFORM PHOTOMETRY Run this even if your target has just one ObsID. In a single screen window, nice idl |& tee -a merge.log .run acis_extract, 'all.srclist', MERGE_NAME='photometry', /MERGE_OBSERVATIONS, /MERGE_FOR_PHOTOMETRY, OVERLAP_LIMIT=0.10, MIN_QUALITY=0.5 acis_extract, 'all.srclist', MERGE_NAME='photometry', COLLATED_FILENAME='tables/photometry.collated' exit end egrep -i "WARNING|ERROR|halted" merge.log | egrep -v "DISPLAY variable|arithmetic error|no in-band|No HDUNAME|different value|reading FILTER|fitsio" ===================================================================== REVIEW MULTI-ObsID EXTRACTION RESULTS Run this review step even if you only have a single ObsID in this pointing. When only one ObsID has been extracted, a few plots are uninformative. idl acis_extract, '', /MERGE_OBSERVATIONS, /PLOT, COLLATED_FILENAME='tables/photometry.collated' * The 'Flux1' vs. 'Flux2' plot (saved as flux2_vs_flux1.ps) reminds us that AE estimates photon flux in two ways. * The 'SCAL_MAX/SCAL_MIN' plot (saved as backscal_vs_theta.ps) shows the ratio between the largest and smallest BACKSCAL values among the merged ObsIDs. AE tries to keep this ratio < 2, for reasons discussed in the manual. Sources observed at radically different off-axis angles may not be able to achive that goal. Uncrowded sources are usually able to achieve the same BACKSCAL value in every extraction (a ratio of ~1 in this plot). * The 'log Net Counts' histogram (saved as netcts_vs_theta.ps) is best viewed with log scaling of the Y-axis, so that a powerlaw distribution will appear as a line. * The 'Background Surface Brightness' vs. 'Off-axis Angle' plot (saved as bkgd_sb_vs_theta.ps) serves as a sanity check on the scaling of the background extractions. Any correlation with off-axis angle should have an astrophysical explanation. Outliers above the main locus should be sources suffering significant background from the wings of a neighboring source. * The 'Background Surface Brightness' vs. 'Source Flux' plot (saved as bkgd_sb_vs_flux.ps) reminds us that the background algorithms do not do a perfect job of excluding power from the point source itself---in this plot you may see a small correlation between the background surface brightness and the flux of the corresponding source. * In the 'Mean ARF value' histogram (saved as meanarf.ps), the low tail is caused by reduced apertures and by vignetting. * The 'Scaling of Bkg Extraction' histogram (saved as backscal_composite_hist.ps) shows the relative size of the merged bkg region and src aperture (BACKSCAL). You may want to display this as a log-log plot. BACKSCAL will often range from <10 to >1000. * In the 'Counts in Background Region' histogram (saved as bkg_counts_composite_hist.ps), the outliers below the peak near 100 represent sources for which our recipe and background algorithms were not able to achive the goal of 100 counts in the merged background. ===================================================================== REVIEW EXTREME RESULTS FROM THE "BETTER BACKGROUNDS" ALGORITHM Because the "better backgrounds" algorithm is complex, and because statistical concerns do not allow background region SIZES to be independently chosen for each extraction of a specific source, a few sources may have unusual background regions. We believe the observer should visually review some of the extreme cases, with two goals in mind. First, as with all visual reviews, an experienced eye can spot mistakes in data processing and flaws in algorithms or code. Second, examining extreme cases helps to remind the observer that some sources have inherently poor quality background estimates due to difficult observational circumstances. For example: * When a source aperture overlaps with the wings of a WEAK neighbor, the background algorithm will be unable to obtain very many counts sampling that neighbor because there were only a few counts observed from it. The background region will thus be unable to grow very large, and will end up with very few counts. * When the background level in the observation is very low (e.g. in a short observation), the background region will be forced to grow very large, and will thus not be very "local". The code below uses the SHOW stage of AE to review a handful of sources with the most extreme backgrounds. In the IDL terminal window, the first line of each source report includes a phrase describing why the source's background region is of interest. In the ds9 display, blue plusses mark the adopted background region (defined as a set of pixels in the exposure map). Red plusses mark pixels inside the source aperture, which are not allowed to become part of the background region. This review is more educational than operational, meaning that you will seldom if ever decide to invest time in constructing alternative background regions. Often, the human eye cannot find an alternative region that would have been obviously superior. idl acis_extract, COLLATED_FILENAME='tables/photometry.collated', SRCLIST_FILENAME='bkg_review.srclist', /SHOW_REGIONS exit ############################################################################# #################### 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 ############################################################################# ===================================================================== FIT A VARIETY OF SPECTRAL MODELS TO EACH SOURCE Copy the XSPEC scripts distributed with AE to a local directory. cp -R /usr/common/itt/lib/apps/TARA/code/ae/xspec_scripts . Revise the plausible ranges for model parameters found in these XSPEC scripts to be appropriate for your sources. The *_min and *_max variables in the fitting scripts are used as "soft parameter limits" in XSPEC, and are used to identify "limit violations" in our table generator. NH in these fitting scripts is in units of 1X10^22 cm^(-2). xspec_scripts/thermal/tbabs_vapec.xcm: xspec_scripts/powerlaw/tbabs_pow.xcm xspec_scripts/thermal_2T/tbabs_2vapec.xcm --------------------------------------------------------------------- Make a list (xspec.srclist) of the sources we want to fit. In the example below a cut is made on the SRC_SIGNIF source property. nice idl |& tee -a fit_spectra.log .run ; Read summary FITS table. bt = mrdfits('tables/photometry.collated', 1) ; 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 [band_full] BKG_CNTS = bt.BKG_CNTS [band_full] SRC_SIGNIF = bt.SRC_SIGNIF[band_full] fit_ind = where(SRC_SIGNIF GE 2, COMPLEMENT=nofit_ind) forprint, TEXTOUT='xspec.srclist' , bt.CATALOG_NAME, SUBSET= fit_ind, /NoCOMMENT forprint, TEXTOUT='noxspec.srclist', bt.CATALOG_NAME, SUBSET=nofit_ind, /NoCOMMENT dataset_1d, id1, NET_CNTS[fit_ind], XTIT='NET_CNTS', DENSITY_TITLE='sources fit to spectral model (SNR > 2)', BINSIZE=2, BIN_LOCATION=0.5 dataset_1d, id2, BKG_CNTS[fit_ind], XTIT='BKG_CNTS', DENSITY_TITLE='sources fit to spectral model (SNR > 2)', BINSIZE=2, BIN_LOCATION=0.5 ; List the 20 sources in the XSPEC list with the fewest background counts. print, 'Source Name LABEL BKG_CNTS in full-band' forprint, bt.CATALOG_NAME, bt.LABEL, BKG_CNTS, SUBSET=(fit_ind[sort(BKG_CNTS[fit_ind])])[0:19], F='(%"%s %s %d")' end --------------------------------------------------------------------- Fit a 1-temperature thermal plasma model to every source using a variety of initial parameter values. In the xspec_scripts/ directory we have 5 "model changes" scripts (A.xcm, B.xcm, C.xcm, D.xcm, E.xcm) that override the default initial values for NH and kT. We will run fits with each of these initial parameter sets in order to explore the sensitivity of the result on the initial values. We recommend trying some frozen models as well to help you judge whether the data are informative with regard to plasma temperature. For example, a figure in Preibish et al. 2005, ApJS 160,401 suggests two typical temperatures for young stars: 0.86 keV and 2.6 keV. In the xspec_scripts/ directory we have 2 "model changes" scripts (std1.xcm, std2.xcm) that freeze kT at these temperatures. The FIT stage in AE includes a lock-file mechanism that attempts to prevent simultaneous AE sessions from colliding. When one session detects that another is working on the source it will pause and print a message like this: Blocked by lock file 181948.97-161633.2/spectral_models/ae_lock ... However, I have seen deadlock situations arise with this mechanism, e.g. all AE sessions block indefinitely. If that happens, then remove the lock file by hand and normal processing will probably resume. Create a new 'screen' session named after your target, and create a 'screen' window for each fitting session: setenv TARGET screen -S ${TARGET}_XSPEC -d -m -t fitting foreach model (A B C D E std1 std2) screen -S ${TARGET}_XSPEC -X setenv MODEL ${model} screen -S ${TARGET}_XSPEC -X setenv PROMPT "${model} %% " screen -S ${TARGET}_XSPEC -X screen -t ${model} sleep 1.5 end screen -S ${TARGET}_XSPEC -r Remove any lingering AE lock files and run fitting sessions for these model configurations: A B C D E std1 std2 rm */photometry/spectral_models/ae_lock AE fitting sessions can be run in parallel in the 7 screen windows (corresponding to 7 spectral models) created above. If you have more than 7 processor cores available, you may wish to split up the source list into segments and run them in parallel. For example, if you have 3 machines each with 7 cores available then each machine could have a screen session that works with 1/3 of the source list, as shown below. idl .r ae ae_split_srclist, 3, SRCLIST_FILENAME='xspec.srclist' , 'xspec' ; setenv SEGMENT ".00" ; nice idl |& tee -a fit_${MODEL}${SEGMENT}.log ; setenv SEGMENT ".01" ; nice idl |& tee -a fit_${MODEL}${SEGMENT}.log ; setenv SEGMENT ".02" ; nice idl |& tee -a fit_${MODEL}${SEGMENT}.log .run model_name = getenv('MODEL') segment = getenv('SEGMENT') model_changes_filename = 'xspec_scripts/'+model_name+'.xcm' if ~file_test(model_changes_filename) then message, 'CANNOT FIND MODEL_CHANGES FILE!' acis_extract, 'xspec'+segment+'.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm', MODEL_CHANGES_FILENAME=model_changes_filename exit end Any AE lock files remaining indicate sources that were probably not fit with all models: 'ls' -1 */photometry/spectral_models/ae_lock > refit.srclist rm */photometry/spectral_models/ae_lock Examine the AE messages found in the log files to find failures other than an XSPEC timeout, which is handled explictly later. egrep -i "WARNING|ERROR|halted" fit_*.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|status 152|process killed after" If you find any failures, look nearby in the AE log for the path to XSPEC's log file which you can examine to understand the failure. If any non-timeout problems occurred, then reprocess those sources now. The most common thing to go wrong with the fitting scripts is that the search for parameter errors in XSPEC will take a very long time, tripping a timeout built into AE. Make lists for each model of the sources suffering timeout failures, and repeat those fitting runs with the error computation disabled, e.g. foreach model (A B C D E std1 std2) awk -F "[ /]+" '/process killed after/ {print $6}' fit_${model}*.log > ${model}_noerr.srclist end wc -l *_noerr.srclist Run fitting sessions with errors disabled for these model configurations: A B C D E std1 std2 We specify a generous FIT_TIMEOUT since we're willing to spend whatever it takes to perform the fit itself. nice idl |& tee -a fit_${MODEL}_noerr.log .run model_name = getenv('MODEL') model_changes_filename = 'xspec_scripts/'+model_name+'.xcm' if ~file_test(model_changes_filename) then message, 'CANNOT FIND MODEL_CHANGES FILE!' acis_extract, model_name+'_noerr.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, FIT_TIMEOUT=1000, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm', MODEL_CHANGES_FILENAME=[model_changes_filename,'xspec_scripts/noerr.xcm'] exit end egrep -i "WARNING|ERROR|halted" fit_*_noerr.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|Omit calculation" --------------------------------------------------------------------- Collate the 7 models constructed above, determine which sources require special models, and try to predict which model the observer will think is preferred. The tool ae_suggest_spectral_model has a default set of criteria that define when a model is deemed to be astrophysically "extreme". Using the plots created below, and knowlege of sources in your field, ADJUST THE "extreme" DEFINITION AS DESIRED FOR YOUR target. On way a model can be "extreme" is to have NH or kT parameters outside a nominal range, which can be changed by the keyword parameters NH_min, NH_max, kT_min, kT_max. These ranges are generally more restrictive than the parameter ranges used in the fitting scripts, which are intended to define parameter values that we will not accept. nice idl |& tee -a fit_collate.log .r ae .run ae_suggest_spectral_model, /COLLATE ae_suggest_spectral_model, /ANALYZE end --------------------------------------------------------------------- Fit a two-temperature model for bright sources. The ae_suggest_spectral_model tool will create a list of sources that need two-temperature fits. nice idl |& tee -a fit_2T.log .run acis_extract, 'thermal_2T.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, FIT_TIMEOUT=1000, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal_2T/tbabs_2vapec.xcm' exit end egrep -i "WARNING|ERROR|halted" fit_2T*.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|status 152|process killed after" As with other models, the search for parameter errors in XSPEC can take a very long time, tripping a timeout built into AE. We will again make a source list containing the XSPEC failures and repeat the fitting run with the error computation disabled, and refit those sources with errors disabled. awk -F "[ /]+" '/process killed after/ {print $6}' fit_2T*.log > thermal_2T_noerr.srclist wc -l thermal_2T_noerr.srclist nice idl |& tee -a fit_2T_noerr.log .run acis_extract, 'thermal_2T_noerr.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, FIT_TIMEOUT=1000, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal_2T/tbabs_2vapec.xcm', MODEL_CHANGES_FILENAME='xspec_scripts/noerr.xcm' exit end | For any of these models that are fit on only a subset of the sources, | | if you change you mind about the subset you wish to fit then it would | | be useful to first remove the existing models to spare the observer | | from seeing superfluous models during the visual review. For example, | | to remove all the tbabs_2apec models you can use a little shell script | | like this: setenv MODEL nogrp_tbabs_2vapec foreach src (`egrep -v '^;|^[:space:]*$' xspec.srclist | cut -f1 -d " "`) set file=$src/photometry/source.spectra dmlist "${file}[${MODEL}]" block >& /dev/null if ( $status == 0 ) then echo "$file $MODEL" fdelhdu "${file}[${MODEL}]" N Y endif end --------------------------------------------------------------------- ?? Fit the F.xcm model with a very soft initial temperature for sources with very low NH. ?? ?? The ae_suggest_spectral_model tool will create a list of sources that have a very low NH. ?? ?? nice idl |& tee -a fit_NH_low.log ?? .run ?? acis_extract, 'thermal_NH_low.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm', MODEL_CHANGES_FILENAME='xspec_scripts/F.xcm' ?? exit ?? end ?? ?? egrep -i "WARNING|ERROR|halted" fit_NH_low*.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|status 152|process killed after" ?? ?? As with other models, the search for parameter errors in XSPEC can take a very long time, tripping a timeout built into AE. We will again make a source list containing the XSPEC failures and repeat the fitting run with the error computation disabled, and refit those sources with errors disabled. ?? ?? awk -F "[ /]+" '/process killed after/ {print $6}' fit_NH_low*.log > thermal_NH_low_noerr.srclist ?? wc -l thermal_NH_low_noerr.srclist ?? ?? nice idl |& tee -a fit_NH_low_noerr.log ?? .run ?? acis_extract, 'thermal_NH_low_noerr.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, FIT_TIMEOUT=1000, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm', MODEL_CHANGES_FILENAME=['xspec_scripts/F.xcm','xspec_scripts/noerr.xcm'] ?? exit ?? end --------------------------------------------------------------------- Fit a model with kT frozen for very hard sources. The ae_suggest_spectral_model tool will create a list of sources that need a fit with kT frozen at it maximum value. nice idl |& tee -a fit_kT_max.log .run acis_extract, 'thermal_kT_max.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm', MODEL_CHANGES_FILENAME='xspec_scripts/kT_max.xcm' exit end egrep -i "WARNING|ERROR|halted" fit_kT_max*.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|status 152|process killed after" As with other models, the search for parameter errors in XSPEC can take a very long time, tripping a timeout built into AE. We will again make a source list containing the XSPEC failures and repeat the fitting run with the error computation disabled, and refit those sources with errors disabled. awk -F "[ /]+" '/process killed after/ {print $6}' fit_kT_max*.log > thermal_kT_max_noerr.srclist wc -l thermal_kT_max_noerr.srclist nice idl |& tee -a fit_kT_max_noerr.log .run acis_extract, 'thermal_kT_max_noerr.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, FIT_TIMEOUT=1000, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/thermal/tbabs_vapec.xcm', MODEL_CHANGES_FILENAME=['xspec_scripts/kT_max.xcm','xspec_scripts/noerr.xcm'] exit end --------------------------------------------------------------------- Fit a powerlaw model for certain hard sources. The ae_suggest_spectral_model tool will create a list of sources that need a fit with a powerlaw model. nice idl |& tee -a fit_pow.log .run acis_extract, 'pow.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/powerlaw/tbabs_pow.xcm' exit end egrep -i "WARNING|ERROR|halted" fit_pow*.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|status 152|process killed after" As with other models, the search for parameter errors in XSPEC can take a very long time, tripping a timeout built into AE. We will again make a source list containing the XSPEC failures and repeat the fitting run with the error computation disabled, and refit those sources with errors disabled. awk -F "[ /]+" '/process killed after/ {print $6}' fit_pow*.log > pow_noerr.srclist wc -l pow_noerr.srclist nice idl |& tee -a fit_pow_noerr.log .run acis_extract, 'pow_noerr.srclist', MERGE_NAME='photometry', /FIT_SPECTRA, FIT_TIMEOUT=1000, CHANNEL_RANGE=[35,548], /CSTAT, MODEL_FILENAME='xspec_scripts/powerlaw/tbabs_pow.xcm', MODEL_CHANGES_FILENAME='xspec_scripts/noerr.xcm' exit end --------------------------------------------------------------------- Look for problems in the special "skip errors" runs. egrep -i "WARNING|ERROR|halted" fit_*_noerr.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|Omit calculation" NOTE that trash files are left in ~/.xspec/cache/ when XSPEC processes were killed by timeouts. You may wish to remove these files. --------------------------------------------------------------------- Assign the model preferences we've chosen above, and collate those models. WARNING! The ae_default_model_preference tool will not overwrite an existing "best model" decision unless you supply the /FORCE option. However, do NOT use /FORCE if there are existing "best model" decisions that you wish to preserve!!! restore, 'fit_spectra.sav' .run ae .run ;;; REMOVE COMMENTS ON NEXT LINE ONLY IF YOU'RE SURE THERE ARE NO HUMAN-REVIEW RESULTS YOU WANT TO PRESERVE!!!! ;;; ae_default_model_preference, 'xspec.srclist', MERGE_NAME='photometry', chosen_model, /FORCE acis_extract, 'xspec.srclist', MERGE_NAME='photometry', COLLATED_FILENAME='tables/suggested_mdl.collated', MATCH_EXISTING=bt, HDUNAME='BEST_MDL' bt_bm = mrdfits('tables/suggested_mdl.collated',1) if (n_elements(bt_bm) NE num_sources) then message, 'BUG: suggested_mdl.collated is the wrong length!' end egrep -i "WARNING|ERROR|halted" fit_collate.log --------------------------------------------------------------------- Set up a visual review of the models. For convenience we sort the sources by the reason for the review, assigning priorities (0=highest, corresponding to the "bright" reason) for each review reason. Sources which have no known reason to review will have the lowest priority level (6:6.999), and will be sorted by decreasing flux correction. .run class = ['bright','NH_low','kT_max','powerlaw','extreme','frozen','not_best','sensitive','nominal','BUG'] num_classes = n_elements(class) priority = replicate(float(num_classes-1),num_sources) ; Sources in the class "nominal" (those with no "review reason") will be sorted by decreasing flux_correction. flux_correction = bt_bm.FC0P5_8 / bt_bm.F0P5_8 ind = where(~finite(flux_correction), count) if (count GT 0) then flux_correction[ind] = 1E10 ind = where(review_reason EQ '', count) if (count GT 0) then begin sort_ind = ind[reverse(sort(flux_correction[ind]))] priority [sort_ind] = (num_classes-2) + (findgen(count)/count) review_reason[sort_ind] = string(flux_correction[sort_ind], F='(%"flux correction = %7.1e")') endif ; Sources within each other classe will be sorted by SRC_SIGNIF. for level = (num_classes-3), 0, -1 do begin ind = where(strmatch(review_reason, '*'+class[level]+'*'), count) if (count GT 0) then priority[ind] = level if (count GT 0) then begin sort_ind = ind[reverse(sort(SRC_SIGNIF[ind]))] priority [sort_ind] = level + (findgen(count)/count) endif endfor ; Write each priority level to a separate review source list. note = review_reason+'; suggest '+chosen_model for level = 0,(num_classes-1) do begin sourcelist_fn = string(class[level],F='(%"spec_review_%s.srclist")') ind = where(floor(priority) EQ level, count) if (count EQ 0) then file_delete, sourcelist_fn, /QUIET $ else begin ; Sort by real-valued priority. sort_ind = ind[sort(priority[ind])] forprint, TEXTOUT=sourcelist_fn, SUBSET=sort_ind, bt_bm.catalog_name, note, F='(%"%s ; %s")', /NoCOMMENT print, 'Wrote '+sourcelist_fn endelse endfor ; level save, FILE='fit_spectra.sav' end ===================================================================== VISUALLY REVIEW THE MODELS Now the observer has the difficult task of reviewing the models available for each source, possibly choosing a model other than that recommended by the algorithm above. NOTE THAT THE POWERLAW MODEL IS NEVER RECOMMENDED BY THE ALGORITHM. The ae_spectra_viewer tool lets you browse through the sources, examining the spectral models for each source. When you select a model in the list, your decision is recorded in the keyword BEST_MDL in that source's source.spectra file. Regardless of whether the fits were performed with the C-statistic (on un-grouped data) or with the Chi^2 statistic (on grouped data), you are shown two plots: the cumulative un-grouped model and data, plus the grouped model and data. DURING THIS REVIEW, BE SURE ANY NOTES YOU TAKE REFERENCE SOURCES BY THEIR LABELS rather than by the "sequence number" shown in the tool for navigational convenience. Those sequence numbers merely reflect the position of the source in the set being reviewed, not any identifying property of the sources! In the example below we keep a log of the IDL session as a backup record of the observer's choices. wc -l spec_review_bright.srclist spec_review_NH_low.srclist spec_review_kT_max.srclist spec_review_powerlaw.srclist spec_review_extreme.srclist spec_review_frozen.srclist spec_review_not_best.srclist spec_review_sensitive.srclist spec_review_nominal.srclist idl |& tee -a spec_review.log .run ae ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_bright.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_kT_max.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_powerlaw.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_extreme.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_frozen.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_not_best.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_sensitive.srclist' ae_spectra_viewer, 'tables/suggested_mdl.collated', SRCLIST='spec_review_nominal.srclist' If you used the C-statistic (/CSTAT option to AE) you may run into cumulative spectrum plots where the normalization of the model looks "incorrect", i.e. the right-hand end points of the cumulative plots (representing full-band fluxes) are far apart. We do not understand why this happens. It is NOT clear that this behavior is "wrong"; after all the quantity being minimized in a C-stat fit is NOT anything related to these cumulative spectra. We're not comfortable offering scientific advice on how to choose among different models. You might identify counterparts to your X-ray sources in other wavebands, and consider that external information (e.g. JHK magnitudes, colors, Av). If your sources are mostly stars you might adopt a preference for thermal models over powerlaw models. You might attempt to use a statistical model comparison method. Happy fitting! Find the sources whose best fit was marked as "provisional" during the review. acis_extract, 'xspec.srclist', MERGE_NAME='photometry', COLLATED_FILENAME='temp.collated', MATCH_EXISTING='tables/tbabs_vapec_A.collated', HDUNAME='BEST_MDL' bt = mrdfits('temp.collated',1) forprint, TEXTOUT='spec_review_provisional.srclist', bt.CATALOG_NAME, SUBSET=where(bt.PROVISNL), /NoCOMMENT ===================================================================== PERFORM ANY ADDITIONAL FITTING REQUIRED There may be some sources for which all the models run above are not acceptable. * Some may require a different model family, or adjusting parameters normally frozen (e.g. abundances). NOTE THAT THE ABUNDANCES USED IN AE'S THERMAL PLASMA FITTING SCRIPTS ARE SPECIFICALLY CHOSEN TO MODEL LOW-MASS STARS, AND ARE NOT SUITABLE FOR MODELING O-STARS OR OTHER THERMAL SOURCES!!! The "model changes script" xspec_scripts/solar_abundances.xcm can be used to run fits with solar abundances. * Examination of fit logs will sometimes show that a fit is initially reasonable, but "runs off into the weeds" during the error estimation as parameter space is more fully explored. If you're happy accepting the initial local minimum then you may choose to simply disable the error estimation (as discussed earlier). * Fits can sometimes be coaxed to a reasonable result interactively, e.g. usingt the /INTERACTIVE option to the FIT_SPECTRUM stage. Remember that such a fit will not be collated however unless you designate it as the preferred model, via the ae_default_model_preference or ae_spectra_viewer tools. ===================================================================== COLLATE THE CHOSEN MODELS Once a model has been chosen for every source, then collate those "best models": acis_extract, 'all.srclist', MERGE_NAME='photometry', COLLATED_FILENAME='tables/photometry.collated', MATCH_EXISTING=0, HDUNAME='BEST_MDL' You can print the best model for all the sources, laid out 12 to a page, like this: idl bt = mrdfits('tables/photometry.collated',1) findpro, 'acis_extract', DIRLIST=package_dir plot_name = strtrim(bt.CATALOG_NAME,2) + '/spectral_models/' + strtrim(bt.MODEL,2) + '/ldata.ps' ind = where(~file_test(plot_name), count) if (count GT 0) then plot_name[ind] = package_dir[0]+'no_model.eps' cmd = package_dir[0]+'plot_spectra.pl ' + strjoin(plot_name, ' ') spawn, cmd The code above calls a Perl script, plot_spectral.pl, that builds a Latex document (spectra.tex) that presents the XSPEC plots 12 to a page. Latex is run, and the PostScript document spectra.ps is produced. ===================================================================== COMPUTE LUMINOSITIES, EXAMINE DISTRIBUTIONS OF SPECTRAL PARAMETERS At this point it is helpful to compute luminosities from XSPEC fluxes, and to display distributions of spectral parameters. This is easily done by the tool ae_flatten_collation. Examination of the resulting plots and of the messages produced by the tool may suggest spectral models that are suspicious and should be examined more closely. The ae_flatten_collation tool will also identify fit parameters which violate the stated parameter ranges, omitting those from the FITS and ASCII tables is produces. In many such cases you may decide to refit with the offending parameter frozen at its limit. The tool tries to report various problems in the spectral models, for example error flags that XSPEC reported. You may iterate a few times re-fitting problem sources, collating, and running ae_flatten_collation. In the code below, IF you have adopted a specific distance to your target, then supply it via the DISTANCE parameter. cd tables idl |& tee -a flatten.log ; Supply distance in pc. distance = ?? .run ae ae_flatten_collation, COLLATED_FILENAME='photometry.collated', FLAT_TABLE_FILENAME='xray_properties.fits', DISTANCE=distance, /FAST, SORT=0 NOTE that the position error and Pb columns in the table generated here are bogus because the collation we're using here comes only from the 'photometry' merge, not from the Pb and position merges. ############################################################################# #################### SECTION VI: COLLATING RESULTS #################### ############################################################################# ===================================================================== CHECK FOR AFTERGLOW-DOMINATED SOURCES This has to be done again because earlier afterglow scanning involved ObsIDs chosen for Pb-optimized merges, which might be different than the ObsIDs chosen for the spectal analysis we're about to perform. The fraction of in-band afterglow suspects is saved in the keyword AG_FRAC in XXXX/photometry/source.stats. When tables are generated for the catalog, AG_FRAC can be used to flag sources that have significant afterglow contamination. However, care must be taken to ignore AG_FRAC for bright sources, where numerous false positives are expected. Run the tool ae_afterglow_report once to identify the sources that have at least one suspected afterglow. idl ; Identify source extractions suspected to contain afterglow events; sort by Pb calculated with suspect events removed. .r ae .run ae_afterglow_report, 'all.srclist', MERGE_NAME='photometry', BAND_NUM=0, /SORT_BY_FRACTION exit end idl |& tee agr.log ; Re-run tool on just the suspect sources to produce a cleaner report. .r ae .run ae_afterglow_report, 'agr.srclist', MERGE_NAME='photometry', BAND_NUM=0, suspect_name, Pb_revised, ag_count, ag_fraction save, suspect_name, Pb_revised, ag_count, ag_fraction, FILE='agr.sav' exit end egrep -i "WARNING|ERROR|halted" agr.log You don't need to take any action on this now -- the table generator will use the AG_FRAC keyword to flag any substantially afterglow-contaminated sources. ===================================================================== BUILD A FINAL COLLATION OF THE AE DATA PRODUCTS BY COMBINING PB, POSITION, AND PHOTOMETRY/FITTING COLLATIONS The catalog refinement portion of this recipe constructed four separate named merges aimed at estimating different source properties: * The 'Pb_soft_band', 'Pb_hard_band', and 'Pb_full_band' merges estimated Pb in three bands. * The 'position' merge estimated source position. These were earlier combined into the FITS table 'catalog.fits', now archived to pass_final_catalog/. Now, let's collate the photometry and spectral fitting results derived from the 'photometry' merge. Although photometry.collated was built just a few steps ago, DO NOT SKIP THE COLLATION BELOW! We want to pick up the AG_FRAC values stored by the afterglow screening we just did. idl .run acis_extract, 'all.srclist', MERGE_NAME='photometry', COLLATED_FILENAME='tables/photometry.collated', HDUNAME='BEST_MDL' exit end Verify that the same source list, sorted the same way, has been used to produce catalog.fits and photometry.collated. dmlist "pass_final_catalog/catalog.fits[cols CATALOG_NAME]" data > temp1 dmlist "tables/photometry.collated[cols CATALOG_NAME]" data > temp2 diff temp1 temp2 Finally, let's combine information from all five of the AE merges into a final collation. dmpaste pass_final_catalog/catalog.fits "tables/photometry.collated[cols -CATALOG_NAME,-LABEL,-RA,-DEC,-ERX_DATA,-ERY_DATA,-ERR_DATA,-POSNTYPE]" tables/catalog_and_photometry.collation clob+ ===================================================================== CONSTRUCT A MORE CONVENIENT TABLE OF X-RAY PROPERTIES The final collation of AE quantities, catalog_and_photometry.collation, can be improved in several ways: * The table must be "flattened" to be published, and to be easily used by collaborators, i.e. vector columns used to store photometric quantities computed for several energy bands must be converted into multiple scalar columns. * Only three of the 15 energy bands that AE normally reports are commonly of interest. Eliminating the obscure bands simplifies the table. * The photon fluxes that AE estimates have little scientific relevance; adding rough estimates of energy fluxes makes the table far more useful. * The uncertainties on NET_CNTS estimated by AE can be replaced with the more accurate uncertainties estimated by the "aprates" tool in CIAO. IF you have collated some spectral fitting results, then you'll want to supply the distance to the target (in pc) via the optional DISTANCE input to ae_flatten_collation. This tool will produce both FITS and ASCII tables. pushd tables idl |& tee -a flatten.log .run ae .run ae_flatten_collation, COLLATED_FILENAME='catalog_and_photometry.collation', FLAT_TABLE_FILENAME='xray_properties.fits', SORT=0 exit end egrep -i "WARNING|ERROR|halted" flatten.log | egrep -v "arithmetic error|Large number of counts|Maximum number of iterations|source and/or background counts are 0" popd ===================================================================== BUILD LaTeX TABLES for Publication pushd tables idl |& tee -a hmsfr.log .run ae .run ; Supply local path to table templates. findpro, 'acis_extract', DIRLIST=package_dir template_filename = package_dir+'/hmsfr_tables.tex' hmsfr_tables, 'xray_properties.fits', template_filename exit end egrep -i "WARNING|ERROR|halted" hmsfr.log Several LaTeX tables are produced, plus a small LaTeX document "xray_properties.tex" which includes those tables. To view them, make a PDF file, then open it: pdflatex xray_properties open xray_properties.pdf OR gv --orientation=landscape --scale=2 xray_properties.pdf & Both of the above viewing methods are annoying -- the first because you can only rotate one page at a time, the second because gv is just awful. What I got to work is to read xray_properties.pdf into acrobat, rotate the whole document using Tools -> Pages -> Rotate, then save it as xray_properties.pdf again. Then just preview xray_properties.pdf. If you have the misfortune of discovering late in the game that the astrometric frame of your data needs an adjustment, you can use the 2-vector SKY_OFFSET input to the table generator to apply an offset to the coordinates and re-generate the source names. The units of SKY_OFFSET are CIAO sky pixels. The RA & x axes have opposite signs so if you want to increase RA of the sources you supply a negative deltaX value. The DEC and y axes have the same signs. ############################################################################# #################### SECTION VII: STACKING DIM SPECTRA #################### ############################################################################# Given a set of sources on which an AE extraction has been performed, it's easy to ``stack'' them together to produce a multi-ObsID spectrum, ARF, RMF, and photometric quantities. You simply have to ``trick'' AE into thinking that the {\em N} ``observation'' sub-directories you have for the several sources you wish to stack are really {\em N} observations of a single source, e.g.\ one named ``stacked''. Then you can run the MERGE\_OBSERVATIONS stage on the ``stacked'' source to get a combined spectrum, background, ARF, RMF, and photometric quantities. Avoiding photometric bias is very important for stacking work, so we omit the /MERGE\_FOR\_PHOTOMETRY option (\S\ref{sec:discarding_observations}). That leaves the {\em observer} responsible for deciding whether any high-backgrounds extractions (e.g. off-axis) should be omitted! An example of simple shell commands and the AE call to do this is below. Note that back-quote characters are used in several places in the shell commands below. # Create a directory in which to perform your stacking analyses. mkdir stacking cd stacking # Create a symbolic link to the AE point source extraction directory. ln -s /foo/bar/data/extract/point_sources point_sources # Create a subdirectory for a group of sources you wish to stack, and move there. mkdir weak_sources cd weak_sources # By whatever means you choose, list the names of sources you wish to stack in a sourcelist file. For example idl bt = mrdfits('../../point_sources/tables/photometry.collated',1) include_in_stack = ( ... ) forprint, TEXTOUT='stacked.srclist', bt.OBJECT, /NoCOMMENT, SUBSET=where(include_in_stack) # Using symbolic links, make all the extraction directories to be stacked appear in the stacking directory. foreach src (`egrep -v '^;|^[:space:]*$' stacked.srclist`) echo $src foreach spectrum_file (`cd ../point_sources/$src; echo */source.pi`) set obs=`dirname $spectrum_file` ln -s ../point_sources/${src}/${obs} ${src}_${obs} end end # AE expects to find the file "stacked/source.stats" so just copy one from any of the sources you're stacking. cp ../point_sources/${src}/source.stats . # Back in the stacking directory, run an AE merge on the "source" named "weak_sources" and group the spectrum. cd .. echo weak_sources > temp.srclist idl |& tee -a weak_sources/stacked.log .run acis_extract, 'temp.srclist', /MERGE_OBSERVATIONS, OVERLAP_LIMIT=0.10, /SKIP_PSF, /SKIP_NEIGHBORHOOD, /SKIP_TIMING acis_extract, 'temp.srclist', COLLATED_FILENAME='weak_sources/weak_sources.collated' acis_extract, 'temp.srclist', /FIT_SPECTRA, CHANNEL_RANGE=[35,548], MODEL_FILENAME='point_sources/xspec_scripts/plot_gross_and_background_spectra.xcm', SNR_RANGE=[5,10] exit end egrep -i "WARNING|ERROR|halted" weak_sources/stacked.log | egrep -v "DISPLAY variable|arithmetic error|no in-band|No HDUNAME|fitsio|FILTER" Note that in the MERGE\_OBSERVATIONS stage the source spectra are added; the ARFs are averaged (weighted by EXPOSURE); and the EXPOSURE keywords are added (\S \ref{sec:arfs}). For example if you stack 8 sources with 10 counts each from the same 20ks observation then the summed spectrum will have 80 counts, the EXPOSURE keyword will be 160ks, and the multi-ObsID ARF will be the simple average of the 8 ARFs. The flux you derive from fitting with the merged files is the (average) flux {\bf per source}. The merged event list, source.evt, will not be displayed properly in ds9 unless you remove the FITS keywords TDMIN18, TDMAX18, TDMIN19, TDMAX19 ############################################################################# #################### SECTION VIII: DIFFUSE EMISSION #################### ############################################################################# Create a workspace for diffuse analysis: cd /data/extract/ mkdir diffuse_sources.noindex cd diffuse_sources.noindex ===================================================================== CREATE A LIST OF ObsIDS IN THE TARGET, AND A SHELL FOR EACH ObsID As shown in SECTION I (SET UP), list all the ObsIDs in the environment variable $OBS_LIST and create a "screen" window for each. ===================================================================== CONSTRUCT DIFFUSE EVENT DATA FOR EACH OBSERVATION Build suitable event lists named diffuse.evt and corresponding exposure maps named diffuse.emap for EACH observation. After performing a point source extraction (above), you will have point source models for each ObsID (in the files obsXXXX/ae_better_backgrounds.sav) left over from running ae_better_backgrounds. We will use these models and the tool ae_better_masking to perform more intelligent masking than the simple circular mask regions produced by AE. The ae_better_masking tool masks pixels on the sky where the estimated flux from point sources is significant compared to the estimated non-point source flux. Thus, bright stars are masked more heavily than weak ones. Masking is applied to the "validation.evt" event data rather than to the "spectral.evt" data because for diffuse analysis we want to suppress background as much as possible and do not care about any cleaning steps that may be damaging to point sources. --------------------------------------------------------------------- For each ObsID, compare the event data to the point source and readout streak models we have leftover from the point source extraction. You are looking for any obvious features not modeled that should be masked, such as really badly-piled sources and very crowded cores of clusters. foreach obs ($OBS_LIST) echo ${obs} ds9 -title ${obs} ../obs${obs}/validation.evt ../obs${obs}/observation_counts_model.img -scalelims 0 10 -match scales -zoom to fit -match frames wcs & end Any extra masking you desire should be described by drawing regions in ds9 and saving them to extra_maskfile.reg in celestial coordinates. Mask the validation.evt data for EACH ObsID. Using your screen windows, for EACH observation, execute: nice idl |& tee -a ae_better_masking.${OBS}.log .run ae .run obsname = getenv("OBS") obsdir = '../obs' + obsname + '/' ae_better_masking, obsname, EVTFILE_BASENAME='validation.evt', /REUSE_MODELS, /SKIP_EXTRACT_BACKGROUNDS, EXTRA_MASKFILE='extra_maskfile.reg' exit end Cntrl-c Once all ObsID's are done, in one screen window, egrep -i "WARNING|ERROR|halted" ae_better_masking.*.log | egrep -v "DISPLAY variable|arithmetic error|no in-band|error=" This tool first brings up the iarray_1 exposure map with extra_maskfile.reg shown on it (but nothing masked). Later it makes a 4-tile ds9 with the results of the masking. --------------------------------------------------------------------------- Extra masking ****** LOOK at the masked event list and emap that are brought up in ds9! When a source is very bright, its wings will deserve masking all the way out to the edge of the PSF footprint and you will see a square mask on the emap. In such cases, there are presumably pixels beyond the PSF footprint that would have been masked if the PSF footprint was larger. Look closely at the masked event data to see if these remaining wings are bright enough to damage your diffuse analysis. You may also notice crowded regions that probably have more unresolved point source light than diffuse emission. ****** If you find that the masked data contains light that you do not wish to be considered diffuse emission then revise the hand-drawn regions in extra_maskfile.reg and repeat the ae_better_masking runs. If your hand-drawn masking includes detector-oriented features like readout streaks, then you will need independent extra_maskfile.reg files for each ObsID, since aimpoints and roll angles will generally differ; if your masking includes only astrophysical features, then a single extra_maskfile.reg file in celestial coordinates will suffice. If, instead, you implement additional masking by directly modifying the output products of ae_better_masking data, then you MUST do it in a way that ensures that the SAME masking is done to both the emap and event data; we recommend constructing the masked emap and then using dmimgpick to filter the event list to match the emap. FYI, we have a tool named ae_make_streak_regions that will draw rectangles around readout streaks, given a list of sources that produce streaks. Now that ae_better_masking is building models of readout streaks, we hope that drawing masks around streaks will no longer be necessary. Note that the CIAO TOOL acis_streak_map can also identify streaks. --------------------------------------------------------------------------- When you have suitable masked data products, use symbolic links to give them standardized names in the obsXXXX directories: cd /data/extract/ Ê foreach obs ($OBS_LIST) ÊÊÊecho ÊÊÊÊ${obs} ÊÊÊÊ ÊÊÊpushd obs${obs} ÊÊÊln -s background.emap diffuse.emap ÊÊÊln -s background.evt Êdiffuse.evt ÊÊÊpopd Êend ===================================================================== CONSTRUCT "STOWED" EVENT DATA CORRESPONDING TO EACH OBSERVATION Our spectral analysis of diffuse sources makes use of the "stowed background" data sets in CALDB to subtract instrumental background. (See the AE manual.) Such stowed data sets can also be used to subtract instrumental background when smoothing images of the diffuse data. If you are following our standard L1->L2 data reduction recipe then you will already have available in the obsXXXX/ directories symbolic links ("stowed.evt" and "stowed.emap") pointing to the following stowed background data products for EACH observation: 1. An event list covering the CCDs used in your analysis, reprojected onto the tangent plane of the observation and cleaned in a manner corresponding to how your observation was cleaned. 2. An "exposure map" corresponding to the stowed event data. This image has no physical interpretation, but the ratio of the stowed and observation emaps records the scaling of the stowed data required to use it as a model of the instrumental background in our observation. This scaling accounts for differences in exposure time (for each CCD separately) between the observation and the stowed data calibration, and tries to account for differences in background rates in those two data sets (i.e. for "space weather"). --------------------------------------------------------------------- If you do not have these stowed data products, then see Appendix A for suggestions on how to construct them for EACH observation. Note that these stowed event lists and emaps must be UNMASKED (in order to construct diffuse.bkg.scaling below). Be sure that your stowed data products cover the same set of CCDs that appear in diffuse.evt/diffuse.emap! In the end, these data products must be available via the paths "extract/obsXXXX/stowed.evt" and "extract/obsXXXX/stowed.emap". --------------------------------------------------------------------- Apply your diffuse masking to the stowed data products for each observation. 1. The observation's MASKED emap (diffuse.emap) is used to mask the stowed event list (stowed.evt) so it can serve as a background event list (diffuse.bkg.evt) corresponding to the masked observation data (diffuse.evt) when working in image space. 2. The observation's UNMASKED emap (obs.emap) and the stowed data's UNMASKED emap (stowed.emap) are used to build a map (diffuse.bkg.scaling) of the scaling that must be applied to the background data (diffuse.bkg.evt) when working in image space. IT IS VERY IMPORTANT THAT THIS SCALING MAP IS CONTINUOUS, NOT MASKED, BECAUSE IT WILL BE RESAMPLED LATER by build_scene.pro when background images are constructed. If there are holes in this scaling map, then some values in the resampled scaling map will be very small leading to a big spike in the background image! In ONE of your screen windows, cd /data/extract/ nice idl |& tee -a stowed_masking.log .run ae ae_mask_stowed_data exit end ===================================================================== SMOOTH THE DIFFUSE DATA At this point you are set up to use the recipe in tara_smooth.txt (or some other smoothing process that can make use of masked exposure maps) to smooth the diffuse data. ===================================================================== REGION FILES & SOURCE LISTS cd diffuse_sources.noindex By whatever clever means you can find, construct region files that define one or more diffuse ``objects''. The regions must be in celestial coordinates if you have multiple observations, but can by in ``physical'' (Chandra SKY) coordinates if you have only one observation. Each region file can contain multiple or compound components if desired, e.g. a polygon plus a circle. Since you've masked point sources from the event list and exposure map (in the previous step) the diffuse region files do not need to worry about excluding point sources. Construct an ASCII file with two columns to serve as an AE ``diffuse catalog'' that lists all the regions you need to extract: * The first column should be a name for the source (no spaces are allowed). * The second column should be the path to the ds9 region file that defines the source. For example, this simple C-shell script will build such a catalog for all the ``.reg'' files in the current directory, using the filename base as the name of the source: 'rm' diffuse_all.srclist foreach file (*_diffuse.reg) echo "`basename $file .reg` $file" >> diffuse_all.srclist end ===================================================================== CREATE SYMLINKS DEFINING THE EXTRACTION APERTURES FOR THE STOWED BACKGROUND DATA In this diffuse analysis, the stowed background data takes on the role of "background" in AE and in XSPEC. Using symlinks, we force the background region (background.reg) in each extraction to be the same as the observation's aperture (extract.reg). foreach src (`egrep -v '^;|^[:space:]*$' diffuse_all.srclist | cut -d " " -f1`) echo $src foreach obs ($OBS_LIST) echo " $obs" mkdir -p $src/${obs}/ ln -s extract.reg $src/${obs}/background.reg end end ===================================================================== EXTRACT OBSERVED SPECTRA IN DIFFUSE REGIONS EXTRACT STOWED BACKGROUND SPECTRA IN DIFFUSE REGIONS For EACH observation and for each region listed in diffuse_all.srclist, AE uses the diffuse observation data for a "source" extraction and uses the ACIS Stowed Data for a "background" extraction. WARNING! THE SURFACE BRIGHTNESS VALUES YOU DERIVE FROM SPECTRAL FITTING WILL BE WRONG if you do not supply the appropriate EMAP_ENERGY parameter to AE in the calls below! See a discussion on how we define an "area on the sky" for a diffuse source in the Diffuse Sources section of the AE manual and in Broos et al. (2010). We assume that your emap was built for a single energy (supplied as the "monoenergy" parameter to mkinstmap). The shell script below will print that energy for every ObsID---make sure it matches the EMAP_ENERGY value supplied to AE below. foreach obs ($OBS_LIST) echo OBS${obs} dmlist ../obs${obs}/asphist/ccd0.instmap head | egrep "ENERG_LO|ENERG_HI" end You can use multiple IDL sessions in parallel on different processors to speed up the AE runs below, which must be run for EACH observation. The CONSTRUCT_REGIONS stage will open a ds9 session that is used to convert the supplied diffuse regions from celestial to physical (SKY) coordinates. nice idl |& tee -a extract_diffuse.${OBS}.log obs = getenv('OBS') dir = '../obs' + obs + '/' ardlib = dir + 'ardlib.par' pbk_filename = dir + 'obs.pbkfile' msk_filename = dir + 'obs.mskfile' aspect_fn = dir + 'obs.asol' print, obs & print, dir .run acis_extract, 'diffuse_all.srclist', obs, dir+'diffuse.evt', /CONSTRUCT_REGIONS, /DIFFUSE acis_extract, 'diffuse_all.srclist', obs, dir+'diffuse.evt', /EXTRACT_SPECTRA, EMAP_FILENAME=dir+'diffuse.emap', EMAP_ENERGY=1.0, ARDLIB_FILENAME=ardlib, PBKFILE=pbk_filename, MSKFILE=msk_filename, ASPECT_FN=aspect_fn ;WE NEED AN EXPLANATION HERE FOR WHY WE USE THE UNMASKED STOWED DATA PRODUCTS FOR THE BACKGROUND EXTRACTION, RATHER THAN THE MASKED ONES!!! acis_extract, 'diffuse_all.srclist', obs, dir+'stowed.evt', /EXTRACT_BACKGROUNDS, EMAP_FILENAME=dir+'stowed.emap' exit end Look for errors in the logs: egrep -i "WARNING|ERROR|halted" extract_diffuse.*.log | egrep -v "DISPLAY variable|arithmetic error|error=|long time|GRATTYPE|CCD_ID|not observed" This warning message requires action on your part: WARNING! Aborted extraction because WMAP is zero! You should investigate, and then remove inside008/9894/obs.stats so that the MERGE will not see it. For each of these cases, display the ObsID's emap with the source's extraction region. You will probably find that they barely overlap, which explains why the WMAP computed by dmextract was empty. You should discard extractions like this, by removing the extraction directory or removing just the obs.stats file therein. If you have a reasonable number of sources and ObsIDs then you can visually examine the extracted event lists: foreach src (`egrep -v '^;|^[:space:]*$' diffuse_all.srclist | cut -d " " -f1`) echo $src setenv CMD "ds9 -bin factor 8 " foreach obs ($OBS_LIST) if (! -e $src/${obs}/source.evt) continue setenv CMD "$CMD $src/${obs}/source.evt -region $src/${obs}/extract.reg $src/${obs}/background.evt " end eval $CMD & end WE NEED AN EXPLANATION HERE FOR WHY WE USE THE UNMASKED STOWED DATA PRODUCTS FOR THE BACKGROUND EXTRACTION, RATHER THAN THE MASKED ONES!!! ===================================================================== CHECK NORMALIZATION OF STOWED BACKGROUNDS The normalization of the stowed backgrounds was carefully adjusted in the L1->L2 processing for each ObsID (Hickox & Markevitch, 2006, \apj, 645, 95, Section 4.1.3). However, there are plenty of things that can go wrong, so you should verify that the merged background-corrected spectra have flux consistent with zero in the energy band (9-12 keV) on which the scaling was derived. RUN THIS ONCE, not for each ObsID. mkdir tables idl |& tee -a check_scaling.log .run acis_extract, 'diffuse_all.srclist', /MERGE_OBSERVATIONS, EBAND_LO=9.0, EBAND_HI=12.0 acis_extract, 'diffuse_all.srclist', COLLATED_FILENAME='tables/photometry9-12.collated' bt = mrdfits('tables/photometry9-12.collated', 1) sc = bt.SRC_CNTS nc = bt.NET_CNTS nc_error = 0.5 * (bt.NET_CNTS_SIGMA_UP + bt.NET_CNTS_SIGMA_LOW) if (n_elements(bt) EQ 1) then begin help, sc, nc, nc_error endif else begin function_1d, id1, indgen(n_elements(bt)), nc/sc, Y_ERROR=nc_error/sc, PSYM=1, LINE=6, XTIT='source number', YTIT='NET_CNTS/SRC_CNTS (9-12 keV)' function_1d, id2, indgen(n_elements(bt)), nc/nc_error, PSYM=1, LINE=6, XTIT='source number', YTIT='NET_CNTS/NET_CNTS_error (9-12 keV)' dataset_1d , id3, nc/nc_error, XTIT='NET_CNTS/NET_CNTS_error (9-12 keV)' endelse forprint, bt.LABEL end egrep -i "WARNING|ERROR|halted" check_scaling.log | egrep -v "DISPLAY variable|arithmetic error|no in-band|No HDUNAME|keyword not found|FILTER" Drag the function_1d widet wider and change its NSKIP field to 1 (Edit dialog box) so you can see each individual error bar. Ideally, the set of error bars on the 9-12 keV photometry plotted above will be *consistent* with zero (some .... If not, then here are some of the possible problems that come to mind: 1. The diffuse and stowed data have not been cleaned in the same way (e.g. the stowed data was scaled to match validation.evt but you have extracted data derived from spectral.evt, different bad pixel tables were used for the two datasets, etc.). 2. The diffuse and stowed data have not had the same point source masking applied. 3. The diffuse data somehow contains a very heavily piled source that is producing X-ray events in the 9:12 keV band. 4. Other reasons??? ===================================================================== REVIEW SINGLE-ObsID EXTRACTION RESULTS AE has lots of ways to visualize the source properties it estimates. From a data processing perspective it's reasonable to expect that in many cases flaws in the catalog, recipe, recipe execution, AE code, CIAO code, or CALDB might show up as outliers in some of these plots. From a science perspective these plots may facilitate identification of unusual sources and characterization of the source population. For EACH observation, execute: idl acis_extract, 'diffuse_all.srclist', getenv("OBS"), /EXTRACT_SPECTRA, /PLOT acis_extract, 'diffuse_all.srclist', getenv("OBS"), /EXTRACT_BACKGROUNDS, /PLOT ===================================================================== MERGE PRODUCTS ACROSS ALL ObsIDS AND PERFORM PHOTOMETRY mkdir tables nice idl |& tee -a merge.log .run acis_extract, 'diffuse_all.srclist', /MERGE_OBSERVATIONS acis_extract, 'diffuse_all.srclist', COLLATED_FILENAME='tables/photometry.collated', REGION_FILE='all.reg', LABEL_FILENAME='label.txt' exit end egrep -i "WARNING|ERROR|halted" merge.log | egrep -v "DISPLAY variable|arithmetic error|no in-band|No HDUNAME|fitsio|FILTER" ===================================================================== REVIEW MULTI-ObsID EXTRACTION RESULTS idl acis_extract, /SHOW, COLLATED_FILENAME='tables/photometry.collated' acis_extract, '', /MERGE_OBSERVATIONS, /PLOT, COLLATED_FILENAME='tables/photometry.collated' Two styles of plots are produced. Some show the distribution of a single source property, and others are scatter plots for two source properties that may reveal useful insights into the data. (When only one ObsID has been extracted, a few plots are uninformative.) (Need to add plot showing SCAL_MAX/SCAL_MIN) !!!!!!!!!!!! May 2011 !!!!!!!!!!!!!!!! Development version of AE recently suffered significant changes (version 3923) which have not been adequately tested. Be sure do to the SHOW stage above! !!!!!!!!!!!! May 2011 !!!!!!!!!!!!!!!! ===================================================================== CALCULATE THE GEOMETRIC AREAS OF THE DIFFUSE REGIONS We must compute the geometric areas of the diffuse regions in order to convert surface brightness estimates to fluxes integrated over the diffuse regions. AE cannot calculate these geometric areas for two reasons: 1. The observer is free to define a diffuse region that is compound, i.e. an ellipse "minus" a circle, since CIAO will accept such regions as spatial filters. AE does not have the ability to parse such compound region files, and thus cannot calculate an area analytically. 2. IF the diffuse region was guaranteed to lie fully within some ObsID, then AE could calculate the geometric area of the exposure map that lies withing the region. However, in a multi-pointing target, a diffuse region could easily be larger than the ACIS field of view. Thus, at no point in the extraction will AE have a FITS image that fills the diffuse region. Thus, our only choice is to force the observer to supply a "project scene image", i.e. some image that spans all the diffuse regions. Each diffuse region can then be applied as a spatial filter on this image, using dmcopy, and the area of the result can be calculated. A Chandra SKY (PHYSICAL in ds9) coordinate system (corresponding to any ObsID) must be defined on this scene image. In our own workflow, we generate a project-spanning multi-ObsID exposure map corresponding to the diffuse event data (i.e. an emap that represents the masks that were applied to remove the point sources), which serves nicely for computing region areas below. idl |& tee geometric_areas.log .run ae .run ; Specify the name of a "scene" image that has a FOV the encompasses all the diffuse regions. ; The pixel values in this image are irrelevant. scene_fn = '../adaptive_smoothing/iarray.diffuse.emap' ; Read the 2-column sourcelist, which as both source names and paths to defining region files in celestial coordinates. readcol, 'diffuse_all.srclist', sourcename, catalog_region_fn, FORMAT='A,A', COMMENT=';' ; Trim whitespace and remove blank lines. sourcename = strtrim(sourcename,2) ind = where(sourcename NE '', num_sources) sourcename = sourcename[ind] print, num_sources, F='(%"\n %d sources found in catalog.")' if (num_sources EQ 0) then retall ;; Create a randomly named scratch directory that will not conflict with another instance of AE. repeat begin session_name = string(random()*1E4, F='(I4.4)') tempdir = 'AE' + session_name +'.noindex/' tempdir = filepath(tempdir, /TMP) endrep until (NOT file_test(tempdir)) file_mkdir, tempdir help, tempdir run_command, /INIT, PARAM_DIR=tempdir flat_map_fn = tempdir + 'flat.img' temp_region_fn = tempdir + 'temp.reg' temp_image_fn = tempdir + 'temp.img' ; Create a map of the scene filled with 1's. scene = readfits(scene_fn, header) flat_map = make_array(/INTEGER, DIM=size(scene, /DIM), VALUE=1) empty_map = make_array(/FLOAT , DIM=size(scene, /DIM), VALUE=!VALUES.F_NAN) writefits, flat_map_fn, flat_map, header ;; Load observation data into ds9. print, 'Spawning ds9 to perform coordinate conversions ...' ae_send_to_ds9, my_ds9, NAME='geometric_areas_'+session_name ae_send_to_ds9, my_ds9, scene_fn ; Create a data structure to store the list of image indexes belonging to each tesselate. pixels_in_tesselate = ptrarr(num_sources, /ALLOC) for ii=0,num_sources-1 do begin sourcedir = sourcename[ii] + '/' unnamed_src_stats_fn = sourcedir + 'source.stats' ;; Fail if the region file uses the "field()" construct, which will corrupt the algorithm we use in recipe.txt to compute the geometric area of the extraction region. ae_ds9_to_ciao_regionfile, catalog_region_fn[ii], '/dev/null', FIELD_SYNTAX_FOUND=field_syntax_found if field_syntax_found then begin print, 'ERROR: the syntax "field()" is not allowed in region files defining diffuse sources.' retall endif ;; Load region file into ds9 and resave in PHYSICAL coordinates. cmd = strarr(6) cmd[0] = string(my_ds9, F='(%"xpaset -p ''%s'' regions delete all")') cmd[1] = string(my_ds9, F='(%"xpaset -p ''%s'' regions format ds9")') cmd[2] = string(my_ds9, catalog_region_fn[ii], F='(%"xpaset -p ''%s'' regions load %s")') cmd[3] = string(my_ds9, F='(%"xpaset -p ''%s'' regions system physical")') cmd[4] = string(my_ds9, F='(%"xpaset -p ''%s'' regions format ciao")') cmd[5] = string(my_ds9, temp_region_fn, F='(%"xpaset -p ''%s'' regions save %s")') run_command, cmd, /QUIET ; Use CIAO to find the scene pixels that are inside the region. cmd = string(flat_map_fn, temp_region_fn, temp_image_fn, F="(%'dmcopy ""%s[sky=region(%s)][opt full]"" %s clobber+')") run_command, cmd mask = readfits(temp_image_fn) *pixels_in_tesselate[ii] = where(mask GT 0, num_pixels_in_tesselate) geometric_area = num_pixels_in_tesselate * (3600*sxpar(header,'CDELT2'))^2 print, sourcename[ii], geometric_area, F='(%"%30s %0.5g arcsec^2")' ; Write area to source.stats. unnamed_src_stats = headfits(unnamed_src_stats_fn, ERRMSG=error) if keyword_set(error) then message, 'ERROR: could not read '+unnamed_src_stats_fn fxaddpar, unnamed_src_stats, 'SRC_AREA', geometric_area, '[arcsec**2] geometric area of region' writefits, unnamed_src_stats_fn, 0, unnamed_src_stats endfor ; Save the map template and the data structure holding the list of image indexes belonging to each tesselate. save, num_sources, sourcename, empty_map, pixels_in_tesselate, header, FILE='pixels_in_tesselate.sav' ; Clean up temp files. if file_test(tempdir) then begin list = reverse(file_search(tempdir,'*',/MATCH_INITIAL_DOT,COUNT=count)) if (count GT 0) then file_delete, list file_delete, tempdir endif exit end ===================================================================== GROUP THE SPECTRA PLOT THE EXTRACTED AND STOWED BACKGROUND SPECTRA Copy the XSPEC scripts distributed with AE to a local directory. cp -R /usr/common/itt/lib/apps/TARA/code/ae/xspec_scripts . idl .run acis_extract, 'diffuse_all.srclist', /FIT_SPECTRA, CHANNEL_RANGE=[35,480], MODEL_FILENAME='xspec_scripts/plot_gross_and_background_spectra.xcm', SNR_RANGE=[5,10] acis_extract, 'diffuse_all.srclist', /FIT_SPECTRA, CHANNEL_RANGE=[35,480], MODEL_FILENAME='xspec_scripts/plot_gross_and_background_spectra.xcm', SNR_RANGE=[5,20] acis_extract, 'diffuse_all.srclist', /FIT_SPECTRA, CHANNEL_RANGE=[35,480], MODEL_FILENAME='xspec_scripts/plot_gross_and_background_spectra.xcm', SNR_RANGE=[5,30] end plot_name = file_search('*/spectral_models/*plot_gross_and_background_spectra/ldata.ps', COUNT=num_plots) findpro, 'acis_extract', DIRLIST=package_dir cmd = package_dir[0]+'plot_spectra.pl ' + strjoin(plot_name, ' ') spawn, cmd ===================================================================== FIT A VARIETY OF SPECTRAL MODELS TO EACH SOURCE --------------------------------------------------------------------- Devise a Model for the Celestial Background (IF YOU CHOOSE TO USE THIS STRATEGY) In any diffuse analysis it is important to choose a suitable model for the background components not represented in the stowed data (e.g. emission from extragalactic sources, the local hot bubble, and solar wind charge exchange). Some investigators try to physically model each of these components. Instead, we generally define a "celestial background" region within our observation, fit an arbitrary model to that spectrum, and use that frozen celestial background model in subsequent fits of our diffuse sources. If you have several candidate regions you're considering for the role of celestial background, then it is convenient to extract each of them separately, and then try to build a celestial background model in XSPEC that fits all of them simultaneously. During that process, plotting the fit residuals should help you decide which of the candidate regions have spectra that are similar enough to be adopted as the celestial background region. As you build this model, make sure that you are not using any XSPEC configurations that are different from those used in the AE fitting scripts, for example the setting of the "abund" command. You can manually apply AE's grouping algorithm to the spectra you have extracted so you can work with them in XSPEC: idl acis_extract, 'diffuse_all.srclist', /FIT_SPECTRA, CHANNEL_RANGE=[35,480], MODEL_FILENAME='xspec_scripts/nofit.xcm' The grouped spectrum produced above can be directly loaded into XSPEC and fit with whatever model you dream up. We normally freeze all parameters of this celestial background model, then save the model out of XSPEC, e.g. save model background.xcm Then edit the .xcm file you just created as follows: 1. As documentation for how the model was derived comment out (with the # character) commands that configure physical models in XSPEC, including: xsect abund 2. Remove all commands not required to define the model, including lines with these XSPEC commands: statistic method cosmo xset systematic 3. Specify a suitable name for this celestial background model by adding "3:sky" to the "model" command, e.g. model 3:sky TBabs(apec + apec + apec) + powerlaw Look in the AE fitting script(s) you will be using, just before the section where the source model is defined, for an optional section that shows how to load this celestial background model; remove the comment character from the 3 lines of code in that section, changing the name of your model .xcm file as required. --------------------------------------------------------------------- Configure Fitting Scripts Revise the plausible ranges for model parameters found in whatever XSPEC scripts you will be using to be appropriate for your sources. The *_min and *_max variables in the fitting scripts are used as "soft parameter limits" in XSPEC, and are used to identify "limit violations" in our table generator. NH in these fitting scripts is in units of 1X10^22 cm^(-2). Since the instrumental background includes an emission line above 7 keV, we usually modify the point source fitting process in two ways: 1. In the FIT_SPECTRA call to AE, we specify a fitting band of 0.5:7.0 keV (channels 35:480]. 2. We edit the fitting scripts (variables EbandLo, EbandHi, and BandName) to specify that full-band and hard-band fluxes are calculated on bands that stop at 7 keV instead of 8 keV. --------------------------------------------------------------------- Perform Fitting Runs An EXAMPLE AE fitting and collation run is shown below. The fitting script tbabs_2vapec+bkg.xcm is a copy of AE's 2-temperature script that has been modified to use flux bands that stop at 7 keV and to load a celestial background model. The "changes" script solar.xcm configures the vapec source models to use solar abundances. nice idl |& tee -a fit_tbabs_2vapec+bkg.log .run acis_extract, 'diffuse_all.srclist', /FIT_SPECTRA, CHANNEL_RANGE=[35,480], MODEL_FILENAME='tbabs_2vapec+bkg.xcm', MODEL_CHANGES_FILENAME='xspec_scripts/solar.xcm' acis_extract, 'diffuse_all.srclist', COLLATED_FILENAME='tables/tbabs_vapec+bkg_solar.collated', HDUNAME='*tbabs_vapec+bkg_solar*' exit end Examine the AE messages found in the log files to find failures other than an XSPEC timeout, which is handled explictly later. egrep -i "WARNING|ERROR|halted" fit*.log |egrep -v "DISPLAY variable|arithmetic error|unignored|skip_errors|variable_name|status 152|process killed after" Refit with parameter error estimation disabled any runs where XSPEC timed out (as shown in the point source section of this recipe): awk -F "[ /]+" '/process killed after/ {print $6}' fit_tbabs_2vapec+bkg.log > noerr.srclist wc -l noerr.srclist You may need to hand-fit some spectra (using AE's /INTERACT option), change various parameter soft limits, thaw various abundances, try other models, etc. When you're finished fitting you need to collate the preferred model for each source. If your workflow was such that the _last_ fit run on each source is the one you like, then just collate like this: idl acis_extract, 'diffuse_all.srclist', COLLATED_FILENAME='tables/best_model.collated' Otherwise, you'll have to explicitly declare which model is preferred for each source, using AE's "BEST_MDL" mechanishm (a FITS keyword BEST_MDL placed into the source.spectra files). This can be done interactively with the ae_spectra_viewer tool, or code you write can do it with the ae_default_model_preference tool. Once a BEST_MDL has been recorded for every source, then collate like this: acis_extract, 'diffuse_all.srclist', COLLATED_FILENAME='tables/best_model.collated', HDUNAME='BEST_MDL' ===================================================================== CREATE A "FACEBOOK" OF THE XSPEC PLOTS FROM THE DIFFUSE REGIONS The code below calls a Perl script, plot_spectral.pl, that builds a Latex document (spectra.tex) that presents the XSPEC plots 12 to a page. Latex is run, and the PostScript document spectra.ps is produced. idl bt = mrdfits('tables/best_model.collated',1) findpro, 'acis_extract', DIRLIST=package_dir plot_name = strtrim(bt.CATALOG_NAME,2) + '/spectral_models/' + strtrim(bt.MODEL,2) + '/ldata.ps' ind = where(~file_test(plot_name), count) if (count GT 0) then plot_name[ind] = package_dir[0]+'no_model.eps' cmd = package_dir[0]+'plot_spectra.pl ' + strjoin(plot_name, ' ') spawn, cmd ===================================================================== COMPUTE LUMINOSITIES, EXAMINE DISTRIBUTIONS OF SPECTRAL PARAMETERS At this point it is helpful to compute luminosities from XSPEC fluxes, and to display distributions of spectral parameters. This is easily done by the tool ae_flatten_collation. Examination of the resulting plots and of the messages produced by the tool may suggest spectral models that are suspicious and should be examined more closely. The ae_flatten_collation tool will also identify fit parameters which violate the stated parameter ranges, omitting those from the FITS and ASCII tables is produces. In many such cases you may decide to refit with the offending parameter frozen at its limit. The tool tries to report various problems in the spectral models, for example error flags that XSPEC reported. You may iterate a few times re-fitting problem sources, collating, and running ae_flatten_collation. In the code below, specify the distance to your target. cd tables idl |& tee -a flatten.log ; Supply distance in pc. distance = ?? ; pc .run ae ae_flatten_collation, COLLATED_FILENAME='best_model.collated', FLAT_TABLE_FILENAME='xray_properties.fits', /DIFFUSE, DISTANCE=distance, /FAST, SORT=0 This tool will produce both FITS and ASCII tables. It may be convenient for Co-I's to format the ASCII table in PDF. I cannot get Apple's TextEdit to do this, but TextWrangler can: * Use File:Page Setup to define a "wide" page (perhaps 40" x 12"). * Use File:Print to bring up print dialog box, where you can "Save as PDF". Keep in mind that for diffuse sources AE scales the calibration (ARF) so that the spectral model derived by XSPEC is on a "per arcsec^2" basis. See Section 5.12 in the AE manual. All ``flux'' quantities computed by XSPEC should then be understood to be surface brightness quantities with units of (erg /s /cm**2 /arcsec**2). The ae_flatten_collation tool uses the distance to the target to produce "FitLuminosity" quantities in units of log(erg /s /pc**2). The ae_flatten_collation tool produces histograms of various luminosities. If you'd like to see the distribution of some combination of luminosities, you'll have to do the work yourself, for example: idl flat_table_fn = 'xray_properties.fits' bt = mrdfits(flat_table_fn, 1) dataset_1d, id1, alog10(10.^(bt.FitLuminosity_1t) + 10.^(bt.FitLuminosity_2t) + 10.^bt.FitLuminosity_3t ), DATASET='SB1+2+3H7', BINSIZE=0.2, XTIT='log(erg /s /pc**2)' dataset_1d, id1, alog10(10.^(bt.FitLuminosity_1tc) + 10.^(bt.FitLuminosity_2tc) + 10.^bt.FitLuminosity_3tc), DATASET='SB1+2+3CH7', BINSIZE=0.2 ===================================================================== CREATE LATEX TABLES Edit hmsfr_tables.tex (used below) so that the column definitions for the diffuse table correspond to the set of elements that are represented in the Abundance_CI column of xray_properties.fits. idl |& tee -a hmsfr.log .r ae flat_table_fn = 'xray_properties.fits' ; Supply distance in pc. distance = ?? nominal_effective_area = 1 ; Not applicable to diffuse sources. ; Supply local path to table templates. findpro, 'acis_extract', DIRLIST=package_dir template_filename = package_dir+'/hmsfr_tables.tex' hmsfr_tables, flat_table_fn, template_filename, nominal_effective_area, /DIFFUSE, DISTANCE=distance, THERMAL_MODEL_PATTERN='*diffuse*' Edit xray_properties.tex to remove empty tables. Remove abundance entries that are frozen to 1.0. 'mv' diffuse_spectroscopy.tex temp sed -e 's:\$1.0\*\\phd\$: :g' temp > diffuse_spectroscopy.tex pdflatex xray_properties; open xray_properties.pdf OR pdflatex xray_properties; gv --orientation=landscape --scale=2 xray_properties.pdf & ===================================================================== CREATE MAPS FOR INTERESTING PARAMETERS OF THE SPECTRAL MODELS. Create maps of interesting quantities. Adjust the code below as appropriate for the model you have fit. The file pixels_in_tesselate.sav in the main diffuse directory (created earlier in the section CALCULATE THE GEOMETRIC AREAS OF THE DIFFUSE REGIONS) stores the indices for image pixels enclosed by each tesselate. mkdir maps cd maps idl |& tee -a build_maps.log .r ae restore, /V, '../pixels_in_tesselate.sav' ; Specify the location of the collation holding fit parameters we want to map. flat_table_fn = '../tables/xray_properties.fits' bt = mrdfits(flat_table_fn, 1, theader) num_sources = n_elements(bt) NAME = strtrim(bt.NAME,2) LABEL = strtrim(bt.LABEL ,2) restore, /V, '../pixels_in_tesselate.sav' .run if (total(/INT, NAME NE sourcename) GT 0) then message, 'ERROR: sources in collation do not correspond to sources in pixels_in_tesselate.sav' ; --------------------------------------------------------------- ; Build maps of sb. chi_map = empty_map sb_soft_observed_map = empty_map sb_observed_map = empty_map sb1_observed_map = empty_map sb2_observed_map = empty_map sb3_observed_map = empty_map sb4_observed_map = empty_map sb5_observed_map = empty_map sb0_corrected_map = empty_map sb1_corrected_map = empty_map sb2_corrected_map = empty_map sb3_corrected_map = empty_map sb4_corrected_map = empty_map sb5_corrected_map = empty_map kt_avg12_map = empty_map for ii=0L,num_sources-1 do begin ind = *pixels_in_tesselate[ii] chi_map [ind] = (bt.ReducedChiSq) [ii] sb_soft_observed_map [ind] = (bt.FitLuminosity_s) [ii] sb_observed_map [ind] = (bt.FitLuminosity_t) [ii] sb1_observed_map [ind] = (bt.FitLuminosity_1t) [ii] sb2_observed_map [ind] = (bt.FitLuminosity_2t) [ii] sb3_observed_map [ind] = (bt.FitLuminosity_3t) [ii] sb4_observed_map [ind] = (bt.FitLuminosity_4t) [ii] sb5_observed_map [ind] = (bt.FitLuminosity_5t) [ii] sb1_corrected_map[ind] = (bt.FitLuminosity_1tc)[ii] sb2_corrected_map[ind] = (bt.FitLuminosity_2tc)[ii] sb3_corrected_map[ind] = (bt.FitLuminosity_3tc)[ii] sb4_corrected_map[ind] = (bt.FitLuminosity_4tc)[ii] sb5_corrected_map[ind] = (bt.FitLuminosity_5tc)[ii] ; kt_avg12_map [ind] = ((bt.NORM1*bt.KT1 + bt.NORM2*bt.KT2) / (bt.NORM1 + bt.NORM2))[ii] endfor ; Add some FITS keywords to map headers and write the maps. fxaddpar, header, 'HDUNAME', 'reduced chi^2 map' writefits, 'chi_map.img', chi_map, header fxaddpar, header, 'BUNIT', 'log(erg /s /pc**2)' fxaddpar, header, 'HDUNAME', 'SBH2 map' writefits, 'SBH2_map.img', sb_soft_observed_map, header fxaddpar, header, 'HDUNAME', 'SBH7 map' writefits, 'SBH7_map.img', sb_observed_map, header fxaddpar, header, 'HDUNAME', 'SB1H7 map' writefits, 'SB1H7_map.img', sb1_observed_map, header fxaddpar, header, 'HDUNAME', 'SB2H7 map' writefits, 'SB2H7_map.img', sb2_observed_map, header fxaddpar, header, 'HDUNAME', 'SB3H7 map' writefits, 'SB3H7_map.img', sb3_observed_map, header fxaddpar, header, 'HDUNAME', 'SB1+2+3H7 map' writefits, 'SB1+2+3H7_map.img', alog10(10.^sb1_observed_map+10.^sb2_observed_map+10.^sb3_observed_map), header fxaddpar, header, 'HDUNAME', 'SB4H7 map' writefits, 'SB4H7_map.img', sb4_observed_map, header fxaddpar, header, 'HDUNAME', 'SB5H7 map' writefits, 'SB5H7_map.img', sb5_observed_map, header ; fxaddpar, header, 'HDUNAME', 'SB0CH7 map' ; writefits, 'SB0CH7_map.img', sb0_corrected_map, header fxaddpar, header, 'HDUNAME', 'SB1CH7 map' writefits, 'SB1CH7_map.img', sb1_corrected_map, header fxaddpar, header, 'HDUNAME', 'SB2CH7 map' writefits, 'SB2CH7_map.img', sb2_corrected_map, header fxaddpar, header, 'HDUNAME', 'SB1+2+3CH7 map' writefits, 'SB1+2+3CH7_map.img', alog10(10.^sb1_corrected_map+10.^sb2_corrected_map+10.^sb3_corrected_map), header ; fxaddpar, header, 'HDUNAME', 'KT_avg12_map' ; writefits, 'KT_avg12_map.img', kt_avg12_map, header fxaddpar, header, 'HDUNAME', 'SB3CH7 map' writefits, 'SB3CH7_map.img', sb3_corrected_map, header fxaddpar, header, 'HDUNAME', 'SB4CH7 map' writefits, 'SB4CH7_map.img', sb4_corrected_map, header fxaddpar, header, 'HDUNAME', 'SB5CH7 map' writefits, 'SB5CH7_map.img', sb5_corrected_map, header end .run ; --------------------------------------------------------------- ; Loop over each fit parameter, making both a grayscale map showing the best-fit values ; and a color-code map that tries to depict both best-fit and confidence interval information. ; --------------------------------------------------------------- param_names = ['NH1','NH2','NH3','KT1','KT2','KT3','TAU1','TAU2','TAU3','O','Neon','MG','S','SI','FE', 'NH4','NH5','KT4','KT5'] for jj=0,n_elements(param_names)-1 do begin if ~tag_exist(bt, param_names[jj], index=colnum) then begin print, 'WARNING! Could not find parameter '+param_names[jj] continue endif cmd = ['param_value','upper_confidence_limit','lower_confidence_limit','param_range'] + $ ' = bt.' + param_names[jj] + $ ['','_HiLim','_LoLim','_CI'] for kk=0,3 do begin if ~execute(cmd[kk]) then message, 'Cannot execute: '+cmd[kk] endfor param_was_frozen = strmatch(param_range, '*\**') tunit_kywd = string(1+colnum,F='(%"TUNIT%d")') tunit = sxpar(theader, tunit_kywd) ;; ================================================= ;; Encode best-fit and confidence interval information using color. ;; ================================================= min_hue = 0 max_hue = 260 legend_xsize = (size(empty_map,/DIM))[0] / 3 legend_ysize = 60 saturation_floor = 0.3 value_floor = 0.3 ;---------------------------------------------------- ; Strategy I for using color to depict confidence interval information ; ; We choose to present the best-fit values as a percentage offset from a typical value, e.g. the median value, ; so that abundances of different elements can be more appropriately compared. reference_value = median(param_value) param_percentage_offset = 100 * (param_value - reference_value) / reference_value ; Encode these parameter value ratios using the "hue" axis of the HSV color model. ; We have to pick some scaling: ; param_percentage_offset = -50% will be red (hue=min_hue) and ; param_percentage_offset = +50% will be purple (hue=max_hue). hue = min_hue > ( min_hue + max_hue * (param_percentage_offset - (-50)) / (+50 - (-50)) ) < max_hue print, param_names[jj], reference_value, reference_value, min_hue, max_hue, F='(%"\nMapping the fractional offset from a typical value, (%s - %0.2f)/%0.2f, with HSV color model.\n Hue=%d depicts a fractional offset of -50%%, Hue=%d depicts a fractional offset of +50%%")' ; Just as parameter values were normalized to be unitless quantities, we choose to present confidence interval information ; as the average percentage offset from the best-fit parameter value. upperlimit_percentage_offset = 100 * (upper_confidence_limit - param_value) / param_value lowerlimit_percentage_offset = 100 * (lower_confidence_limit - param_value) / param_value ; For display purposes, we average the upper and lower uncertainty sizes. legend_limit_percentage_offset = [-90,-80,-70,-60,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,70,80,90] limit_percentage_offset = replicate(!values.f_nan, num_sources) for ii=0L,num_sources-1 do $ limit_percentage_offset[ii] = mean(/NAN, abs([upperlimit_percentage_offset[ii], lowerlimit_percentage_offset[ii]])) ; When both upper and lower intervals are not known, assume they are large. ind = where(~finite(limit_percentage_offset), count) if (count GT 0) then limit_percentage_offset[ind] = 200 ; In order to encode this limit_percentage_offset via the saturation or value axes of the HSV model, we must ; convert to metric that is 1.0 for "good" tiles (tight confidence intervals) and is smaller for large ones. ; Let's arbitrarily scale so that ; when limit_percentage_offset= 0 (zero errors) we display a saturation/value of 1.0 ; when limit_percentage_offset=50 (50% errors) we display a saturation/value of 0.5 uncertainty_metric = 0 > (1.0 - (0.5/50)*abs(limit_percentage_offset) ) legend_uncertainty_metric = 0 > (1.0 - (0.5/50)*abs(legend_limit_percentage_offset) ) print, F='(%" ''Value'' coordinate encodes average percentage offset from the best-fit parameter value to the upper/lower limits:\n %%offset Value\n ---------------")' forprint, legend_limit_percentage_offset, legend_uncertainty_metric > value_floor, F='(%" %4d %4.2f")' ; Encode the uncertainty_metric using either the "saturation" or "value" axes of the Hue-Saturation-Value color model. saturation = replicate(1.0, num_sources) value = uncertainty_metric ; Build a corresponding colorbar to serve as a legend. hue_ramp = min_hue + (max_hue-min_hue)*(findgen(legend_xsize)/(legend_xsize-1)) value_ramp = value_floor + (1.0-value_floor)*(findgen(legend_ysize)/(legend_ysize-1)) make_2d, hue_ramp, intarr(legend_ysize), legend_h, junk make_2d, intarr(legend_xsize), value_ramp, junk , legend_v legend_s = replicate(1.0, legend_xsize, legend_ysize) ;---------------------------------------------------- ; Build the grayscale and color-coded maps. param_map = empty_map hmap = empty_map smap = empty_map vmap = empty_map for ii=0L,num_sources-1 do begin ind = *pixels_in_tesselate[ii] ; MARK TILES THAT ARE FROZEN by leaving some pixels set to NaN if param_was_frozen[ii] then begin index_to_point, ind, xind, yind, size(empty_map) xind -= median(xind) yind -= median(yind) good = where(xind NE -yind) ind = ind[good] endif param_map[ind] = param_value[ii] hmap [ind] = hue[ii] smap [ind] = saturation[ii] vmap [ind] = value[ii] endfor ;ii ;print, 'Writing maps for ', param_names[jj] fxaddpar, header, 'HDUNAME', param_names[jj]+' map' fxaddpar, header, 'BUNIT', tunit writefits, param_names[jj]+'_map.img', param_map, header ; Convert to RGB color system, adopting a floor on saturation and value so that we retain some color in all tiles. color_convert, legend_h, saturation_floor > legend_s < 1.0, value_floor > legend_v < 1.0, legend_r, legend_g, legend_b, /HSV_RGB color_convert, hmap, saturation_floor > smap < 1.0, value_floor > vmap < 1.0, rmap, gmap, bmap, /HSV_RGB ; Set the NaN pixels to a specific color ind = where(~finite(hmap), count) if (count GT 0) then begin rmap[ind] = 255 gmap[ind] = 255 bmap[ind] = 255 endif ; Paste in the color bar legend. x0 = (size(empty_map,/DIM))[0] - legend_xsize y0 = 0 x0 = 126 y0 = 923 rmap[x0,y0] = legend_r gmap[x0,y0] = legend_g bmap[x0,y0] = legend_b writefits, param_names[jj]+'_rmap.img', rmap, header writefits, param_names[jj]+'_gmap.img', gmap, header writefits, param_names[jj]+'_bmap.img', bmap, header endfor ;jj end View the greyscale maps in ds9 with linear, minmax scaling ds9 -scale mode minmax -linear SB*_map.img & ds9 -scale mode minmax -linear NH*_map.img & ds9 -scale mode minmax -linear KT*_map.img& ds9 -scale mode minmax -linear O_map.img Neon_map.img MG_map.img S_map.img SI_map.img FE_map.img & You MUST view these 3-color maps in ds9 with linear, minmax scaling to produce the intended colors!! You must NOT RESCALE these RGB images in ds9! foreach name (NH1 NH2 NH3 KT1 KT2 KT3 TAU1 TAU2 TAU3 O Neon Mg S SI FE) ds9 -scale mode minmax -linear -rgb -rgb lock colorbar yes -rgb lock scale yes -red ${name}_rmap.img -green ${name}_gmap.img -blue ${name}_bmap.img & end ; ;---------------------------------------------------- ; ; Strategy II for using color to depict confidence interval information ; ; ; ; We choose to normalize the best-fit values by a typical value, e.g. the "bisector" computed by ; ; validate_xspec_confidence_interval, so that abundances of different elements can be more appropriately compared. ; ; Log scaling seems appropriate, and it's convenient to use log base 2 so that a ratio of 0.5 maps to -1 ; ; and a ratio of 2.0 maps to +1. ; param_percentage_offset = alog(param.param_value / param.bisector) / alog(2.0) ; ; ; Encode these parameter value ratios using the "hue" axis of the HSV color model. ; ; We have to pick some scaling: ; ; a ratio of 0.5 (param_percentage_offset==-1) will be red (hue=0) and ; ; a ratio of 2.0 (param_percentage_offset==+1) will be purple (hue=max_hue). ; ; hue = min_hue > ( min_hue + max_hue * (param_percentage_offset + 1.0) / 2.0 ) < max_hue ; ; ; ; ; Compute the significance of each tile's deviation from a canonical "bisector" value, based on where the bisector value lies with respect to the confidence interval. ; ; Encode this deviation significance using either the "saturation" or "value" axes of the Hue-Saturation-Value color model. ; ; ; ; Tiles with typical parameter values (near the bisector) will be greenish, and almost always have low saturation or value since the bisector will be well within their confidence intervals. ; ;---------------------------------------------------- ; ; bisector_is_excluded = (param.lower_confidence_limit GE param.bisector) OR (param.upper_confidence_limit LE param.bisector) ; ; print, param_names[jj], F='(%"\nMapping %s ...")' ; print, total(/INT, bisector_is_excluded), param.bisector, F='(%"These %d sources have confidence intervals that exclude the bisector value %0.2f ")' ; ; ; We have to pick some scaling for the confidence interval information---we choose to map the location of the bisector value within the parameter's confidence interval and to scale such that this value is 1.0 when the bisector is outside the interval. ; ; deviation_significance = fltarr(num_sources) ; ; ; Deal with cases where param_value is less than the bisector. ; location = ((param.bisector - param.param_value) / (param.upper_confidence_limit - param.param_value)) ; ; ind = where(location GT 0, count) ; if (count GT 0) then deviation_significance[ind] = location[ind] ; ; ; Deal with cases where param_value is greater than the bisector. ; location = ((param.bisector - param.param_value) / (param.lower_confidence_limit - param.param_value)) ; ; ind = where(location GT 0, count) ; if (count GT 0) then deviation_significance[ind] = location[ind] ; ; ; Encode the deviation_significance using either the "saturation" or "value" axes of the Hue-Saturation-Value color model. ; saturation = replicate(1.0, num_sources) ; value = deviation_significance #################### APPENDIX A: STOWED BACKGROUND DATA #################### The "stowed background" data in CALDB is discussed in these documents: Hickox, R.~C., \& Markevitch, M.\ 2006, \apj, 645, 95 http://cxc.harvard.edu/caldb/calibration/acis.html (Blank Field Event Files section) http://cxc.harvard.edu/cal/Links/Acis/acis/Cal_prods/bkgrnd/current/ http://cxc.harvard.edu/cal/Acis/Cal_prods/bkgrnd/acisbg/index.html http://cxc.harvard.edu/cal/Acis/Cal_prods/bkgrnd/acisbg/data/README The key recommendations for how to perform equivalent processing on the stowed and science data are found in Maxim's cookbook: http://cxc.harvard.edu/contrib/maxim/acisbg/COOKBOOK In this recipe we use the stowed data that was CTI corrected by the CXC (since we've now retired our corrector). --------------------------------------------------------------------- FIND THE APPROPRIATE STOWED DATA (for EACH observation) Build a list of the CALDB stowed data corresponding to the CCDs used in your observation. The epoch of CALDB files you want depends on when your observation was taken! See description of stowed data epochs "D" and "E" at http://cxc.harvard.edu/contrib/maxim/acisbg/data/README Both epochs are combined for FI devices; one epoch is selected for BI devices using the observation date. set ccdlist=` dmkeypar acis.astrometric.calibrated.subpix.evt1 DETNAM echo=yes | sed -e 's/ACIS-//'` set date_obs=`dmkeypar acis.astrometric.calibrated.subpix.evt1 DATE-OBS echo=yes | cut -c1-4,6,7,9,10` printf "# ACIS-${ccdlist}\n" > CALDB_stowed.txt printf "# DATE-OBS = $date_obs \n" >> CALDB_stowed.txt 'ls' -1 $CALDB/data/chandra/acis/bkgrnd/acis[01234689]D200[05]*bgstow_cti* | egrep ".+acis[${ccdlist}].+" >> CALDB_stowed.txt if ($date_obs > 20050901) then 'ls' -1 $CALDB/data/chandra/acis/bkgrnd/acis[57]D2005*bgstow_cti* | egrep ".+acis[${ccdlist}].+" >> CALDB_stowed.txt else 'ls' -1 $CALDB/data/chandra/acis/bkgrnd/acis[57]D2000*bgstow_cti* | egrep ".+acis[${ccdlist}].+" >> CALDB_stowed.txt endif cat CALDB_stowed.txt --------------------------------------------------------------------- MERGE THE STOWED DATA FILES (for EACH observation). AS MUCH AS POSSIBLE, APPLY THE SAME CLEANING STEPS USED ON THE DIFFUSE DATA. We don't have the courage to worry about the fact that the stowed background and our observation used different levels of background flare cleaning. We don't have the courage to worry about the differences between the bad pixel table used by Maxim and the table we are using on our data. punlearn dmmerge dmcopy dmmerge infile="@CALDB_stowed.txt" outfile=acis.stowed.calibrated.evt1 clobber=yes IF THE DATAMODE of your observation is NOT VFAINT, then you must CLEAR the Clean55 bits in the stowed data so that the observation and stowed data sets will have equivalent filtering applied (status=0 filter below)!!! if (`dmkeypar acis.astrometric.calibrated.subpix.evt1 DATAMODE echo=yes` != 'VFAINT') then echo "CLEARING THE CLEAN55 BIT IN THE STOWED DATA" mv -f acis.stowed.calibrated.evt1 temp.evt dmtcalc infile=temp.evt outfile=acis.stowed.calibrated.evt1 expression="status=status,status=X8F" clobber=yes endif Find out what energy filter has been applied to your diffuse event data, and change the energy filter below to match! dmstat "diffuse.evt[cols energy]" dmcopy "acis.stowed.calibrated.evt1[grade=0,2,3,4,6,energy<12000,status=0]" stowed.evt --------------------------------------------------------------------- REPROJECT THE STOWED DATA (for EACH observation) We must add a fake TIME column to the stowed data with random times from the real observation so that reproject_events can compute appropriate sky coordinates for the stowed events. Randomization of event positions (controlled by the parameter "random") should not be important for background data, which have no point sources. We choose to turn it off anyway (random=-1 in call to reproject_events). cd obsXXXX/ punlearn dmhedit dmtcalc reproject_events # Copy range of observation's first GTI table to TSTART and TSTOP in the stowed data. dmstat "diffuse.evt[3][cols START]" set tstart=`pget dmstat out_min` dmstat "diffuse.evt[3][cols STOP]" set tstop=`pget dmstat out_max` dmhedit "stowed.evt[2]" filelist=none operation=add key=TSTART value=$tstart dmhedit "stowed.evt[2]" filelist=none operation=add key=TSTOP value=$tstop # Add to the stowed data a random TIME column with values on [TSTART:TSTOP]. dmtcalc stowed.evt temp.evt expression="TIME=TSTART+((TSTOP-TSTART)*(#rand))" clobber=yes # Reproject the stowed data to the tangent plane of the observation. reproject_events infile=temp.evt aspect=obs.asol match=diffuse.evt outfile=stowed.evt random=-1 clobber=yes --------------------------------------------------------------------- MAKE "BACKGROUND EXPOSURE MAPS" FOR THE STOWED DATA (for EACH observation) When extracting diffuse sources, we will want to model the instrumental background by extracting the "stowed" ACIS data falling within each source aperture. Those "stowed" spectra will have to be appropriately scaled to match the "depth" of the ACIS observation within that aperture. A background spectrum in AE is scaled by the ratio of two integrals---the integral of a "source" exposure map over the source aperture, and the integral of a "background" exposure map over the background aperture. Recall that a "source" exposure map conflates three separate effects contributing to the "depth" of the observation at each position on the sky: 1. The effective area of the HRMA and ACIS is non-uniform across the ACIS detector. 2. Since the pointing of Chandra is dithered, the amount of time a point on the sky fell onto the detector varies strongly with position. 3. Since each CCD's telemetry stream is independent, each CCD can have a different total observing time. Thus, we must construct a similar "background" exposure map that, when integrated over the same aperture, produces the appropriate scaling of the stowed spectrum. Since effect #1 (effective area variation) is unrelated to the instrumental background, we would like it to appear in the background exposure map unaltered. Since effect #2 (dither effects) is the same for both X-rays and instrumental background we would again like it to appear in the background exposure map unaltered. In contrast, effect #3 (total observing time on each CCD) is clearly very different between the science observation and the calibration observations used to build the stowed dataset in CALDB. Thus, the "background" exposure map that we desire should be constructed by appropriately scaling and then summing the single-CCD exposure maps for the observation. Approximate scalings for the single-ObsID exposure maps could be computed from the single-CCD EXPOSUR? values found in the observation and stowed event lists. However, the ACIS instrumental background is known to be time-varying, and thus the proper scaling cannot be known only from observing times. Hickox & Markevitch (2006, \apj, 645, 95, Section 4.1.3) have shown that good results can be obtained by scaling a stowed spectrum so that it matches the observed spectrum over the energy range 9-12 keV, where there are no X-rays. Thus, the code below compares single-ObsID 9-12 keV spectra from the observation and stowed data to derive an adjustment for each single-CCD exposure map. IT IS VITAL THAT YOU SCALE THE STOWED DATA WITH RESPECT TO WHATEVER OBSERVATION DATA YOU WILL BE USING FOR DIFFUSE ANALYSIS!! 1. The two event lists must cover the SAME field of view. For example, you must not compare a masked observation event list to an unmasked stowed event list. 2. The two event lists must be similarly cleaned. For example, on VFAINT data, an event list that was not cleaned by the acis_detect_afterglow and Clean55 cleaning algorithms (spectral.evt in our recipes) may have a count rate in the 9-12 keV band that is 1.7 times that obtained when those algorithms are applied (validation.evt in our recipes). In the code below, obsdata_filename should be the name of an UNMASKED event list that has cleaning similar to that applied to diffuse.evt. nice idl |& tee -a ae_make_emap.log .run ae .run obsdata_filename = 'validation.evt' stowed_filename = 'stowed.evt' bt_obs = mrdfits(obsdata_filename, 1, obs_header) bt_stowed = mrdfits(stowed_filename , 1, stowed_header) ; Scaling of stowed data uses the energy range 9-12 keV (Hickox & Markevitch, 2006, \apj, 645, 95, Section 4.1.3). emax = max(bt_obs.energy) < max(bt_stowed.energy) < 12000 emin = 9000 obs_ccd_hist = histogram(MIN=0,MAX=9,BIN=1, ( bt_obs.ccd_id)[where(( bt_obs.energy GE emin) AND ( bt_obs.energy LE emax))]) stowed_ccd_hist = histogram(MIN=0,MAX=9,BIN=1, (bt_stowed.ccd_id)[where((bt_stowed.energy GE emin) AND (bt_stowed.energy LE emax))]) ccd_scaling = stowed_ccd_hist / float(obs_ccd_hist) ; Plot stowed and observed spectra in the scaling energy band. iarray_scaling = mean(ccd_scaling[0:3]) dataset_1d, id, ( bt_obs.energy)[where( bt_obs.ccd_id LE 3)], DATASET='I-array observed', XTITLE='energy (eV)' dataset_1d, id, (bt_stowed.energy)[where(bt_stowed.ccd_id LE 3)], DATASET='I-array stowed (scaled)', NORM_ABSC=[0,15000], NORM_VAL=[iarray_scaling,iarray_scaling],COLOR='red',LINE=0 for ii=0,9 do begin if ~finite(ccd_scaling[ii]) then continue ; Record the effective LIVETIMEs of the stowed data for each CCD. keyname = string(ii, F="(%'LIVTIME%d')") livtime = sxpar( obs_header, keyname ) cmd = string(stowed_filename, keyname, ccd_scaling[ii] * livtime, F="(%'dmhedit %s filelist=none operation=add key=%s value=%0.1f')") print, cmd spawn, cmd ; Display a map of the mismatch significance between the stowed and observed data, for each CCD, to look for patterns. observ_image_spec = string(obsdata_filename, ii, emin, emax, F="(%'%s[ccd_id=%d,energy=%d:%d][bin chip=::128][opt type=r4]')") stowed_image_spec = string( stowed_filename, ii, emin, emax, F="(%'%s[ccd_id=%d,energy=%d:%d][bin chip=::128][opt type=r4]')") outfile = string(ii, F="(%'stowed_error_significance_ccd%d.img')") cmd = string(observ_image_spec+","+stowed_image_spec, outfile, ccd_scaling[ii], (ccd_scaling[ii])^2, $ F="(%'dmimgcalc infile=""%s"" infile2=none outfile=%s operation=""imgout=(img1-(img2/%0.3f))/sqrt(img1+(img2/%0.3f))"" verbose=1 clob+')") print, cmd spawn, cmd spawn, "ds9 -linear "+outfile+" -zoom to fit &" endfor print, emin/1000., emax/1000., obsdata_filename, F="(%'Stowed data is scaled as shown below so the energy range %0.1f:%0.1f keV matches %s:')" forprint, indgen(10), ccd_scaling, F="(%' CCD%d %0.2f')" ; Use the existing observation's emap (obs.emap) as a scene template to build single-CCD emaps, scale those by the values in ccd_scaling, and then sum. ae_make_emap, obsdata_filename, ['trash'], CCD_LIST=['012367'], ARDLIB_FILENAME='ardlib.par', /REUSE_ASPHIST, /REUSE_INSTMAP, MATCHFILE='obs.emap', ASPECT_FN='obs.asol', CCD_SCALING=ccd_scaling, SCALED_PREFIX='stowed.' end mv stowed.trash.emap stowed.emap ds9 -bin factor 8 -log diffuse.evt diffuse.emap stowed.evt stowed.emap -zoom to fit -match frames wcs &