Download the notebook.
CIAO Workshop Exercises¶
This Jupyter notebook uses the bash_kernel. Users will need to have this kernel installed before they can run the commadns in this notebook.
CIAO 4.17 users with pip3 can use the following commands to install bash_kernel into their CIAO installtion:
pip3 install bash_kernel
python -m bash_kernel.install --sys-prefix
jupyter kernelspec list
Introduction¶
These exercises are designed to familiarize students with Chandra data and some of the basics of X-ray data analysis using CIAO with Chandra data.
Students will use a variety of tools and scripts to perform typical Chandra data analysis tasks. These exercises are not intended to be a Cookbook for all Chandra data analysis, but rather are meant to expose students to basic X-ray analysis techniques.
Students should already have CIAO installed on their systems. These exercises have been developed with CIAO 4.17 using CALDB 4.12.0.
In [1]:
### conda activate ciao
### PATH="${ASCDS_INSTALL}/bin:${PATH}"
ciaover -v
The current environment is configured for: CIAO : CIAO 4.17.0 Friday, December 06, 2024 Contrib : Package release 1 Tuesday, April 29, 2025 bindir : /export/miniforge/envs/ciao-4.17/bin Python path : /export/miniforge/envs/ciao-4.17/bin CALDB : 4.12.0 System information: Linux lenin.cfa.harvard.edu 4.18.0-425.3.1.el8.x86_64 #1 SMP Fri Sep 30 11:45:06 EDT 2022 x86_64 x86_64 x86_64 GNU/Linux
Getting to know Chandra data¶
In this section students will obtain Chandra data and learn some techniques to inspect the quality of the data and how to reprocess their dataset.
Download Dataset¶
Exercise 1¶
Obtain the standard data distribution for Chandra OBS_ID 13858.
Q: How did you obtain the data (Chaser, download_chandra_obsid, other)?
In [2]:
/bin/rm -rf 13858
download_chandra_obsid 13858
Downloading files for ObsId 13858, total size is 66 Mb.
Type Format Size 0........H.........1 Download Time Average Rate
---------------------------------------------------------------------------
vvref pdf 31 Mb #################### < 1 s 59647.0 kb/s
evt1 fits 22 Mb #################### < 1 s 74503.7 kb/s
evt2 fits 5 Mb #################### < 1 s 45131.8 kb/s
asol fits 3 Mb #################### < 1 s 58422.8 kb/s
mtl fits 637 Kb #################### < 1 s 30261.0 kb/s
bias fits 495 Kb #################### < 1 s 25360.9 kb/s
bias fits 451 Kb #################### < 1 s 21750.1 kb/s
bias fits 437 Kb #################### < 1 s 27102.3 kb/s
bias fits 423 Kb #################### < 1 s 14718.9 kb/s
osol fits 368 Kb #################### < 1 s 20471.1 kb/s
osol fits 368 Kb #################### < 1 s 18189.3 kb/s
stat fits 361 Kb #################### < 1 s 17987.8 kb/s
osol fits 356 Kb #################### < 1 s 20627.4 kb/s
eph1 fits 308 Kb #################### < 1 s 15816.7 kb/s
eph1 fits 303 Kb #################### < 1 s 8951.9 kb/s
eph1 fits 283 Kb #################### < 1 s 11565.8 kb/s
aqual fits 250 Kb #################### < 1 s 15315.5 kb/s
cntr_img jpg 139 Kb #################### < 1 s 6587.8 kb/s
osol fits 129 Kb #################### < 1 s 5225.4 kb/s
vv pdf 43 Kb #################### < 1 s 2736.3 kb/s
full_img jpg 36 Kb #################### < 1 s 1980.1 kb/s
full_img fits 33 Kb #################### < 1 s 2099.4 kb/s
cntr_img fits 29 Kb #################### < 1 s 1810.2 kb/s
oif fits 20 Kb #################### < 1 s 1422.8 kb/s
bpix fits 12 Kb #################### < 1 s 338.5 kb/s
readme ascii 10 Kb #################### < 1 s 240.0 kb/s
eph1 fits 8 Kb #################### < 1 s 617.8 kb/s
fov fits 8 Kb #################### < 1 s 327.1 kb/s
flt fits 6 Kb #################### < 1 s 199.1 kb/s
msk fits 5 Kb #################### < 1 s 177.4 kb/s
pbk fits 4 Kb #################### < 1 s 217.3 kb/s
Total download size for ObsId 13858 = 66 Mb
Total download time for ObsId 13858 = 2 s
Students may want to uncompress, gunzip, all of the files in the top-level directory, along with the primary and secondary directories at this point.
In [3]:
gunzip 13858/primary/*gz
gunzip 13858/secondary/*gz
In [4]:
download_obsid_caldb 13858 ./CALDB
download_obsid_caldb
infile = 13858
outdir = ./CALDB
background = no
missing = no
clobber = no
verbose = 1
mode = ql
Retrieving files for CALDB_VER = 4.12.0
Retrieving CALDB index files
Processing infile=13858/primary/acisf13858N002_evt2.fits
Retrieving CALDB data files
Filename: 0------------------1
telD1999-07-23geomN0007.fits ####################
telD1999-07-23aimptsN0002.fits ####################
telD1999-07-23tdetN0001.fits ####################
telD1999-07-23skyN0002.fits ####################
telD1999-07-23sgeomN0001.fits ####################
hrmaD1996-12-20axeffaN0008.fits ####################
hrmaD1996-12-20vignetN0003.fits ####################
acisD1997-04-17qeN0006.fits ####################
acisD2010-02-01qeuN0007.fits ####################
acisD2000-11-28badpixN0004.fits ####################
acisD1999-08-13contamN0015.fits ####################
acisD1996-11-01gradeN0004.fits ####################
acisD2000-01-29grdimgN0001.fits ####################
acisD2000-01-29gain_ctiN0008.fits ####################
acisD2005-07-01evtspltN0002.fits ####################
acisD2010-01-01ctiN0012.fits ####################
acisD2012-05-01t_gainN0008.fits ####################
acisD1999-07-22subpixN0001.fits ####################
acisD2000-01-29fef_pha_ctiN0004.fits ####################
acisD1999-09-16dead_areaN0001.fits ####################
hrmaD1996-12-20reefN0001.fits ####################
acisD2000-01-29p2_respN0009_105-107.fits####################
acisD2000-01-29p2_respN0009_107-109.fits####################
acisD2000-01-29p2_respN0009_109-111.fits####################
acisD2000-01-29p2_respN0009_111-113.fits####################
acisD2000-01-29p2_respN0009_113-115.fits####################
acisD2000-01-29p2_respN0009_115-117.fits####################
acisD2000-01-29p2_respN0009_117-119.fits####################
acisD2000-01-29p2_respN0009_119-120.fits####################
acisD2002-11-15gtilimN0004.fits ####################
Be sure to source the new setup scripts to
use these CALDB files.
(t)csh:
source ./CALDB/software/tools/caldbinit.unix
bash:
source ./CALDB/software/tools/caldbinit.sh
In [5]:
export CALDB=`pwd`/CALDB
export CALDBCONFIG=`pwd`/CALDB/software/tools/caldb.config
export CALDBALIAS=none
Review V&V Report¶
Exercise 2¶
Read the Chandra Verification and Validation report in the top level directory of the data distribution: axaff13858N001_VV001_vv2.pdf
In [6]:
# okular 13858/axaff13858N002_VV001_vv2.pdf
Q: Summarize the V&V comments: Joint Proposal: NRAO
Q: What is the target of this observation? SDSS J091449.05+085321.1
Q: What is the sequence number SEQ_NUM of this observation? 702584
Display data in ds9¶
Locate the event files. There are two event files included in the standard data distribution. The Level 1 event file, evt1, is in the secondary/ directory, the Level 2 event file, evt2, is in the primary/ directory.
In [7]:
/bin/ls -1 13858/*/*evt*.fits
13858/primary/acisf13858N002_evt2.fits 13858/secondary/acisf13858_000N002_evt1.fits
Exercise 3¶
Display the Level 2 event file using ds9. Make sure to use Log scale.
In [8]:
ds9 13858/primary/acisf13858N002_evt2.fits -scale log &
sleep 10
[1] 268025
Q: Using the cursor, record the coordinate of the bright source in the center of the image:
In [9]:
xpaset -p ds9 crosshair 4091 4076
xpaset -p ds9 saveimage png ds9_01.png
display < ds9_01.png
In [10]:
xpaget ds9 crosshair wcs fk5 sexagesimal
xpaset -p ds9 mode region
9:14:49.0692 +8:53:20.454
Exercise 4¶
By default, ds9 only shows part of the Chandra dataset. Use the Bin menu to bin the data by a factor of 4 and then by a factor of 8.
Binning is different than Zooming. Zooming is done to the image (so after the event file has been binned). Set bin=1, and then zoom to ⅛.
In [11]:
xpaset -p ds9 bin factor 4
xpaset -p ds9 saveimage png ds9_exercise04_a.png
display < ds9_exercise04_a.png
In [12]:
xpaset -p ds9 bin factor 8
xpaset -p ds9 saveimage png ds9_exercise04_b.png
display < ds9_exercise04_b.png
In [13]:
xpaset -p ds9 bin factor 1
xpaset -p ds9 zoom 0.125
xpaset -p ds9 saveimage png ds9_exercise04_c.png
display < ds9_exercise04_c.png
Q: Describe why the Bin 8 image is different than the image Zoom'ed by ⅛:
**Zooming takes the default 1024x1024 image and samples every 8th row and column to create the display. Binning rebins the original event file with 8x8 pixels.***
Exercise 5¶
ds9 bins all the events in the event file into an image. That includes all events for all time and with all energies.
Use the Bin → Binning Parameters menu to add an energy=500:7000 as a Bin Filter.
In [14]:
xpaset -p ds9 zoom to 1
xpaset -p ds9 bin filter 'energy=500:7000'
xpaset -p ds9 saveimage png ds9_exercise05_a.png
display < ds9_exercise05_a.png
Q: Describe the difference in the images compared to the unfiltered image: Less background noise, max pixel value is lower (see colorbar scale compared to above)
Q: What unit is energy in? eV
Extra Credit¶
Try different energy ranges
In [15]:
xpaset -p ds9 bin filter 'energy=2000:10000'
xpaset -p ds9 saveimage png ds9_exercise05_extra1.png
display < ds9_exercise05_extra1.png
Try different time ranges
In [16]:
dmlist 13858/primary/acisf13858N002_evt2.fits header | grep TST
0039 TSTART 456518974.4871399999 [s] Real8 Observation start time (MET) 0040 TSTOP 456536836.1381000280 [s] Real8 Observation end time (MET) 0084 STARTBEP 301280942 Int4 BEP timer value at TSTART 0085 STOPBEP 1808980142 Int4 BEP timer value at TSTOP
In [17]:
xpaset -p ds9 bin filter 'time=456520000:456530000'
xpaset -p ds9 saveimage png ds9_exercise05_extra2.png
display < ds9_exercise05_extra2.png
Inspect headers¶
Most Chandra data products are FITS files, specifically FITS binary tables, with extensive headers that fully describe the dataset. Students should become familiar with some of the basic keywords and tool used to examine those keywords
Exercise 7¶
Open the Level 2 event file in prism
The EVENTS extension is automatically selected. The Header Keywords are shown in the top, right pane.
In [18]:
xpaset -p ds9 prism
In [19]:
#sleep 10
# import prism_exercise7.png
# display < prism_exercise7.png
MISSION: AXAFTELESCOP: CHANDRAINSTRUME: ACISDETNAM: ACIS-5678GRATING: NONEDATAMODE: FAINTREADMODE: TIMEDDATE-OBS: 2012-06-19T18:49:34OBSERVER: Dr. Kayhan GultekinONTIME: 15069.1001159
Q: Record the value for the following keywords:
Q: What units are the ACIS focal plane temperature, FP_TEMP, recorded in? K
In [20]:
xpaset -p ds9 quit
[1]+ Done ds9 13858/primary/acisf13858N002_evt2.fits -scale log
Exercise 8¶
Use dmlist with the header option to display the header.
In [21]:
dmlist 13858/primary/acisf13858N002_evt2.fits header,clean | \
egrep '^._TARG|^.*_NOM|^.*_PNT|^SIM_|^DTCOR|^ASCDSVER|^CALDBVER'
ASCDSVER 10.9.2 Processing system revision SIM_X -0.68282252473119 [mm] SIM focus pos SIM_Y 0 [mm] SIM orthogonal axis pos SIM_Z -190.1400660499 [mm] SIM translation stage pos RA_PNT 138.7036860339 [deg] Pointing RA DEC_PNT 8.8918212418 [deg] Pointing Dec ROLL_PNT 241.1580974095 [deg] Pointing Roll RA_NOM 138.7036860339 [deg] Nominal RA DEC_NOM 8.8918212418 [deg] Nominal Dec ROLL_NOM 241.1580974095 [deg] Nominal Roll DTCOR 0.98693426381071 Dead time correction CALDBVER 4.9.4
Q: Record the value for the following header keywords: see above
Q: This is an "ACIS-S" observation. How can you tell this from the event file header information? There is nothing in the header that says this is an 'ACIS-S' observation directly. The value of the
SIM_Zkeyword indicates that the aim point is located on the ACIS-S array
Extra Credit¶
In [22]:
dmkeypar 13858/primary/acisf13858N002_evt2.fits TIMEDEL echo+
3.14104
In [23]:
dmmakepar 13858/primary/acisf13858N002_evt2.fits[events] \
dmmakepar_exercise08.par clob+
In [24]:
pdump dmmakepar_exercise08.par | \
egrep -i '^._TARG|^.*_NOM|^.*_PNT|^SIM_|^DTCOR|^ASCDSVER|^CALDBVER'
ascdsver='10.9.2' sim_x='-0.68282252473119' sim_y='0' sim_z='-190.14006604987' ra_pnt='138.70368603385' dec_pnt='8.891821241756199' roll_pnt='241.15809740947' ra_nom='138.70368603385' dec_nom='8.891821241756199' roll_nom='241.15809740947' dtcor='0.98693426381071' caldbver='4.9.4'
In [25]:
plist dmmakepar_exercise08.par | \
egrep -i '^._TARG|^.*_NOM|^.*_PNT|SIM_|DTCOR|ASCDSVER|CALDBVER'
(ascdsver = 10.9.2) Processing system revision (sim_x = -0.68282252473119) [mm] SIM focus pos (sim_y = 0.0) [mm] SIM orthogonal axis pos (sim_z = -190.14006604987) [mm] SIM translation stage pos (ra_pnt = 138.70368603385) [deg] Pointing RA (dec_pnt = 8.8918212417562) [deg] Pointing Dec (roll_pnt = 241.15809740947) [deg] Pointing Roll (ra_nom = 138.70368603385) [deg] Nominal RA (dec_nom = 8.8918212417562) [deg] Nominal Dec (roll_nom = 241.15809740947) [deg] Nominal Roll (dtcor = 0.98693426381071) Dead time correction (caldbver = 4.9.4)
In [26]:
pget dmmakepar_exercise08.par \
sim_x sim_y sim_z ra_pnt dec_pnt roll_pnt ra_nom dec_nom roll_nom \
dtcor caldbver ascdsver
-0.68282252473119 0 -190.14006604987 138.70368603385 8.891821241756199 241.15809740947 138.70368603385 8.891821241756199 241.15809740947 0.98693426381071 4.9.4 10.9.2
In [27]:
dmhistory 13858/primary/acisf13858N002_evt2.fits acis_process_events
# dmhistory (CIAO 4.17.0): WARNING: Found and corrected "pixlib" library parameters
# dmhistory (CIAO 4.17.0): WARNING: Found and corrected "pixlib" library parameters
acis_process_events infile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/output/acisf13858_000N002_tmp_evt1.fits" outfile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/output/acisf13858_000N002_evt1.fits" acaofffile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/input/pcadf13858_000N001_asol1.fits[time=456518974.4871400:456536836.1381000]" apply_cti="yes" apply_tgain="yes" obsfile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/output/axaff13858_000N001_obs1.par" geompar="geom" logfile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/output/acis_process_events.log" gradefile="CALDB" grade_image_file="CALDB" gainfile="CALDB" badpixfile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/output/acisf13858_000N002_bpix1.fits" threshfile="CALDB" ctifile="CALDB" tgainfile="CALDB" mtlfile="/dsops/repro5/sdp.1/opus/prs_run/tmp//ACIS_F_L1_731033575n387/output/acisf13858_000N002_fptemp_egti1.fits" eventdef="{d:time,s:ccd_id,s:node_id,i:expno,s:chip,s:tdet,f:det,f:sky,s:phas,l:pha,l:pha_ro,f:energy,l:pi,s:fltgrade,s:grade,x:status}" doevtgrade="yes" check_vf_pha="no" trail="0.027" calculate_pi="yes" pi_bin_width="14.6" pi_num_bins="1024" max_cti_iter="15" cti_converge="0.1" clobber="no" verbose="0" stop="sky" rand_seed="1" rand_pha="yes" pix_adj="EDSER" subpixfile="CALDB" stdlev1="{d:time,l:expno,s:ccd_id,s:node_id,s:chip,s:tdet,f:det,f:sky,s:phas,l:pha,l:pha_ro,f:energy,l:pi,s:fltgrade,s:grade,x:status}" grdlev1="{d:time,l:expno,s:ccd_id,s:node_id,s:chip,s:tdet,f:det,f:sky,l:pha,l:pha_ro,s:corn_pha,f:energy,l:pi,s:fltgrade,s:grade,x:status}" cclev1="{d:time,d:time_ro,l:expno,s:ccd_id,s:node_id,s:chip,s:tdet,f:det,f:sky,f:sky_1d,s:phas,l:pha,l:pha_ro,f:energy,l:pi,s:fltgrade,s:grade,x:status}" ccgrdlev1="{d:time,d:time_ro,l:expno,s:ccd_id,s:node_id,s:chip,s:tdet,f:det,f:sky,f:sky_1d,l:pha,l:pha_ro,s:corn_pha,f:energy,l:pi,s:fltgrade,s:grade,x:status}" cclev1a="{d:time,d:time_ro,l:expno,s:ccd_id,s:node_id,s:chip,f:chipy_tg,f:chipy_zo,s:tdet,f:det,f:sky,f:sky_1d,s:phas,l:pha,l:pha_ro,f:energy,l:pi,s:fltgrade,s:grade,f:rd,s:tg_m,f:tg_lam,f:tg_mlam,s:tg_srcid,s:tg_part,s:tg_smap,x:status}" ccgrdlev1a="{d:time,d:time_ro,l:expno,s:ccd_id,s:node_id,s:chip,f:chipy_tg,f:chipy_zo,s:tdet,f:det,f:sky,f:sky_1d,l:pha,l:pha_ro,s:corn_pha,f:energy,l:pi,s:fltgrade,s:grade,f:rd,s:tg_m,f:tg_lam,f:tg_mlam,s:tg_srcid,s:tg_part,s:tg_smap,x:status}"
Reprocess dataset¶
The Chandra calibration database (CALDB) is continually updated. The then most recent CALDB is used for observations as they are taken. Some calibrations, such as the time dependent gain calibrations, can only be definitively computed based on historical observations; thus the file in the current CALDB is always predicted. These calibrations are later updated to be definitive post facto.
The CIAO team strongly advises users to always reprocess data obtained from the archive.
Exercise 9¶
Use chandra_repro to reprocess the dataset.
In [28]:
chandra_repro 13858 out= clob+
Running chandra_repro
version: 07 April 2025
Processing input directory '/lenin1.real/Junk/Workbook/13858'
No boresight correction update to asol file is needed.
Resetting afterglow status bits in evt1.fits file...
Running the destreak tool on the evt1.fits file...
Running acis_build_badpix and acis_find_afterglow to create a new bad pixel file...
Running acis_process_events to reprocess the evt1.fits file...
Output from acis_process_events:
# acis_process_events (CIAO 4.17.0): The following error occurred 2 times:
dsAPEPULSEHEIGHTERR -- WARNING: pulse height is less than split threshold when performing serial CTI adjustment.
Filtering the evt1.fits file by grade and status and time...
Applying the good time intervals from the flt1.fits file...
The new evt2.fits file is: /lenin1.real/Junk/Workbook/13858/repro/acisf13858_repro_evt2.fits
Updating the event file header with chandra_repro HISTORY record
Creating FOV file...
Setting observation-specific bad pixel file in local ardlib.par.
Cleaning up intermediate files
WARNING: Observation-specific bad pixel file set for session use:
/lenin1.real/Junk/Workbook/13858/repro/acisf13858_repro_bpix1.fits
Run 'punlearn ardlib' when analysis of this dataset completed.
Any issues pertaining to data quality for this observation will be listed in the Comments section of the Validation and Verification report located in:
/lenin1.real/Junk/Workbook/13858/repro/axaff13858N002_VV001_vv2.pdf
The data have been reprocessed.
Start your analysis with the new products in
/lenin1.real/Junk/Workbook/13858/repro
In [29]:
dmdiff 13858/primary/acisf13858N002_evt2.fits 13858/repro/acisf13858_repro_evt2.fits data- || echo
Infile 1: 13858/primary/acisf13858N002_evt2.fits Infile 2: 13858/repro/acisf13858_repro_evt2.fits ---------------------------------------------------------------------- Compare Headers ---------------------------------------------------------------------- Compare Key Lists: # dmdiff (CIAO 4.17.0): WARNING: keyword 'CHKVFPHA' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'STARTMJF' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'STARTMNF' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'STOPMJF' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'STOPMNF' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'CTI_TMAP' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'OBI_NUM' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'CLSTBITS' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'STATFILT' not present in infile1. Compare Keyword Details: # dmdiff (CIAO 4.17.0): WARNING: keyword 'ASCDSVER' comments differ. # dmdiff (CIAO 4.17.0): comment1="Processing system revision" # dmdiff (CIAO 4.17.0): comment2="ASCDS version number" # dmdiff (CIAO 4.17.0): WARNING: keyword 'CHECKSUM' comments differ. # dmdiff (CIAO 4.17.0): comment1="HDU checksum updated 2021-03-03T05:58:37" # dmdiff (CIAO 4.17.0): comment2="HDU checksum updated 2025-05-02T16:16:04" # dmdiff (CIAO 4.17.0): WARNING: keyword 'DATASUM' comments differ. # dmdiff (CIAO 4.17.0): comment1="data unit checksum updated 2021-03-03T05:58:27" # dmdiff (CIAO 4.17.0): comment2="data unit checksum updated 2025-05-02T16:16:03" Compare Keyword Values: Keyword: Message: Value(s): Diff: -------- -------------------------------------- ---------------------------------- ------------------------ CREATOR Values are not equal cxc - Version DS10.9 acis_process_events - CIAO 4.17.0 ASCDSVER Values are not equal 10.9.2 CIAO 4.17.0 CHECKSUM Values are not equal UGBQUF9QUFAQUF7Q 5eegAded8dedAded DATASUM Values are not equal 3437043554 3777912322 DATE Values are not equal 2021-03-03T05:58:25 2025-05-02T16:15:45 CTIFILE Values are not equal acisD2010-01-01ctiN0009.fits acisD2010-01-01ctiN0012.fits BPIXFILE Values are not equal acisf13858_000N002_bpix1.fits acisf13858_repro_bpix1.fits HISTNUM Values are not equal 533 672 +139 (+26.1%) ONTIME Values are not equal 15069.100115895 15070.556908429 +1.45679 (+0.00967%) ONTIME7 Values are not equal 15069.100115895 15070.556908429 +1.45679 (+0.00967%) ONTIME5 Values are not equal 15069.100115895 15070.556908488 +1.45679 (+0.00967%) ONTIME6 Values are not equal 15069.100115895 15070.556908429 +1.45679 (+0.00967%) ONTIME8 Values are not equal 15069.100115895 15070.556908488 +1.45679 (+0.00967%) LIVETIME Values are not equal 14872.211229171 14873.648987637 +1.43776 (+0.00967%) LIVTIME7 Values are not equal 14872.211229171 14873.648987637 +1.43776 (+0.00967%) LIVTIME5 Values are not equal 14872.211229171 14873.648987696 +1.43776 (+0.00967%) LIVTIME6 Values are not equal 14872.211229171 14873.648987637 +1.43776 (+0.00967%) LIVTIME8 Values are not equal 14872.211229171 14873.648987696 +1.43776 (+0.00967%) EXPOSURE Values are not equal 14872.211229171 14873.648987637 +1.43776 (+0.00967%) EXPOSUR7 Values are not equal 14872.211229171 14873.648987637 +1.43776 (+0.00967%) EXPOSUR5 Values are not equal 14872.211229171 14873.648987696 +1.43776 (+0.00967%) EXPOSUR6 Values are not equal 14872.211229171 14873.648987637 +1.43776 (+0.00967%) EXPOSUR8 Values are not equal 14872.211229171 14873.648987696 +1.43776 (+0.00967%) ---------------------------------------------------------------------- Compare Subspaces ---------------------------------------------------------------------- Compare Subspace Structure: Compare Column Details: # dmdiff (CIAO 4.17.0): WARNING: subspace column 'status' unsupported datatype (Bit), can not compare. # dmdiff (CIAO 4.17.0): WARNING: subspace column 'phas' units differ. # dmdiff (CIAO 4.17.0): unit1="" # dmdiff (CIAO 4.17.0): unit2="adu" # dmdiff (CIAO 4.17.0): WARNING: subspace column 'phas' comments differ. # dmdiff (CIAO 4.17.0): comment1="" # dmdiff (CIAO 4.17.0): comment2="array of pixel pulse heights" Compare Subspace Ranges: component 1 Column: Range Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ time range 1 min. Values are not equal 456520730.637154 456520729.942298 -0.694856 (-1.52e-07%) time range 1 max. Values are not equal 456535799.737269 456535800.499206 +0.761937 (+1.67e-07%) expno range 1 min. Values are not equal 0 3 +3 ( +inf%) expno range 1 max. Values are not equal 2147483647 4801 -2.14748e+09 ( -100%) Compare Subspace Ranges: component 2 Column: Range Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ time range 1 min. Values are not equal 456520730.637154 456520729.983338 -0.653816 (-1.43e-07%) time range 1 max. Values are not equal 456535799.737269 456535800.540246 +0.802977 (+1.76e-07%) expno range 1 min. Values are not equal 0 3 +3 ( +inf%) expno range 1 max. Values are not equal 2147483647 4801 -2.14748e+09 ( -100%) Compare Subspace Ranges: component 3 Column: Range Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ time range 1 min. Values are not equal 456520730.637154 456520730.024378 -0.612776 (-1.34e-07%) time range 1 max. Values are not equal 456535799.737269 456535800.581286 +0.844017 (+1.85e-07%) expno range 1 min. Values are not equal 0 3 +3 ( +inf%) expno range 1 max. Values are not equal 2147483647 4801 -2.14748e+09 ( -100%) Compare Subspace Ranges: component 4 Column: Range Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ time range 1 min. Values are not equal 456520730.637154 456520730.065418 -0.571736 (-1.25e-07%) time range 1 max. Values are not equal 456535799.737269 456535800.622326 +0.885057 (+1.94e-07%) expno range 1 min. Values are not equal 0 3 +3 ( +inf%) expno range 1 max. Values are not equal 2147483647 4801 -2.14748e+09 ( -100%) Compare Subspace Regions: component 1 Compare Subspace Regions: component 2 Compare Subspace Regions: component 3 Compare Subspace Regions: component 4
In [30]:
dmkeypar 13858/primary/acisf13858N002_evt2.fits CALDBVER echo+
dmkeypar 13858/repro/acisf13858_repro_evt2.fits CALDBVER echo+
check_ciao_caldb
4.9.4
4.9.4
CALDB environment variable = /lenin1.real/Junk/Workbook/CALDB
CALDB version = 4.12.0
release date = 2025-01-30T16:00:00 UTC
CALDB query completed successfully.
In [31]:
dmlist 13858/primary/acisf13858N002_evt2.fits counts
dmlist 13858/repro/acisf13858_repro_evt2.fits counts
134011 134015
Q: Compare the header keyword values in the _repro_evt2.fits file with the evt2 file obtained from the archive. Discuss the differences: We see newer calibration files are used (
GAINFILE,CTIFILE, andTGAINFIL). We see small diff inONTIME.
Q: What version of the CALDB is installed? What is the value of the CALDBVER keyword? It is unchanged.
CALDBVERis not updated by the tools; it's set in SDP.
Q: Compare the number of events in the reprocessed Level 2 event file with the number of events in the archived Level 2 event file. Why are they same (or not the same)? **Slightly different. New calibrations mean some good events go bad, some bad events go good (grade migration) **
Imaging Basics (Spatial Analysis)¶
In this section students will exercise some of the basic CIAO tools and scripts needed to perform basic imaging tasks.
Detect Sources¶
One of the most common analysis tasks is to detect sources.
Please keep in mind: CIAO detect tools only detect candidate sources and they are not photometric tools. The output from the detect tools should be used to guide further analysis.
See also:
- Wavdetect thread: http://cxc.cfa.harvard.edu/ciao/threads/wavdetect/
- Using detect output: http://cxc.cfa.harvard.edu/ciao/threads/detect_output/
Exercise 10¶
In this exercise students should complete the following steps
Create image, exposure map, and psf map. Run the fluximage script on the dataset. Students should be sure to set the binsize=1 and set bands=default or bands=broad.
Detect Sources. Run wavdetect on the fluximage output counts image, using the exposure map and psfmap. Students should select a set of wavelet scales to use that is appropriate for the dataset being analyzed.
In [32]:
fluximage 13858/repro/acisf13858_repro_evt2.fits out=acisf13858 bin=1 bands=broad clob+ mode=h psfecf=0.393
Running fluximage
Version: 04 November 2021
Using CSC ACIS broad science energy band.
Aspect solution 13858/repro/pcadf13858_000N001_asol1.fits found.
Bad-pixel file 13858/repro/acisf13858_repro_bpix1.fits found.
Mask file 13858/repro/acisf13858_000N002_msk1.fits found.
The output images will have 2940 by 4168 pixels, pixel size of 0.492 arcsec,
and cover x=2720.5:5660.5:1,y=1800.5:5968.5:1.
Running tasks in parallel with 4 processors.
Creating 4 aspect histograms for obsid 13858
Creating 4 instrument maps for obsid 13858
Creating 4 exposure maps for obsid 13858
Combining 4 exposure maps for obsid 13858
Thresholding data for obsid 13858
Exposure-correcting image for obsid 13858
Creating PSF map for obsid 13858
The following files were created:
The clipped counts image is:
acisf13858_broad_thresh.img
The observation FOV is:
acisf13858.fov
The clipped exposure map is:
acisf13858_broad_thresh.expmap
The PSF map is:
acisf13858_broad_thresh.psfmap
The exposure-corrected image is:
acisf13858_broad_flux.img
In [33]:
pset wavdetect \
infile=acisf13858_broad_thresh.img \
expfile=acisf13858_broad_thresh.expmap \
psffile=acisf13858_broad_thresh.psfmap \
scales='1.4 2 4 8 12 16 24 32' \
outfile=acisf13858_wav.src \
scellfile=acisf13858_wav.cell \
imagefile=acisf13858_wav.recon \
defnbkg=acisf13858_wav.nbkg \
clob+
wavdetect mode=h
In [34]:
ds9 acisf13858_broad_thresh.psfmap -scale linear -cmap standard -zoom to fit \
-saveimage png ds9_exercise10_a.png -quit
display < ds9_exercise10_a.png
In [35]:
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region acisf13858_wav.src \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_b.png -quit
display < ds9_exercise10_b.png
Q: Display the PSF map in ds9. Why does it look the way it does? Chandra PSF increases in size the further away from optical axis.
Q: Record choice of wavelet scales and why that set was selected: above. Pseduo sqrt(2) used to detect srcs across field
Q: Open the broad-band image in ds9 and load the source list as a region file. Comment on the detected sources and any missed detections: src detects look good ; regions are a little small; maybe some faint stuff missed
Q: Are the ellipses the position error or the apparent size of the source? They are a measure of the observed source size. They are not the intrinsic (deconvolved) source size, nor are they the position error.
Extra Credit¶
Try differenct wavelet scales
In [36]:
wavdetect scales='1 4 16 64' mode=h clobber=yes outfile=fewscales.src
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region fewscales.src \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_ec1.png -quit
display < ds9_exercise10_ec1.png
Try different significance thresholds, sigthresh.
In [37]:
wavdetect sigthresh=1e-7 mode=h clobber=yes outfile=lowersigthresh.src
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region lowersigthresh.src \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_ec2.png -quit
display < ds9_exercise10_ec2.png
Try different ecf values when making the PSF map.
In [38]:
mkpsfmap acisf13858_broad_thresh.img ecf90.psfmap \
energy=1.4967 ecf=0.90 mode=h clob+
wavdetect mode=h clobber=yes outfile=ecf90.src psffile=ecf90.psfmap
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region ecf90.src \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_ec3.png -quit
display < ds9_exercise10_ec3.png
Try different energy bands.
In [39]:
mkpsfmap acisf13858_broad_thresh.img 5keV.psfmap \
energy=5.0 ecf=0.90 mode=h clob+
wavdetect mode=h clobber=yes outfile=5keV.src psffile=5keV.psfmap
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region 5keV.src \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_ec4.png -quit
display < ds9_exercise10_ec4.png
Try using the celldetect tool.
In [40]:
celldetect acisf13858_broad_thresh.img \
expstk=acisf13858_broad_thresh.expmap \
psffile=ecf90.psfmap \
out=acisf13858_cell.src \
clob+ mode=h
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region acisf13858_cell.src \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_ec5.png -quit
display < ds9_exercise10_ec5.png
Try using the vtpdetect tool.
In [41]:
vtpdetect acisf13858_broad_thresh.img \
exp=acisf13858_broad_thresh.expmap \
out=acisf13858_vtp.src \
clob+ mode=h
ds9 acisf13858_broad_thresh.img -block 2 -scale log -region acisf13858_vtp.src[src_region] \
-scale limits 0 20 -smooth -saveimage png ds9_exercise10_ec6.png -quit
display < ds9_exercise10_ec6.png
Create regions in ds9¶
The detect tools are one way to create regions automatically. However, often users analyzing a single source in an observation will simply create their own region in ds9 and then use those regions with the CIAO tools.
Exercise 11¶
- Open the broad-band counts image in ds9.
- Create a source region. Draw a circular source region around the bright source in the center of the image.
- Save the region, ds9_src.reg. Use ds9 format and celestial (world) coordinates.
- Delete the source region.
- Create a background region. Draw an annulus around the same source. The inner annulus region should be larger than the source region. The outer region should be large enough to get a large number of counts, but should be small enough to avoid any nearby source or any instrumental features (like the gaps between chips).
- Save the region, ds9_bkg.reg. Use CIAO format and physical coordinates.
In [47]:
ds9 acisf13858_broad_thresh.img -scale log -pan to 4096.5 4096.5 physical &
sleep 10
[1] 274880
In [48]:
echo "physical;circle(4090.9,4077.6,10)" | xpaset ds9 regions -format ds9
xpaset -p ds9 saveimage png ds9_exercise11_a.png
xpaget ds9 regions -format ds9 -system wcs -skyformat sexagesimal > ds9_src.reg
xpaset -p ds9 regions delete all
display < ds9_exercise11_a.png
cat ds9_src.reg
# Region file format: DS9 version 4.1 global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1 fk5 circle(9:14:49.0725,+8:53:21.241,4.920")
In [49]:
echo "physical;annulus(4090.9,4077.6,45,125) # background" | xpaset ds9 region -format ds9
xpaset -p ds9 saveimage png ds9_exercise11_b.png
xpaget ds9 regions -format ciao -system physical > ds9_bkg.reg
xpaset -p ds9 regions delete all
display < ds9_exercise11_b.png
cat ds9_bkg.reg
annulus(4090.9,4077.6,45,125)
In [50]:
xpaset -p ds9 quit
[1]+ Done ds9 acisf13858_broad_thresh.img -scale log -pan to 4096.5 4096.5 physical
Q: Record the source location used: see above
Q: Describe the differences between the CIAO format and the ds9 format regions: ds9 format includes meta-data: color, linestyle, etc.
Exercise 12¶
Compute source centroid. In this exercise students will use the dmstat tool to compute the source centroid in their region.
- Use dmcopy to filter the flux image thresh.img using the source region file.
- Use dmlist with the blocks option to display info about the filtered output image.
- Use dmstat with the centroid=yes option to compute the centroid of the counts in the filtered image.
- Use dmcopy to filter the Level 2 event file using the source region file and to bin it into an image with binsize=1.
- Use dmlist with the blocks option to display info about the filtered and binned image.
- Use dmstat with centroid=yes option to compute the centroid of the counts in the filtered and binned image.
In [51]:
dmcopy "acisf13858_broad_thresh.img[sky=region(ds9_src.reg)]" dmcopy_e12_src.fits clob+
In [52]:
dmlist dmcopy_e12_src.fits blocks,cols
--------------------------------------------------------------------------------
Dataset: dmcopy_e12_src.fits
--------------------------------------------------------------------------------
Block Name Type Dimensions
--------------------------------------------------------------------------------
Block 1: EVENTS_IMAGE Image Int4(21x21)
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
ColNo Name Unit Type Range
1 EVENTS_IMAGE[21,21] Int4(21x21) -
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
Group# Axis#
1 1,2 sky(x) = (+4080.50) +(+1.0)* ((#1)-(+0.50))
(y) (+4067.50) (+1.0) ((#2) (+0.50))
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
Group# Axis#
1 1,2 EQPOS(RA ) = (+138.7037)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
(DEC) (+8.8918 ) (+0.000136667) ( (y) (+4096.50))
In [53]:
dmstat dmcopy_e12_src.fits cen+ sig- med+
EVENTS_IMAGE(x, y)
min: 0 @: ( 4088 4068 )
max: 382 @: ( 4091 4078 )
cntrd[log] : ( 10.870538838 10.561931421 )
cntrd[phys]: ( 4090.8705388 4077.5619314 )
good: 316
null: 125
In [54]:
dmcopy 13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_src.reg)][bin sky=1]" \
dmcopy_e12_evtsrc.fits clob+
In [55]:
dmlist dmcopy_e12_evtsrc.fits blocks,cols
--------------------------------------------------------------------------------
Dataset: dmcopy_e12_evtsrc.fits
--------------------------------------------------------------------------------
Block Name Type Dimensions
--------------------------------------------------------------------------------
Block 1: EVENTS_IMAGE Image Int2(20x20)
Block 2: GTI7 Table 2 cols x 1 rows
Block 3: GTI5 Table 2 cols x 1 rows
Block 4: GTI6 Table 2 cols x 1 rows
Block 5: GTI8 Table 2 cols x 1 rows
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
ColNo Name Unit Type Range
1 EVENTS_IMAGE[20,20] Int2(20x20) -
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
Group# Axis#
1 1,2 sky(x) = (+4080.8416) +(+1.0)* ((#1)-(+0.50))
(y) (+4067.5661) (+1.0) ((#2) (+0.50))
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
Group# Axis#
1 1,2 EQPOS(RA ) = (+138.7037)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
(DEC) (+8.8918 ) (+0.000136667) ( (y) (+4096.50))
In [56]:
dmstat dmcopy_e12_evtsrc.fits cen+ sig- med+
EVENTS_IMAGE(x, y)
min: 0 @: ( 4088.3415984 4068.0661232 )
max: 361 @: ( 4090.3415984 4077.0661232 )
cntrd[log] : ( 10.459535937 10.506508206 )
cntrd[phys]: ( 4090.8011343 4077.5726314 )
good: 316
null: 84
Q: Is the filtered image the same size as the fluximage output image size? Why? The CXCDM shrinks the image to the size of the region (by default)
Q: Record the centroid of fluximage output filtered image: 4104.4755245, 4115.8174825
Q: What are the good and null values that dmstat reports? The number of pixels outside the circle
Q: Is the image created by filtering and binning the event file the same size as the other images? Why? No. (20x20 vs 21x21). Filter image vs. Filter table then bin into an image
Q: Record the centroid of the event file filter and binned image: 4104.4519153, 4115.808491
Q: Is the centroid the same for the two images? Why? (Hint: what energies are being used?) Opps, no energy filter on the event file -- but image was made with broad band (500:7000eV).
Aperture Photometry¶
Obtaining the counts (or count rate) is the first step in computing the source flux.
Exercise 13¶
In this exercise students will use the CIAO analysis menu, aka dax, to get the net counts in their regions (from Exercise 11).
- Open the fluximage thresh.img in ds9.
- Load the source region file created in Exercise 11 step 3.
- Load the background region file created in Exercise 11 step 6.
- Double click on the background annulus to display the region properties.
- Under the Property menu, select Background. The region will now be drawn with a dashed line. Close the properties window.
- In the main ds9 window, goto Analysis →CIAO →Statistics → Net Counts.
- A text window will be display containing the information about the current regions.
- Try adjusting the source and background radii and repeating the analysis.
In [57]:
ds9 acisf13858_broad_thresh.img -scale log -pan to 4096.5 4096.5 physical &
sleep 10
[1] 274946
In [58]:
xpaset -p ds9 regions ds9_bkg.reg
xpaset -p ds9 regions select all
xpaset -p ds9 regions background
In [59]:
xpaset -p ds9 regions ds9_src.reg -format ds9
xpaset -p ds9 saveimage png exercise_13_ds9a.png
display < exercise_13_ds9a.png
In [60]:
xpaset -p ds9 regions select all
xpaset -p ds9 analysis task "{Net Counts}"
sleep 10
In [61]:
# sleep 10
# import exercise_13_a.png
# display < exercise_13_a.png
In [62]:
xpaset -p ds9 regions delete all
echo "physical; circle(4090.9,4077.6,17.145457)" | xpaset ds9 regions -format ds9
echo "physical; annulus(4090.9,4077.6,26.305893,125) # background" | xpaset ds9 regions -format ds9
xpaset -p ds9 saveimage png exercise_13_ds9b.png
display < exercise_13_ds9b.png
In [63]:
xpaset -p ds9 regions select all
xpaset -p ds9 analysis task "{Net Counts}"
sleep 10
In [64]:
# sleep 10
# import exercise_13_b.png
# display < exercise_13_b.png
In [65]:
xpaset -p ds9 quit
[1]+ Done ds9 acisf13858_broad_thresh.img -scale log -pan to 4096.5 4096.5 physical
- Net counts with error: 1428 +/- 37.85
- Source counts: 1430
- Background counts: 270
Q: From step 7, record the following information from the Net Counts task:
Q: From step 8, describe how the net counts and net rates vary with changes to the source and background regions: net rate doesn't change much as bkg gets bigger/smaller (ie background is flat)
Exercise 14¶
In this exercise students will use the srcflux script to get various estimates of the flux of the source.
- Obtain the source location in celestial coordinates from the source region file created in Exercise 11 step 3.
- Run the srcflux script on the reprocessed level 2 event file, using the position, pos, from step 1. All other parameters should remain at their defaults.
- Repeat step 2, using a different outroot and psfmethod=arfcorr
- Repeat step 3, using a different outroot, model="xsbbody.black_body" and paramvals="black_body.kT=1".
- Repeat step 4, using bands="hard".
In [66]:
punlearn dmcoords
dmcoords 13858/repro/acisf13858_repro_evt2.fits op=sky x=4090.8011343 y=4077.5726314 celfmt=hms verb=0
pget dmcoords ra dec
09:14:49.073 +08:53:21.24
In [67]:
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.073, +08:53:21.24" exercise14_step2 clob+ mode=h
srcflux (17 April 2024)
infile = 13858/repro/acisf13858_repro_evt2.fits
pos = 09:14:49.073, +08:53:21.24
outroot = exercise14_step2
bands = default
regions = simple
srcreg =
bkgreg =
bkgresp = yes
psfmethod = ideal
psffile =
conf = 0.9
binsize = 1
rmffile =
arffile =
model = xspowerlaw.pow1
paramvals = pow1.PhoIndex=2.0
absmodel = xsphabs.abs1
absparams = abs1.nH=%GAL%
abund = angr
pluginfile =
fovfile =
asolfile =
mskfile =
bpixfile =
dtffile =
ecffile = CALDB
marx_root =
parallel = yes
nproc = INDEF
tmpdir = /tmp
random_seed = -1
clobber = yes
verbose = 1
mode = h
Processing OBI 001
Extracting counts
Setting Ideal PSF : alpha=1 , beta=0
Getting net rate and confidence limits
Getting model independent fluxes
Getting model fluxes
Getting photon fluxes
Getting variability
Running tasks in parallel with 4 processors.
Running aprates for exercise14_step2_0001_broad_rates.par
Running eff2evt for exercise14_step2_broad_0001_src.dat
Running eff2evt for exercise14_step2_broad_0001_bkg.dat
Making response files for exercise14_step2_0001
Running fluximage for exercise14_step2_0001
Making Lightcurve for source 1
Running modeflux for region 1
Using GAL=0.0426 for source 1
Adding net rates to output
Appending flux results onto output
Appending photflux results onto output
Computing Net fluxes
Adding model fluxes to output
Scaling model flux confidence limits
Appending variability results onto output
Summary of source fluxes
Position 0.5 - 7.0 keV
Value 90% Conf Interval
#0001|9 14 49.07 +8 53 21.2 Rate 0.0831 c/s (0.0792,0.087)
Flux 5.6E-13 erg/cm2/s (5.34E-13,5.87E-13)
Mod.Flux 5.82E-13 erg/cm2/s (5.55E-13,6.1E-13)
Unabs Mod.Flux 6.25E-13 erg/cm2/s (5.95E-13,6.54E-13)
In [68]:
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.073, +08:53:21.24" exercise14_step3 \
psfmethod=arfcorr \
clob+ mode=h verbose=0
cat exercise14_step3_summary.txt
Summary of source fluxes
Position 0.5 - 7.0 keV
Value 90% Conf Interval
#0001|9 14 49.07 +8 53 21.2 Rate 0.0967 c/s (0.0922,0.101)
Flux 6.52E-13 erg/cm2/s (6.22E-13,6.83E-13)
Mod.Flux 6.78E-13 erg/cm2/s (6.46E-13,7.1E-13)
Unabs Mod.Flux 7.27E-13 erg/cm2/s (6.93E-13,7.61E-13)
In [69]:
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.073, +08:53:21.24" exercise14_step4 \
psfmethod=arfcorr \
model="xsbbody.black_body" paramvals="black_body.kT=1" \
clob+ mode=h verbose=0
cat exercise14_step4_summary.txt
Summary of source fluxes
Position 0.5 - 7.0 keV
Value 90% Conf Interval
#0001|9 14 49.07 +8 53 21.2 Rate 0.0967 c/s (0.0922,0.101)
Flux 6.52E-13 erg/cm2/s (6.22E-13,6.83E-13)
Mod.Flux 1.02E-12 erg/cm2/s (9.68E-13,1.06E-12)
Unabs Mod.Flux 1.03E-12 erg/cm2/s (9.84E-13,1.08E-12)
In [70]:
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.073, +08:53:21.24" exercise14_step5 \
psfmethod=arfcorr \
model="xsbbody.black_body" paramvals="black_body.kT=1" \
band="hard" \
clob+ mode=h verbose=0
cat exercise14_step5_summary.txt
Summary of source fluxes
Position 2.0 - 7.0 keV
Value 90% Conf Interval
#0001|9 14 49.07 +8 53 21.2 Rate 0.0186 c/s (0.0166,0.0206)
Flux 2.7E-13 erg/cm2/s (2.4E-13,3E-13)
Mod.Flux 2.92E-13 erg/cm2/s (2.6E-13,3.24E-13)
Unabs Mod.Flux 2.94E-13 erg/cm2/s (2.62E-13,3.26E-13)
Q: Record the net count rate, model independent flux, model flux from steps 2 through 5
.
| Estimate | Step 2 | Step 3 | Step 4 | Step 5 |
|---|---|---|---|---|
| Count Rate | 0.827 | 0.962 | 0.962 | 0.0185 |
| Flux | 5.6E-13 | 7.3E-13 | 7.3E-13 | 3.12E-13 |
| Model Flux | 5.8E-13 | 6.74E-13 | 1.0E-12 | 2.91E-13 |
.
Q: Discuss the differences in the estimated fluxes obtained in steps 2 through 5.
Extra Credit¶
Try using psfmethod=quick and psfmethod=marx
In [71]:
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.088, +08:53:21.16" exercise14_step_ec1 \
psfmethod=quick \
clob+ mode=h verbose=0
cat exercise14_step_ec1_summary.txt
Summary of source fluxes
Position 0.5 - 7.0 keV
Value 90% Conf Interval
#0001|9 14 49.08 +8 53 21.1 Rate 0.0928 c/s (0.0884,0.0973)
Flux 6.27E-13 erg/cm2/s (5.97E-13,6.57E-13)
Mod.Flux 6.51E-13 erg/cm2/s (6.2E-13,6.82E-13)
Unabs Mod.Flux 6.98E-13 erg/cm2/s (6.65E-13,7.31E-13)
In [72]:
#source $ASCDS_INSTALL/marx-5.5.0/setup_marx.sh
export MARX_ROOT=$ASCDS_INSTALL
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.088, +08:53:21.16" exercise14_step_ec2 \
psfmethod=marx \
clob+ mode=h verbose=0
cat exercise14_step_ec2_summary.txt
Summary of source fluxes
Position 0.5 - 7.0 keV
Value 90% Conf Interval
#0001|9 14 49.08 +8 53 21.1 Rate 0.0932 c/s (0.0888,0.0976)
Flux 6.29E-13 erg/cm2/s (5.99E-13,6.59E-13)
Mod.Flux 6.53E-13 erg/cm2/s (6.22E-13,6.85E-13)
Unabs Mod.Flux 7.01E-13 erg/cm2/s (6.67E-13,7.34E-13)
| Estimate | ideal | quick | arfcorr | marx |
|---|---|---|---|---|
| Count Rate | 0.827 | 0.962 | 0.955 | 0.0985 |
| Flux | 5.6E-13 | 7.3E-13 | 6.5E-13 | 7.84E-13 |
| Model Flux | 5.8E-13 | 6.74E-13 | 6.7E-12 | 6.90E-13 |
Try other spectral model and other model parameters.
In [ ]:
Try other energy bands (2-10 keV is a common band in the literature).
In [73]:
punlearn srcflux
srcflux 13858/repro/acisf13858_repro_evt2.fits "09:14:49.073, +08:53:21.24" exercise14_step_ec4 \
psfmethod=quick band="2.0:10.0:3.0" \
clob+ mode=h verbose=0
cat exercise14_step_ec4_summary.txt
Summary of source fluxes
Position 2.0 - 10.0 keV
Value 90% Conf Interval
#0001|9 14 49.07 +8 53 21.2 Rate 0.0183 c/s (0.0163,0.0203)
Flux 3.19E-13 erg/cm2/s (2.84E-13,3.53E-13)
Mod.Flux 3.66E-13 erg/cm2/s (3.27E-13,4.06E-13)
Unabs Mod.Flux 3.68E-13 erg/cm2/s (3.28E-13,4.08E-13)
Simulate PSF¶
The Chandra Point Spread Function (PSF) varies spatially across the field of view as well with energy. There is no analytic model of the PSF. The only way to obtain an estimate of the Chandra PSF for detailed spatial analysis is by simulation.
The psfmethod used by srcflux are only appropriate when integrating over a region. The pixel-to-pixel variations used in those methods is insufficient for any type of analysis that requires detailed spatial information.
There are two different simulators users can use.
- SAOTrace is the definitive mirror model and is made available to users via the ChaRT web interface.
- MARX is the end-to-end Chandra simulator used to calibrate the HETG gratings and includes a fairly accurate mirror model and detector models.
Exercise 15¶
In this exercise Students will simulate the PSF using ChaRT and with MARX. The simulate_psf script is used to create PSF images that are appropriate for further data analysis.
- Ensure that MARX is installed and that MARX_ROOT is set.
- Run ChaRT. http://cxc.cfa.harvard.edu/ciao/PSFs/chart2/index.html To run ChaRT users need to know the source location (celestial coordinates), an estimate of the photon flux (from the srcflux .flux output files), and an estimate of the energy. Students may be provided with a set of ChaRT output ray files.
- Run the simulate_psf script. Use the reprocessed level 2 event file as the infile and the coordinates of the source, in decimal degrees, for the ra and dec parameters. Use simulator=file and provide the stack of ray files obtained from ChaRT for the rayfile parameter. Set minsize=128 to obtain a reasonable size output image.
- Display the reprocessed level 2 event file in one ds9 frame and the simulate_psf output PSF in a 2nd frame. Tile the frames horizontally and save the image as a PNG file.
- Run simulate_psf with a different outroot, this time using simulator=marx and set the numiter parameter to the same number of iterations as were used with ChaRT. monoenergy=1.0 and flux=1e-4 can be used for the spectral model.
- Repeat step 4 using the MARX generated PSF.
In [74]:
#source $ASCDS_INSTALL/marx-5.5.0/setup_marx.sh
In [75]:
dmkeypar exercise14_step2_broad.flux net_photflux_aper echo+
0.000226808610552178
In [76]:
curl \
-F email=kglotfelty@cfa.harvard.edu \
-F coords=cel \
-F ra=09:14:49.073 \
-F dec=+08:53:21.24 \
-F asol=obi \
-F obsid=13858 \
-F obinum=0\
-F how_many=numiter \
-F niter=10\
-F randseed=32767 \
-F energy=mono \
-F mono=2.3 \
-F flux=0.000221 \
https://saotrace.cfa.harvard.edu/cgi-bin/runwrapper
<html>
<head><title>ChaRT Job Submitted</title></head>
<body>
<h1>ChaRT v2</h1>
Job submitted. You will recieve an email from cxc_rays@head.cfa.harvard.edu when the job is complete or if there are any errors.
</body>
</html>
In [77]:
curl -o rays.tar.gz https://saotrace.cfa.harvard.edu/pickup/HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered.tar.gz
tar xvfz rays.tar.gz
% Total % Received % Xferd Average Speed Time Time Time Current
Dload Upload Total Spent Left Speed
100 1871k 100 1871k 0 0 46.8M 0 --:--:-- --:--:-- --:--:-- 46.8M
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0001_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0001_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0001_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0005_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0007_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0003_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0009_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0004_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0006_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0000_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0000_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0000_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0008_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0004_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0005_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0003_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0008_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0005_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0007_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0002_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0002_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0002_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0008_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0006_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0009_rays.tot_wt.in.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0003_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0007_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0006_rays.tot_wt.out.log
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0004_rays.fits
HRMA_ra138.70447_dec8.88923_en2.3_flux0.000221_dithered_i0009_rays.tot_wt.out.log
In [78]:
simulate_psf 13858/repro/acisf13858_repro_evt2.fits \
chart_sim ra=138.70453 dec=8.88921 \
simulator=file rayfile="HRMA*rays.fits" \
minsize=128 mode=h
simulate_psf
infile = 13858/repro/acisf13858_repro_evt2.fits
outroot = chart_sim
ra = 138.70453
dec = 8.88921
spectrumfile =
monoenergy = INDEF
flux = INDEF
simulator = file
rayfile = HRMA*rays.fits
projector = marx
random_seed = -1
blur = 0.07000000000000001
readout_streak = no
pileup = no
ideal = yes
extended = yes
binsize = 1
numsig = 7
minsize = 128
maxsize = INDEF
numiter = 1
numrays = INDEF
keepiter = no
asolfile =
marx_root = /export/miniforge/envs/ciao-4.17
verbose = 1
mode = h
Started check_setup
Finished check_setup
10 rayfiles provided; ignoring numiter=1 parameter
Performing iteration 1 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 2 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 3 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 4 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 5 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 6 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 7 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 8 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 9 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 10 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Started create_average_image
Finished create_average_image
Final output PSF image is : chart_sim.psf[PSF]
In [79]:
ds9 acisf13858_broad_thresh.img -scale log -pan to 4096.5 4096.5 physical \
-zoom to 4 \
chart_sim.psf -pan to 4096.5 4096.5 physical \
-view colorbar no \
-saveimage png ds9_exercise15_a.png -quit
display < ds9_exercise15_a.png
In [80]:
simulate_psf 13858/repro/acisf13858_repro_evt2.fits \
marx_sim ra=138.70453 dec=8.88921 \
simulator=marx \
minsize=128 mode=h \
mono=2.3 flux=0.000221 numiter=10
simulate_psf
infile = 13858/repro/acisf13858_repro_evt2.fits
outroot = marx_sim
ra = 138.70453
dec = 8.88921
spectrumfile =
monoenergy = 2.3
flux = 0.000221
simulator = marx
rayfile = @foo.lis
projector = marx
random_seed = -1
blur = 0.07000000000000001
readout_streak = no
pileup = no
ideal = yes
extended = yes
binsize = 1
numsig = 7
minsize = 128
maxsize = INDEF
numiter = 10
numrays = INDEF
keepiter = no
asolfile =
marx_root = /export/miniforge/envs/ciao-4.17
verbose = 1
mode = h
Started check_setup
Finished check_setup
Performing iteration 1 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 2 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 3 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 4 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 5 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 6 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 7 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 8 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 9 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Performing iteration 10 of 10
Started run_marx
Finished run_marx
Started create_psf_image
Finished create_psf_image
Started create_average_image
Finished create_average_image
Final output PSF image is : marx_sim.psf[PSF]
In [81]:
ds9 acisf13858_broad_thresh.img -scale log -pan to 4096.5 4096.5 physical \
-zoom to 4 \
marx_sim.psf -pan to 4096.5 4096.5 physical \
-view colorbar no \
-saveimage png ds9_exercise15_b.png -quit
display < ds9_exercise15_b.png
Q: Comment on a visual comparison of the ChaRT/SAOTrace PSF to the observation:
Q: Comment on a visual comparison of the MARX PSF to the observation:
Q: Describe any differences between the MARX and ChaRT PSFs:
Q: Is this a point source?
Extra Credit¶
Try setting simulate_psf readout_streak=yes. Try pileup=yes. Try blur=0.20
In [82]:
simulate_psf 13858/repro/acisf13858_repro_evt2.fits \
marxA_sim ra=138.70453 dec=8.88921 \
simulator=marx \
minsize=128 mode=h \
mono=2.3 flux=0.000221 numiter=10 verb=0 readout+
In [83]:
simulate_psf 13858/repro/acisf13858_repro_evt2.fits \
marxB_sim ra=138.70453 dec=8.88921 \
simulator=marx \
minsize=128 mode=h \
mono=2.3 flux=0.000221 numiter=10 verb=0 pileup+ ext-
In [84]:
simulate_psf 13858/repro/acisf13858_repro_evt2.fits \
marxC_sim ra=138.70453 dec=8.88921 \
simulator=marx \
minsize=128 mode=h \
mono=2.3 flux=0.000221 numiter=10 verb=0 blur=0.2
In [85]:
ds9 marx_sim.psf -scale log -pan to 4090.8 4077.6 physical \
-zoom to 4 \
marxA_sim.psf -pan to 4090.8 4077.6 physical \
marxB_sim.psf -pan to 4090.8 4077.6 physical \
marxC_sim.psf -pan to 4090.8 4077.6 physical \
-view colorbar no \
-saveimage png ds9_exercise15_ec1.png -quit
display < ds9_exercise15_ec1.png
Run the Lucy-Richardson deconvolution tool arestore using the broad band flux image and the PSFs simulated here. Try with different numbers of iterations.
In [86]:
arestore "acisf13858_broad_thresh.img[sky=bounds(region(ds9_bkg.reg))]" \
chart_sim.psf arestore.out num=150 clob+
In [87]:
ds9 acisf13858_broad_thresh.img -scale log -pan to 4090.8 4077.6 physical \
-zoom to 4 \
chart_sim.psf -pan to 4090.8 4077.6 physical \
arestore.out -pan to 4090.8 4077.6 physical \
-view colorbar no \
-tile mode column \
-saveimage png ds9_exercise15_ec2.png -quit
display < ds9_exercise15_ec2.png
Create a 3-color image¶
Three color (aka "true color" or "tri-color") images can be useful to help guide analysis by providing visual clues about changes in spectra.
Exercise 16¶
- Obtain the data for OBS_ID 13736. These data will only be used for this exercise.
- Reprocess the data using chandra_repro
- Run fluximage on the reprocessed level 2 event file. Use binsize=8 and bands=csc.
- Load the soft, medium, and hard band flux.img images into ds9. This can be done in the GUI by creating a new RGB frame and then loading the files individually or on the command line
ds9 -rgb -red root_soft_flux.img \
-green root_medium_flux.img \
-blue root_hard_flux.img
- Log scale each of the images.
- Use the Analysis → Smooth on each of the images.
- Save the image in PNG format.
- Exit ds9
- Use dmimg2jpg to create a similar 3-color image.
In [88]:
/bin/rm -rf 13736
download_chandra_obsid 13736 evt1,bpix,flt,mtl,msk,dtf,bias,pbk,asol,stat
Downloading files for ObsId 13736, total size is 319 Mb.
Type Format Size 0........H.........1 Download Time Average Rate
---------------------------------------------------------------------------
evt1 fits 286 Mb #################### 4 s 80481.9 kb/s
asol fits 23 Mb #################### < 1 s 86769.0 kb/s
mtl fits 4 Mb #################### < 1 s 27505.4 kb/s
stat fits 3 Mb #################### < 1 s 58979.4 kb/s
bias fits 489 Kb #################### < 1 s 31375.4 kb/s
bias fits 436 Kb #################### < 1 s 27719.7 kb/s
bias fits 429 Kb #################### < 1 s 29922.9 kb/s
bias fits 427 Kb #################### < 1 s 29488.4 kb/s
bias fits 425 Kb #################### < 1 s 30084.3 kb/s
bpix fits 81 Kb #################### < 1 s 7408.7 kb/s
flt fits 7 Kb #################### < 1 s 692.8 kb/s
msk fits 5 Kb #################### < 1 s 503.0 kb/s
pbk fits 4 Kb #################### < 1 s 410.9 kb/s
Total download size for ObsId 13736 = 319 Mb
Total download time for ObsId 13736 = 4 s
In [89]:
download_obsid_caldb 13736 ./CALDB
download_obsid_caldb
infile = 13736
outdir = ./CALDB
background = no
missing = no
clobber = no
verbose = 1
mode = ql
Retrieving files for CALDB_VER = 4.12.0
Retrieving CALDB index files
Processing infile=13736/secondary/acisf13736_000N002_evt1.fits.gz
Retrieving CALDB data files
Filename: 0------------------1
telD1999-07-23geomN0007.fits .................... (skipped)
telD1999-07-23aimptsN0002.fits .................... (skipped)
telD1999-07-23tdetN0001.fits .................... (skipped)
telD1999-07-23skyN0002.fits .................... (skipped)
telD1999-07-23sgeomN0001.fits .................... (skipped)
hrmaD1996-12-20axeffaN0008.fits .................... (skipped)
hrmaD1996-12-20vignetN0003.fits .................... (skipped)
acisD1997-04-17qeN0006.fits .................... (skipped)
acisD2010-02-01qeuN0007.fits .................... (skipped)
acisD2000-11-28badpixN0004.fits .................... (skipped)
acisD1999-08-13contamN0015.fits .................... (skipped)
acisD1996-11-01gradeN0004.fits .................... (skipped)
acisD2000-01-29grdimgN0001.fits .................... (skipped)
acisD2000-01-29gain_ctiN0008.fits .................... (skipped)
acisD2005-07-01evtspltN0002.fits .................... (skipped)
acisD2010-01-01ctiN0012.fits .................... (skipped)
acisD2012-02-01t_gainN0008.fits ####################
acisD1999-07-22subpixN0001.fits .................... (skipped)
acisD2000-01-29fef_pha_ctiN0004.fits .................... (skipped)
acisD1999-09-16dead_areaN0001.fits .................... (skipped)
hrmaD1996-12-20reefN0001.fits .................... (skipped)
acisD2000-01-29p2_respN0009_105-107.fits.................... (skipped)
acisD2000-01-29p2_respN0009_107-109.fits.................... (skipped)
acisD2000-01-29p2_respN0009_109-111.fits.................... (skipped)
acisD2000-01-29p2_respN0009_111-113.fits.................... (skipped)
acisD2000-01-29p2_respN0009_113-115.fits.................... (skipped)
acisD2000-01-29p2_respN0009_115-117.fits.................... (skipped)
acisD2000-01-29p2_respN0009_117-119.fits.................... (skipped)
acisD2000-01-29p2_respN0009_119-120.fits.................... (skipped)
acisD2002-11-15gtilimN0004.fits .................... (skipped)
In [90]:
chandra_repro 13736 out=
Running chandra_repro
version: 07 April 2025
Processing input directory '/lenin1.real/Junk/Workbook/13736'
No boresight correction update to asol file is needed.
Resetting afterglow status bits in evt1.fits file...
Running the destreak tool on the evt1.fits file...
Running acis_build_badpix and acis_find_afterglow to create a new bad pixel file...
Running acis_process_events to reprocess the evt1.fits file...
Filtering the evt1.fits file by grade and status and time...
Applying the good time intervals from the flt1.fits file...
The new evt2.fits file is: /lenin1.real/Junk/Workbook/13736/repro/acisf13736_repro_evt2.fits
Updating the event file header with chandra_repro HISTORY record
Creating FOV file...
Setting observation-specific bad pixel file in local ardlib.par.
Cleaning up intermediate files
WARNING: Observation-specific bad pixel file set for session use:
/lenin1.real/Junk/Workbook/13736/repro/acisf13736_repro_bpix1.fits
Run 'punlearn ardlib' when analysis of this dataset completed.
The data have been reprocessed.
Start your analysis with the new products in
/lenin1.real/Junk/Workbook/13736/repro
In [91]:
fluximage 13736 out=tricolor band=csc bin=8 clob+
Running fluximage
Version: 04 November 2021
Found 13736/repro/acisf13736_repro_evt2.fits
Using event file 13736/repro/acisf13736_repro_evt2.fits
Using CSC ACIS soft science energy band.
Using CSC ACIS medium science energy band.
Using CSC ACIS hard science energy band.
Aspect solution 13736/repro/pcadf13736_000N001_asol1.fits found.
Bad-pixel file 13736/repro/acisf13736_repro_bpix1.fits found.
Mask file 13736/repro/acisf13736_000N002_msk1.fits found.
The output images will have 303 by 394 pixels, pixel size of 3.936 arcsec,
and cover x=3544.5:5968.5:8,y=2808.5:5960.5:8.
Running tasks in parallel with 4 processors.
Creating 5 aspect histograms for obsid 13736
Creating 15 instrument maps for obsid 13736
Creating 15 exposure maps for obsid 13736
Combining 5 exposure maps for 3 bands (obsid 13736)
Thresholding data for obsid 13736
Exposure-correcting 3 images for obsid 13736
The following files were created:
The clipped counts images are:
tricolor_soft_thresh.img
tricolor_medium_thresh.img
tricolor_hard_thresh.img
The observation FOV is:
tricolor_13736.fov
The clipped exposure maps are:
tricolor_soft_thresh.expmap
tricolor_medium_thresh.expmap
tricolor_hard_thresh.expmap
The exposure-corrected images are:
tricolor_soft_flux.img
tricolor_medium_flux.img
tricolor_hard_flux.img
In [92]:
ds9 -rgb -view colorbar no \
-red tricolor_soft_flux.img -scale log -smooth yes -smooth radius 1 \
-green tricolor_medium_flux.img -scale log -smooth yes -smooth radius 1 \
-blue tricolor_hard_flux.img -scale log -smooth yes -smooth radius 1 \
-saveimage png ds9_exercise_16.png -quit
display < ds9_exercise_16.png
In [93]:
dmimg2jpg \
infile=tricolor_soft_flux.img \
greenfile=tricolor_medium_flux.img \
bluefile=tricolor_hard_flux.img \
outfile=dmimg2jpg_exercise13.jpg \
clob+
# display < dmimg2jpg_exercise13.jpg
Q: Discuss the differences in the output from ds9 and dmimg2jpg. ds9 interactive
Extra Credit¶
Smooth the images with aconvolve before displaying
In [94]:
# Should convolve image and exposure map separately!
aconvolve tricolor_soft_flux.img tricolor_soft_flux_gaus.img "lib:gaus(2,5,1,1,1)" meth=slide clob+
aconvolve tricolor_medium_flux.img tricolor_medium_flux_gaus.img "lib:gaus(2,5,1,1,1)" meth=slide clob+
aconvolve tricolor_hard_flux.img tricolor_hard_flux_gaus.img "lib:gaus(2,5,1,1,1)" meth=slide clob+
In [95]:
dmimg2jpg \
infile=tricolor_soft_flux_gaus.img \
greenfile=tricolor_medium_flux_gaus.img \
bluefile=tricolor_hard_flux_gaus.img \
outfile=dmimg2jpg_exercise13_aconvolve.jpg \
clob+
# eog dmimg2jpg_exercise13_aconvolve.jpg
#display < dmimg2jpg_exercise13_aconvolve.jpg
Smooth the images with csmooth before displaying
In [96]:
# Csmooth really, really wants integer counts, so we smooth the counts image not the flux'ed image.
csmooth tricolor_soft_thresh.img none tricolor_soft_flux_csm.img sigmin=3 sclmax=20 mode=h clob+
csmooth tricolor_medium_thresh.img none tricolor_medium_flux_csm.img sigmin=3 sclmax=20 mode=h clob+
csmooth tricolor_hard_thresh.img none tricolor_hard_flux_csm.img sigmin=3 sclmax=20 mode=h clob+
# WARNING: Remainder will be smoothed on scale of 20.000000 # WARNING: Remainder will be smoothed on scale of 20.000000 # WARNING: Remainder will be smoothed on scale of 20.000000
In [97]:
dmimg2jpg \
infile=tricolor_soft_flux_csm.img \
greenfile=tricolor_medium_flux_csm.img \
bluefile=tricolor_hard_flux_csm.img \
outfile=dmimg2jpg_exercise13_csm.jpg \
clob+
# display < dmimg2jpg_exercise13_csm.jpg
Smooth the images with dmimgadapt before displaying
In [98]:
dmimgadapt tricolor_soft_thresh.img tricolor_soft_flux_cone.img cone min=1 max=20 num=100 radscale=linear counts=25 clob+ mode=h
dmimgadapt tricolor_medium_thresh.img tricolor_medium_flux_cone.img min=1 max=20 num=100 radscale=linear counts=25 clob+ mode=h
dmimgadapt tricolor_hard_thresh.img tricolor_hard_flux_cone.img min=1 max=20 num=100 radscale=linear counts=25 clob+ mode=h
In [99]:
dmimg2jpg \
infile=tricolor_soft_flux_cone.img \
greenfile=tricolor_medium_flux_cone.img \
bluefile=tricolor_hard_flux_cone.img \
outfile=dmimg2jpg_exercise13_cone.jpg \
clob+
# display < dmimg2jpg_exercise13_cone.jpg
Obtain overlapping multi-spectral images of this region from other archives. Use reproject_image to project the other datasets to the same tangent plane as the Chandra dataset. Create tri-color image using red for the lowest energy, and blue for the highest energy datasets.
In [ ]:
Obtain the data for OBS_ID 635 and reprocess it with chandra_repro. Split the event file into 3 equal time intervals by filtering on the TIME column. Run fluximage on each time interval using the broad energy band. Use the 3 time-slices to create a 3 color image.
In [100]:
/bin/rm -rf 635
download_chandra_obsid 635 evt1,bpix,flt,mtl,msk,dtf,bias,pbk,asol,stat
Downloading files for ObsId 635, total size is 204 Mb.
Type Format Size 0........H.........1 Download Time Average Rate
---------------------------------------------------------------------------
evt1 fits 173 Mb #################### 3 s 57166.8 kb/s
asol fits 22 Mb #################### < 1 s 59810.0 kb/s
mtl fits 4 Mb #################### < 1 s 45182.5 kb/s
stat fits 3 Mb #################### < 1 s 63736.0 kb/s
bias fits 436 Kb #################### < 1 s 20444.9 kb/s
bias fits 433 Kb #################### < 1 s 18497.3 kb/s
bias fits 433 Kb #################### < 1 s 24453.6 kb/s
bias fits 428 Kb #################### < 1 s 23479.0 kb/s
bias fits 423 Kb #################### < 1 s 28458.2 kb/s
bpix fits 83 Kb #################### < 1 s 4491.1 kb/s
flt fits 6 Kb #################### < 1 s 482.9 kb/s
msk fits 5 Kb #################### < 1 s 345.7 kb/s
pbk fits 4 Kb #################### < 1 s 341.3 kb/s
Total download size for ObsId 635 = 204 Mb
Total download time for ObsId 635 = 4 s
In [101]:
download_obsid_caldb 635 ./CALDB
download_obsid_caldb
infile = 635
outdir = ./CALDB
background = no
missing = no
clobber = no
verbose = 1
mode = ql
Retrieving files for CALDB_VER = 4.12.0
Retrieving CALDB index files
Processing infile=635/secondary/acisf00635_000N005_evt1.fits.gz
Retrieving CALDB data files
Filename: 0------------------1
telD1999-07-23geomN0007.fits .................... (skipped)
telD1999-07-23aimptsN0002.fits .................... (skipped)
telD1999-07-23tdetN0001.fits .................... (skipped)
telD1999-07-23skyN0002.fits .................... (skipped)
telD1999-07-23sgeomN0001.fits .................... (skipped)
hrmaD1996-12-20axeffaN0008.fits .................... (skipped)
hrmaD1996-12-20vignetN0003.fits .................... (skipped)
acisD1997-04-17qeN0006.fits .................... (skipped)
acisD2000-01-29qeuN0007.fits ####################
acisD2000-01-29badpixN0004.fits ####################
acisD1999-08-13contamN0015.fits .................... (skipped)
acisD1996-11-01gradeN0004.fits .................... (skipped)
acisD2000-01-29grdimgN0001.fits .................... (skipped)
acisD2000-01-29gain_ctiN0008.fits .................... (skipped)
acisD1996-11-01evtspltN0002.fits ####################
acisD2000-01-29ctiN0009.fits ####################
acisD2000-01-29t_gainN0008.fits ####################
acisD1999-07-22subpixN0001.fits .................... (skipped)
acisD2000-01-29fef_pha_ctiN0004.fits .................... (skipped)
acisD1999-09-16dead_areaN0001.fits .................... (skipped)
hrmaD1996-12-20reefN0001.fits .................... (skipped)
acisD2000-01-29p2_respN0009_105-107.fits.................... (skipped)
acisD2000-01-29p2_respN0009_107-109.fits.................... (skipped)
acisD2000-01-29p2_respN0009_109-111.fits.................... (skipped)
acisD2000-01-29p2_respN0009_111-113.fits.................... (skipped)
acisD2000-01-29p2_respN0009_113-115.fits.................... (skipped)
acisD2000-01-29p2_respN0009_115-117.fits.................... (skipped)
acisD2000-01-29p2_respN0009_117-119.fits.................... (skipped)
acisD2000-01-29p2_respN0009_119-120.fits.................... (skipped)
acisD1999-07-22gtilimN0007.fits ####################
In [102]:
chandra_repro 635 out=
Running chandra_repro
version: 07 April 2025
Processing input directory '/lenin1.real/Junk/Workbook/635'
No boresight correction update to asol file is needed.
Resetting afterglow status bits in evt1.fits file...
Running acis_build_badpix and acis_find_afterglow to create a new bad pixel file...
Running acis_process_events to reprocess the evt1.fits file...
Filtering the evt1.fits file by grade and status and time...
Applying the good time intervals from the flt1.fits file...
The new evt2.fits file is: /lenin1.real/Junk/Workbook/635/repro/acisf00635_repro_evt2.fits
Updating the event file header with chandra_repro HISTORY record
Creating FOV file...
Setting observation-specific bad pixel file in local ardlib.par.
Cleaning up intermediate files
WARNING: Observation-specific bad pixel file set for session use:
/lenin1.real/Junk/Workbook/635/repro/acisf00635_repro_bpix1.fits
Run 'punlearn ardlib' when analysis of this dataset completed.
The data have been reprocessed.
Start your analysis with the new products in
/lenin1.real/Junk/Workbook/635/repro
In [103]:
dmstat "635/repro/acisf00635_repro_evt2.fits[cols time]" cen- sig- med-
time[s]
min: 72039166.717 @: 1
max: 72141141.711 @: 451325
mean: 72091543.997
sum: 3.2536788186e+13
good: 451326
null: 0
In [104]:
python -c 'a=72039166;b=72141142;dt3=(b-a)/3.0;print(a+0.0*dt3,a+1.0*dt3,a+2.0*dt3,a+3.0*dt3)'
72039166.0 72073158.0 72107150.0 72141142.0
In [105]:
dmcopy "635/repro/acisf00635_repro_evt2.fits[time=72039166.0:72073158.0]" 635/repro/begin_evt.fits clob+
dmcopy "635/repro/acisf00635_repro_evt2.fits[time=72073158.0:72107150.0]" 635/repro/middle_evt.fits clob+
dmcopy "635/repro/acisf00635_repro_evt2.fits[time=72107150.0:72141142.0]" 635/repro/end_evt.fits clob+
In [106]:
fluximage 635/repro/begin_evt.fits rhooph_begin bin=6 band=broad clob+ verb=0
fluximage 635/repro/middle_evt.fits rhooph_middle bin=6 band=broad clob+ verb=0
fluximage 635/repro/end_evt.fits rhooph_end bin=6 band=broad clob+ verb=0
In [107]:
dmimg2jpg \
infile=rhooph_begin_broad_thresh.img \
greenfile=rhooph_middle_broad_thresh.img \
bluefile=rhooph_end_broad_thresh.img \
outfile=rhooph_time_slice.jpg \
clob+
#display < rhooph_time_slice.jpg
In [108]:
e=635/repro/acisf00635_repro_evt2.fits
f=635/repro/acisf00635_repro_fov1.fits
dmcopy "${e}[sky=region($f),energy=500:7000][bin sky=6,time=::#3]" 635_time_cube.fits clob+
ds9 -frame delete -scale log -rgbcube 635_time_cube.fits -saveimage ds9_ex16_timecube.png -quit
display < ds9_ex16_timecube.png
Extract Light Curve¶
References:
- http://cxc.cfa.harvard.edu/ciao/threads/lightcurve/
- http://cxc.cfa.harvard.edu/ciao/threads/variable/
- http://cxc.cfa.harvard.edu/ciao/why/lightcurve.html
In X-ray astronomy the term "light curve" is used to describe the flux of a source vs time.
Exercise 17¶
- Create a light curve using the dmextract tool with the reprocessed level 2 event file filtered with the source and background regions created in Exercise 11 steps 3 and 6. Use a 100 second bin width: [bin time=::100].
- Repeat 1, using the following bin widths: 200, 500, 1000, 1.
- Plot each of the light curves.
In [109]:
dmextract 13858/repro/acisf13858_repro_evt2.fits"[sky=circle(4090.8,4077.6,10),energy=500:7000][bin time=::100]" \
out=dme_100s.lc op=ltc1 \
bkg=13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_bkg.reg),energy=500:7000]" \
mode=h clob+
In [110]:
cat << EOM > mpl_ex17a.py
import matplotlib.pylab as plt
from pycrates import read_file
tab =read_file("dme_100s.lc[cols dt,net_rate]")
x = tab.get_column("dt").values
y = tab.get_column("net_rate").values
plt.plot(x,y, drawstyle="steps")
plt.xlabel("delta time")
plt.ylabel("net rate c/s")
plt.savefig("out_bin100.png")
EOM
python mpl_ex17a.py
display < out_bin100.png
Qt: Session management error: None of the authentication protocols specified are supported
In [111]:
for tt in 200 500 1000 1
do
dmextract 13858/repro/acisf13858_repro_evt2.fits"[sky=circle(4090.8,4077.6,10),energy=500:7000][bin time=::${tt}]" \
out=dme_${tt}s.lc op=ltc1 \
bkg=13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_bkg.reg),energy=500:7000]" \
mode=h clob+
done
# dmextract (CIAO 4.17.0): WARNING: Binning size=1 smaller than TIMEDEL=3.14104 # dmextract (CIAO 4.17.0): WARNING: Binning size=1 smaller than TIMEDEL=3.14104
In [112]:
cat << EOM > mpl_ex17b.py
import matplotlib.pylab as plt
from pycrates import read_file
plt.subplot(4,1,1)
tab =read_file("dme_200s.lc[cols dt,net_rate]")
x = tab.get_column("dt").values
y = tab.get_column("net_rate").values
plt.plot(x,y, drawstyle="steps")
plt.ylabel("net rate c/s")
plt.subplot(4,1,2)
tab =read_file("dme_500s.lc[cols dt,net_rate]")
x = tab.get_column("dt").values
y = tab.get_column("net_rate").values
plt.plot(x,y, drawstyle="steps")
plt.ylabel("net rate c/s")
plt.subplot(4,1,3)
tab =read_file("dme_1000s.lc[cols dt,net_rate]")
x = tab.get_column("dt").values
y = tab.get_column("net_rate").values
plt.plot(x,y, drawstyle="steps")
plt.ylabel("net rate c/s")
plt.subplot(4,1,4)
tab =read_file("dme_1s.lc[cols dt,net_rate]")
x = tab.get_column("dt").values
y = tab.get_column("net_rate").values
plt.plot(x,y, drawstyle="steps")
plt.ylabel("net rate c/s")
plt.xlabel("delta time")
plt.savefig("out_bin_strip.png")
EOM
python mpl_ex17b.py
display < out_bin_strip.png
Qt: Session management error: None of the authentication protocols specified are supported
In [113]:
dmlist dme_100s.lc cols
-------------------------------------------------------------------------------- Columns for Table Block LIGHTCURVE -------------------------------------------------------------------------------- ColNo Name Unit Type Range 1 TIME_BIN channel Int4 1:179 S/C TT corresponding to mid-exposure 2 TIME_MIN s Real8 456518974.4871399999:456536836.1381000280 Minimum Value in Bin 3 TIME s Real8 456518974.4871399999:456536836.1381000280 S/C TT corresponding to mid-exposure 4 TIME_MAX s Real8 456518974.4871399999:456536836.1381000280 Maximum Value in Bin 5 COUNTS count Int4 - Counts 6 STAT_ERR count Real8 0:+Inf Statistical error 7 AREA pixel**2 Real8 -Inf:+Inf Area of extraction 8 EXPOSURE s Real8 -Inf:+Inf Time per interval 9 COUNT_RATE count/s Real8 0:+Inf Rate 10 COUNT_RATE_ERR count/s Real8 0:+Inf Rate Error 11 BG_COUNTS count Real8 -Inf:+Inf Background Counts 12 BG_ERR count Real8 -Inf:+Inf Error on Background counts 13 BG_AREA pixel**2 Real8 -Inf:+Inf Background Area of Extraction 14 BG_EXPOSURE s Real8 -Inf:+Inf Exposure time of background file 15 BG_RATE count/s Real8 -Inf:+Inf Background Rate 16 NORM_BG_COUNTS count Real8 -Inf:+Inf Background Counts 17 NORM_BG_ERR count Real8 -Inf:+Inf Error on Background counts 18 NET_COUNTS count Real8 -Inf:+Inf Net Counts 19 NET_ERR count Real8 -Inf:+Inf Error on Net Counts 20 NET_RATE count/s Real8 -Inf:+Inf Net Count Rate 21 ERR_RATE count/s Real8 -Inf:+Inf Error Rate -------------------------------------------------------------------------------- World Coord Transforms for Columns in Table Block LIGHTCURVE -------------------------------------------------------------------------------- ColNo Name 2: DT_MIN = +0 [s] +1.0 * (TIME_MIN -456519024.487140) 3: DT = +0 [s] +1.0 * (TIME -456519024.487140) 4: DT_MAX = +0 [s] +1.0 * (TIME_MAX -456519024.487140)
In [114]:
dmlist 13858/repro/acisf13858_repro_evt2.fits blocks,subspace
--------------------------------------------------------------------------------
Dataset: 13858/repro/acisf13858_repro_evt2.fits
--------------------------------------------------------------------------------
Block Name Type Dimensions
--------------------------------------------------------------------------------
Block 1: PRIMARY Null
Block 2: EVENTS Table 16 cols x 134015 rows
Block 3: GTI7 Table 2 cols x 1 rows
Block 4: GTI5 Table 2 cols x 1 rows
Block 5: GTI6 Table 2 cols x 1 rows
Block 6: GTI8 Table 2 cols x 1 rows
--------------------------------------------------------------------------------
Data subspace for block EVENTS: Components: 4 Descriptors: 16
--------------------------------------------------------------------------------
--- Component 1 ---
1 time Real8 TABLE GTI7
456520729.9422979951:456535800.4992064238
2 expno Int4 3:4801
3 ccd_id Int2 7:7
4 node_id Int2 0:3
5 chip [ 1] chipx 1:1024
5 chip [ 2] chipy 1:1024
6 tdet [ 1] tdetx 1:8192
6 tdet [ 2] tdety 1:8192
7 det [ 1] detx 0.50: 8192.50
7 det [ 2] dety 0.50: 8192.50
8 sky [ 1] x 0.50: 8192.50
8 sky [ 2] y 0.50: 8192.50
9 phas Int2 -4096:4095
10 pha Int4 0:36855
11 pha_ro Int4 0:36855
12 energy Real4 0: 1000000.0
13 pi Int4 1:1024
14 fltgrade Int2 0:255
15 grade Int2 0:0,2:2,3:3,4:4,6:6
16 status Bit
--- Component 2 ---
1 time Real8 TABLE GTI5
456520729.9833379984:456535800.5402464867
2 expno Int4 3:4801
3 ccd_id Int2 5:5
4 node_id Int2 0:3
5 chip [ 1] chipx 1:1024
5 chip [ 2] chipy 1:1024
6 tdet [ 1] tdetx 1:8192
6 tdet [ 2] tdety 1:8192
7 det [ 1] detx 0.50: 8192.50
7 det [ 2] dety 0.50: 8192.50
8 sky [ 1] x 0.50: 8192.50
8 sky [ 2] y 0.50: 8192.50
9 phas Int2 -4096:4095
10 pha Int4 0:36855
11 pha_ro Int4 0:36855
12 energy Real4 0: 1000000.0
13 pi Int4 1:1024
14 fltgrade Int2 0:255
15 grade Int2 0:0,2:2,3:3,4:4,6:6
16 status Bit
--- Component 3 ---
1 time Real8 TABLE GTI6
456520730.0243780017:456535800.5812864304
2 expno Int4 3:4801
3 ccd_id Int2 6:6
4 node_id Int2 0:3
5 chip [ 1] chipx 1:1024
5 chip [ 2] chipy 1:1024
6 tdet [ 1] tdetx 1:8192
6 tdet [ 2] tdety 1:8192
7 det [ 1] detx 0.50: 8192.50
7 det [ 2] dety 0.50: 8192.50
8 sky [ 1] x 0.50: 8192.50
8 sky [ 2] y 0.50: 8192.50
9 phas Int2 -4096:4095
10 pha Int4 0:36855
11 pha_ro Int4 0:36855
12 energy Real4 0: 1000000.0
13 pi Int4 1:1024
14 fltgrade Int2 0:255
15 grade Int2 0:0,2:2,3:3,4:4,6:6
16 status Bit
--- Component 4 ---
1 time Real8 TABLE GTI8
456520730.0654180050:456535800.6223264933
2 expno Int4 3:4801
3 ccd_id Int2 8:8
4 node_id Int2 0:3
5 chip [ 1] chipx 1:1024
5 chip [ 2] chipy 1:1024
6 tdet [ 1] tdetx 1:8192
6 tdet [ 2] tdety 1:8192
7 det [ 1] detx 0.50: 8192.50
7 det [ 2] dety 0.50: 8192.50
8 sky [ 1] x 0.50: 8192.50
8 sky [ 2] y 0.50: 8192.50
9 phas Int2 -4096:4095
10 pha Int4 0:36855
11 pha_ro Int4 0:36855
12 energy Real4 0: 1000000.0
13 pi Int4 1:1024
14 fltgrade Int2 0:255
15 grade Int2 0:0,2:2,3:3,4:4,6:6
16 status Bit
Q: Why is there a gap at the beginning and at the end of the light curves with count_rate=0? outside of GTI but inside of TSTART-TSTOP
Q: Is the delta-time, dt, column actually a column in the dmextract output file? No, it's a coord on time column
Q: What can be learned from the light curve with 1 second time resolution? nothing
Q: What energy events were used to make these light curves? Does using an energy filter change things and if so how? **used energy=500:7000, if not then would have used all energies **
Q: ACIS event files contain multiple good time intervals (one for each chip). Which GTI was use? The first one, which in this case is for CCD=7 which is correct
Exercise 18¶
- Create a light curve using the glvary tool using the reprocessed level 2 event file with the source region.
- Repeat 1 using the background region.
- Plot the glvary source and background lightcurves using chips
In [115]:
glvary 13858/repro/acisf13858_repro_evt2.fits"[sky=circle(4090.8,4077.6,10),energy=500:7000]" \
out=glvary_src.prob lc=glvary_src.lc clob+ mode=h
cat << EOM > mpl_ex18a.py
import matplotlib.pylab as plt
from pycrates import read_file
tab =read_file("glvary_src.lc[cols time,count_rate]")
x = tab.get_column("time").values
y = tab.get_column("count_rate").values
plt.plot(x,y, drawstyle="steps")
plt.xlabel("time")
plt.ylabel("count rate c/s")
plt.savefig("glvary_src.png")
EOM
python mpl_ex18a.py
display < glvary_src.png
Qt: Session management error: None of the authentication protocols specified are supported
In [116]:
glvary 13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_bkg.reg),energy=500:7000]" \
out=glvary_bkg.prob lc=glvary_bkg.lc clob+ mode=h
cat << EOM > mpl_ex18b.py
import matplotlib.pylab as plt
from pycrates import read_file
tab =read_file("glvary_bkg.lc[cols time,count_rate]")
x = tab.get_column("time").values
y = tab.get_column("count_rate").values
plt.plot(x,y, drawstyle="steps")
plt.xlabel("time")
plt.ylabel("count rate c/s")
plt.savefig("glvary_bkg.png")
EOM
python mpl_ex18b.py
display < glvary_bkg.png
Qt: Session management error: None of the authentication protocols specified are supported
Q: Describe how the glvary lightcurves differ from the dmextract light curves. it picked binning, no
dtcolumn. No net-rate
Q: Could any variability in the source region be due to variability in the background? No variability, background is not variable either
Extra Credit¶
Use the dither_region tool to compute the time resolved area fraction of the source and separately the background regions. Try using those as input to dmextract and glvary as exposure corrections. Compare to the results without those corrections.
In [117]:
punlearn ardlib
acis_set_ardlib 13858/repro/acisf13858_repro_bpix1.fits abs-
Updated ardlib parameter file: /home/kjg/cxcds_param4/ardlib.par AXAF_ACIS0_BADPIX_FILE -> CALDB AXAF_ACIS1_BADPIX_FILE -> CALDB AXAF_ACIS2_BADPIX_FILE -> CALDB AXAF_ACIS3_BADPIX_FILE -> CALDB AXAF_ACIS4_BADPIX_FILE -> CALDB AXAF_ACIS5_BADPIX_FILE -> 13858/repro/acisf13858_repro_bpix1.fits[BADPIX5] AXAF_ACIS6_BADPIX_FILE -> 13858/repro/acisf13858_repro_bpix1.fits[BADPIX6] AXAF_ACIS7_BADPIX_FILE -> 13858/repro/acisf13858_repro_bpix1.fits[BADPIX7] AXAF_ACIS8_BADPIX_FILE -> 13858/repro/acisf13858_repro_bpix1.fits[BADPIX8] AXAF_ACIS9_BADPIX_FILE -> CALDB
In [118]:
dither_region \
13858/repro/pcadf13858_000N001_asol1.fits \
"circle(4090.8,4077.6,10)" \
dr_src.out \
wcs= 13858/repro/acisf13858_repro_evt2.fits \
mask= 13858/repro/acisf13858_000N002_msk1.fits \
mode=h clob+
In [119]:
dmlist dr_src.out cols
-------------------------------------------------------------------------------- Columns for Table Block AREAFRACTION -------------------------------------------------------------------------------- ColNo Name Unit Type Range 1 TIME s Real8 456520092.2497000098:456535804.9880399704 Time 2 EQPOS(RA,DEC) deg Real8 -360.0: 360.0 Sky Position 3 ROLL deg Real8 -Inf:+Inf Roll angle 4 FRACAREA Real8 -Inf:+Inf Fraction area 5 AREA_CHIP_FRAC[10] Real8(10) -Inf:+Inf Region Area Fraction per chip 6 DELTA_T s Real8 -Inf:+Inf Time 7 CHIP_FRAC_TIME[10] Real8(10) -Inf:+Inf Fraction of ontime per chip 8 STATUS Bit[8] Why fraction < 1
In [120]:
cat << EOM > mpl_ex18_drout.py
import matplotlib.pylab as plt
from pycrates import read_file
tab =read_file("dr_src.out[cols time,fracarea]")
x = tab.get_column("time").values
y = tab.get_column("fracarea").values
plt.plot(x,y, drawstyle="steps")
plt.xlabel("time")
plt.ylabel("fraction of area")
plt.savefig("dr_src.png")
EOM
python mpl_ex18_drout.py
display < dr_src.png
Qt: Session management error: None of the authentication protocols specified are supported
In [121]:
dither_region \
13858/repro/pcadf13858_000N001_asol1.fits \
"region(ds9_bkg.reg)" \
dr_bkg.out \
wcs= 13858/repro/acisf13858_repro_evt2.fits \
mask= 13858/repro/acisf13858_000N002_msk1.fits \
mode=h clob+ maxpix=10000
# dither_region (CIAO 4.17.0): WARNING: Too many pixels, resetting resolution from 1.000000e+00 to 2.000000e+00 # dither_region (CIAO 4.17.0): WARNING: Too many pixels, resetting resolution from 2.000000e+00 to 3.000000e+00
In [122]:
cat << EOM > mpl_ex18_droutb.py
import matplotlib.pylab as plt
from pycrates import read_file
tab =read_file("dr_bkg.out[cols time,fracarea]")
x = tab.get_column("time").values
y = tab.get_column("fracarea").values
plt.plot(x,y, drawstyle="steps")
plt.xlabel("time")
plt.ylabel("fraction of area")
plt.savefig("dr_bkg.png")
EOM
python mpl_ex18_droutb.py
display < dr_bkg.png
Qt: Session management error: None of the authentication protocols specified are supported
The background is dithering across bad columns/pixels, but only a tiny fraction of the area is lost.
In [ ]:
Check for background flares¶
Reference:
Background flares may influence results especially when working with large source regions such as supernovae remnants or clusters.
For an on-axis point source, users should carefully consider whether excluding time due to a background flare would be better than accepting the time in order to obtain higher signal to noise.
Exercise 19¶
In this exercise, students will extract a light curve for a background, source-free region and use it to search for any evidence of a background flare.
- Using dmcoords, dmstat, dmlist, or ds9+dax, determine the CCD_ID where the bright source in the center of the field is located.
- Use dmcopy to exclude the wavdetect sources detected in Exercise 10 step 3 from the reprocessed Level 2 event file.
- Use dmextract to create a light curve using the source free event file, use a CCD_ID filter to select events from chip you found in step 1. Use a bin-width of 259.28 seconds.
- Use the deflare tool to search for background flares on the source chip using method=clean.
- Use dmcopy to apply the deflare output good-time interval, GTI, file to the reprocessed level 2 event file.
In [123]:
punlearn dmcoords
dmcoords 13858/repro/acisf13858_repro_evt2.fits op=sky x=4090.8 y=4077.6 verb=0
pget dmcoords x y ra dec chip_id chipx chipy
4090.8 4077.6 09:14:49.073 +08:53:21.25 7 206.901669103904 517.968804802778
In [124]:
python -c 'from region import *;wav=CXCRegion("acisf13858_wav.src");out=field()-wav;out.write("exclude.reg",fits=True, clobber=True)'
dmcopy 13858/repro/acisf13858_repro_evt2.fits"[sky=region(exclude.reg)]" 13858_evt.holes clob+
In [125]:
ds9 13858_evt.holes -scale log -bin factor 8 -saveimage png ds9_ex19.png -quit
display < ds9_ex19.png
In [126]:
dmextract 13858_evt.holes'[ccd_id=7][bin time=::259.28]' bkg_lc op=ltc1 clob+
In [127]:
cat << EOM > mpl_ex19.py
import matplotlib.pylab as plt
from pycrates import read_file
tab =read_file("bkg_lc[cols dt,count_rate]")
x = tab.get_column("dt").values
y = tab.get_column("count_rate").values
plt.plot(x,y, drawstyle="steps")
plt.xlabel("delta time")
plt.ylabel("count_rate")
plt.savefig("example19.png")
EOM
python mpl_ex19.py
display < example19.png
Qt: Session management error: None of the authentication protocols specified are supported
In [128]:
deflare bkg_lc bkg_lc_flare method=clean save="bkg_lc_flare.png"
#convert bkg_lc_flare.ps -trim bkg_lc_flare.png
display < bkg_lc_flare.png
Parameters used to clean the lightcurve are: script version = 16 May 2023 mean = None clip = 3 scale = 1.2 minfrac = 0.1 outfile = bkg_lc_flare plot = True rateaxis = y color = lime pattern = solid pattern color = red Total number of bins in lightcurve = 69 Max length of one bin = 255.892 s Num. bins with a smaller exp. time = 2 Num. bins with exp. time = 0 = 10 Number of bins with a rate of 0 ct/s = 10 Calculated an initial mean (sigma-clipped) rate of 3.23203 ct/s Lightcurve limits use a scale factor of 1.2 about this mean Filtering lightcurve between rates of 2.69336 and 3.87844 ct/s Number of good time bins = 58 Rate filter: 2.6933609672818477 <= count_rate < 3.8784397928858603 Mean level of filtered lightcurve = 3.232033160738217 ct/s Qt: Session management error: None of the authentication protocols specified are supported libGL error: glx: failed to create dri3 screen libGL error: failed to load driver: nouveau GTI limits calculated using a count-rate filter: (count_rate>2.6933609672818477 && count_rate<3.8784397928858603) The corresponding times are: ((time >= 456520530.16714) && (time < 456531160.64713997)) ; 10.29 ksec, bin 1 ((time >= 456531419.92714) && (time < 456535827.68714)) ; 4.32 ksec, bin 2 Exposure time of lightcurve = 14.87 ks Filtered exposure time = 14.62 ks DTCOR value = 0.986934 Creating GTI file Created: bkg_lc_flare Light curve cleaned using the lc_clean routine. Created: bkg_lc_flare.png
In [129]:
dmcopy 13858/repro/acisf13858_repro_evt2.fits"[@bkg_lc_flare]" 13858_deflare.evt clob+
In [130]:
dmlist 13858/repro/acisf13858_repro_evt2.fits counts
dmlist 13858_deflare.evt counts
134015 131495
Q: What CCD_ID was found in step 1? Which tool(s) were used? chip 7,
dmcoords
Q: Was there any evidence of a background flare in this observation, on this chip? If so, describe the flare (when did it occur within the observation, duration, etc.): **Yes, a very short 1-time-bin duration fare in the 2nd half of the observation **
Extra Credit¶
Try different bin-widths when creating the light curve. Does it affect the results?
In [131]:
dmextract 13858_evt.holes'[ccd_id=7][bin time=::50]' bkg_lc_50s op=ltc1 clob+
deflare bkg_lc_50s bkg_lc_flare_50s method=clean save="bkg_lc_flare_50s.png" plot- verb=0
#convert bkg_lc_flare_50s.ps -trim bkg_lc_flare50s.png
display < bkg_lc_flare_50s.png
Qt: Session management error: None of the authentication protocols specified are supported libGL error: glx: failed to create dri3 screen libGL error: failed to load driver: nouveau
Try different deflare method options.
In [132]:
deflare bkg_lc_50s bkg_lc_sigmaclip_50s method=sigma save="bkg_lc_sigmaclip_50s.png" plot- verb=0 < /dev/null
#convert bkg_lc_sigmaclip_50s.ps -trim bkg_lc_sigmaclip50s.png
display < bkg_lc_sigmaclip_50s.png
Qt: Session management error: None of the authentication protocols specified are supported libGL error: glx: failed to create dri3 screen libGL error: failed to load driver: nouveau
Try looking for flares on the other CCD's used in this observation.
In [133]:
dmextract 13858_evt.holes'[ccd_id=6][bin time=::259]' bkg_6lc op=ltc1 clob+
deflare bkg_6lc bkg_6lc_flare method=clean save="bkg_6lc_flare.png" plot- verb=0
#convert bkg_6lc_flare.ps -trim bkg_6lc_flare.png
display < bkg_6lc_flare.png
Qt: Session management error: None of the authentication protocols specified are supported libGL error: glx: failed to create dri3 screen libGL error: failed to load driver: nouveau
Spectral Analysis¶
The energy resolution and sensitivity of the ACIS detector allows for spectral analysis roughly in the 0.3 to 7.0 keV range.
Exercise 20¶
In this exercise Students will use the specextract tool to extract the spectrum, ARF, and RMF files for the source and background regions.
- Run the specextract script using the reprocessed level 2 event file with the source region from Exercise 11 step 3 and the background region from step 6 with the default parameter settings.
- Run specextract as in step 1 with a different outroot and setting weight=no
- Run specextract as in step 2 with a different outroot and setting correctpsf=yes
- Use dmdiff to compare the spectra (.pi) from steps 2 and 3 to the spectrum in step 1.
- Use dmdiff to compare the auxiliary response file (.arf) from steps 2 and 3 to the ARF in step 1.
- Use dmlist to get the BACKSCAL and EXPOSURE keywords from the source and background spectra.
In [134]:
punlearn specextract
specextract 13858/repro/acisf13858_repro_evt2.fits"[sky=circle(4090.8,4077.6,10)]" \
out=specextract_01 \
bkgfile=13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_bkg.reg)]" \
mode=h clob+
Running specextract Version: 26 October 2023 Extracting Spectra:[███████████████████████████████████████████████████████] 2/2 Generating Aspect Histograms & Weights Maps:[██████████████████████████████] 2/2 Generating ARFs & RMFs:[███████████████████████████████████████████████████] 4/4 Grouping Spectra:[█████████████████████████████████████████████████████████] 2/2 Adding Header Keywords:[███████████████████████████████████████████████████] 2/2
In [135]:
punlearn specextract
specextract 13858/repro/acisf13858_repro_evt2.fits"[sky=circle(4090.8,4077.6,10)]" \
out=specextract_02 \
bkgfile=13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_bkg.reg)]" \
weight=no \
mode=h clob+
Running specextract Version: 26 October 2023 Extracting Spectra:[███████████████████████████████████████████████████████] 2/2 Generating Aspect Histograms & Weights Maps:[██████████████████████████████] 2/2 Generating ARFs & RMFs:[███████████████████████████████████████████████████] 4/4 Grouping Spectra:[█████████████████████████████████████████████████████████] 2/2 Adding Header Keywords:[███████████████████████████████████████████████████] 2/2
In [136]:
punlearn specextract
specextract 13858/repro/acisf13858_repro_evt2.fits"[sky=circle(4090.8,4077.6,10)]" \
out=specextract_03 \
bkgfile=13858/repro/acisf13858_repro_evt2.fits"[sky=region(ds9_bkg.reg)]" \
correctpsf=yes \
mode=h clob+
Running specextract Version: 26 October 2023 WARNING: The 'correct' parameter is ignored when 'weight=yes'. Extracting Spectra:[███████████████████████████████████████████████████████] 2/2 Generating Aspect Histograms & Weights Maps:[██████████████████████████████] 2/2 Generating ARFs & RMFs:[███████████████████████████████████████████████████] 4/4 Grouping Spectra:[█████████████████████████████████████████████████████████] 2/2 Adding Header Keywords:[███████████████████████████████████████████████████] 2/2
In [137]:
/bin/ls specextra*pi
specextract_01.pi specextract_02.pi specextract_03.pi specextract_01_bkg.pi specextract_02_bkg.pi specextract_03_bkg.pi specextract_01_grp.pi specextract_02_grp.pi specextract_03_grp.pi
In [138]:
dmdiff specextract_01.pi specextract_02.pi || echo
Infile 1: specextract_01.pi Infile 2: specextract_02.pi ---------------------------------------------------------------------- Compare Headers ---------------------------------------------------------------------- Compare Key Lists: Compare Keyword Details: # dmdiff (CIAO 4.17.0): WARNING: keyword 'CHECKSUM' comments differ. # dmdiff (CIAO 4.17.0): comment1="HDU checksum updated 2025-05-02T17:12:00" # dmdiff (CIAO 4.17.0): comment2="HDU checksum updated 2025-05-02T17:12:22" # dmdiff (CIAO 4.17.0): WARNING: keyword 'DATASUM' comments differ. # dmdiff (CIAO 4.17.0): comment1="data unit checksum updated 2025-05-02T17:11:40" # dmdiff (CIAO 4.17.0): comment2="data unit checksum updated 2025-05-02T17:12:02" Compare Keyword Values: Keyword: Message: Value(s): Diff: -------- -------------------------------------- ---------------------------------- ------------------------ CHECKSUM Values are not equal WKafYKWfWKafWKWf XE1gaE1dRE1dXE1d DATE Values are not equal 2025-05-02T17:11:40 2025-05-02T17:12:02 BACKFILE Values are not equal specextract_01_bkg.pi specextract_02_bkg.pi RESPFILE Values are not equal specextract_01.rmf specextract_02.rmf ANCRFILE Values are not equal specextract_01.arf specextract_02.arf ---------------------------------------------------------------------- Compare Subspaces ---------------------------------------------------------------------- Compare Subspace Structure: Compare Column Details: Compare Subspace Ranges: component 1 Compare Subspace Ranges: component 2 Compare Subspace Ranges: component 3 Compare Subspace Ranges: component 4 Compare Subspace Regions: component 1 Compare Subspace Regions: component 2 Compare Subspace Regions: component 3 Compare Subspace Regions: component 4 ---------------------------------------------------------------------- Compare Tables ---------------------------------------------------------------------- Compare Table Structure: Block name: SPECTRUM Compare Column Details: Compare Virtual Column Details: Compare Column Data:
Few keyword diffs, data are the same
In [139]:
dmdiff specextract_01.pi specextract_03.pi || echo
Infile 1: specextract_01.pi Infile 2: specextract_03.pi ---------------------------------------------------------------------- Compare Headers ---------------------------------------------------------------------- Compare Key Lists: Compare Keyword Details: # dmdiff (CIAO 4.17.0): WARNING: keyword 'CHECKSUM' comments differ. # dmdiff (CIAO 4.17.0): comment1="HDU checksum updated 2025-05-02T17:12:00" # dmdiff (CIAO 4.17.0): comment2="HDU checksum updated 2025-05-02T17:12:47" # dmdiff (CIAO 4.17.0): WARNING: keyword 'DATASUM' comments differ. # dmdiff (CIAO 4.17.0): comment1="data unit checksum updated 2025-05-02T17:11:40" # dmdiff (CIAO 4.17.0): comment2="data unit checksum updated 2025-05-02T17:12:25" Compare Keyword Values: Keyword: Message: Value(s): Diff: -------- -------------------------------------- ---------------------------------- ------------------------ CHECKSUM Values are not equal WKafYKWfWKafWKWf NBFdN9EcNAEcN9Ec DATE Values are not equal 2025-05-02T17:11:40 2025-05-02T17:12:25 BACKFILE Values are not equal specextract_01_bkg.pi specextract_03_bkg.pi RESPFILE Values are not equal specextract_01.rmf specextract_03.rmf ANCRFILE Values are not equal specextract_01.arf specextract_03.arf ---------------------------------------------------------------------- Compare Subspaces ---------------------------------------------------------------------- Compare Subspace Structure: Compare Column Details: Compare Subspace Ranges: component 1 Compare Subspace Ranges: component 2 Compare Subspace Ranges: component 3 Compare Subspace Ranges: component 4 Compare Subspace Regions: component 1 Compare Subspace Regions: component 2 Compare Subspace Regions: component 3 Compare Subspace Regions: component 4 ---------------------------------------------------------------------- Compare Tables ---------------------------------------------------------------------- Compare Table Structure: Block name: SPECTRUM Compare Column Details: Compare Virtual Column Details: Compare Column Data:
Same, few keys, no data diffs.
In [140]:
dmdiff specextract_01.arf"[#row=1:10]" specextract_02.arf"[#row=1:10]" comment- || echo
Infile 1: specextract_01.arf[#row=1:10] Infile 2: specextract_02.arf[#row=1:10] ---------------------------------------------------------------------- Compare Headers ---------------------------------------------------------------------- Compare Key Lists: # dmdiff (CIAO 4.17.0): WARNING: keyword 'CHANTYPE' not present in infile2. # dmdiff (CIAO 4.17.0): WARNING: keyword 'DETCHANS' not present in infile2. # dmdiff (CIAO 4.17.0): WARNING: keyword 'THRSHLD' not present in infile2. # dmdiff (CIAO 4.17.0): WARNING: keyword 'FEFFILE' not present in infile2. # dmdiff (CIAO 4.17.0): WARNING: keyword 'FILTER' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'RESPFILE' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'ASPFILE' not present in infile1. # dmdiff (CIAO 4.17.0): WARNING: keyword 'FRACEXPO' not present in infile1. Compare Keyword Details: # dmdiff (CIAO 4.17.0): WARNING: keyword 'EXPOSURE' units differ. # dmdiff (CIAO 4.17.0): unit1="s" # dmdiff (CIAO 4.17.0): unit2="secs" Compare Keyword Values: Keyword: Message: Value(s): Diff: -------- -------------------------------------- ---------------------------------- ------------------------ CREATOR Values are not equal mkwarf - Version CIAO 4.17.0 mkarf v0.6.2-0 CHECKSUM Values are not equal ZA6ad83YZA3ad73W PaCTQS9TPYATPY9T DATASUM Values are not equal 3394690086 3396984157 DATE Values are not equal 2025-05-02T17:11:46 2025-05-02T17:12:09 DETNAM Values are not equal ACIS-5678 ACIS-7 EXPOSURE Values are not equal 14873.648987637 14873.64956966 +0.000582023 (+3.91e-06%) HISTNUM Values are not equal 121 72 -49 (-40.5%) ---------------------------------------------------------------------- Compare Subspaces ---------------------------------------------------------------------- Compare Subspace Structure: Compare Column Details: Compare Subspace Ranges: component 1 # dmdiff (CIAO 4.17.0): WARNING: subspace column 'ENERG_LO' range list sizes differ, skipping comparison. # dmdiff (CIAO 4.17.0): infile1=0 # dmdiff (CIAO 4.17.0): infile2=1 # dmdiff (CIAO 4.17.0): WARNING: subspace column 'ENERG_HI' range list sizes differ, skipping comparison. # dmdiff (CIAO 4.17.0): infile1=0 # dmdiff (CIAO 4.17.0): infile2=1 Compare Subspace Regions: component 1 ---------------------------------------------------------------------- Compare Tables ---------------------------------------------------------------------- Compare Table Structure: Block name: SPECRESP Compare Column Details: Compare Virtual Column Details: Compare Column Data: Column: Row: Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ SPECRESP 1 Values are not equal 1.84683752059937 1.8468154668808 -2.20537e-05 (-0.00119%) SPECRESP 2 Values are not equal 5.99422359466553 5.99429273605347 +6.91414e-05 (+0.00115%) SPECRESP 3 Values are not equal 8.91802883148193 8.91824913024902 +0.000220299 (+0.00247%) SPECRESP 4 Values are not equal 12.0882625579834 12.0886631011963 +0.000400543 (+0.00331%) SPECRESP 5 Values are not equal 14.9924325942993 14.9929971694946 +0.000564575 (+0.00377%) SPECRESP 6 Values are not equal 17.5202159881592 17.520923614502 +0.000707626 (+0.00404%) SPECRESP 7 Values are not equal 19.7120895385742 19.7129306793213 +0.000841141 (+0.00427%) SPECRESP 8 Values are not equal 22.3341236114502 22.3351154327393 +0.000991821 (+0.00444%) SPECRESP 9 Values are not equal 26.0509738922119 26.0522270202637 +0.00125313 (+0.00481%) SPECRESP 10 Values are not equal 30.0562591552734 30.0578117370605 +0.00155258 (+0.00517%)
In [141]:
dmstat specextract_01.arf"[cols specresp]" | egrep 'min|max'
min: 0.56783735752 @: 1070 max: 636.17321777 @: 125
In [142]:
dmstat specextract_02.arf"[cols specresp]" | egrep 'min|max'
min: 0.56927418709 @: 1070 max: 636.18145752 @: 125
In [143]:
dmstat specextract_03.arf"[cols specresp]" | egrep 'min|max'
min: 0.56783735752 @: 1070 max: 636.17321777 @: 125
In [144]:
dmlist specextract_01.pi,specextract_01_bkg.pi header,clean | egrep 'BACKSCAL|EXPOSURE'
BACKSCAL 4.681337854E-06 [pixel] Fractional area EXPOSURE 14873.6489876370 [s] Exposure time BACKSCAL 0.0006366619481 [pixel] Fractional area EXPOSURE 14873.6489876370 [s] Exposure time
Q: Record the BACKSCAL and EXPOSURE keywords. Why are these keywords important? Used to scale background to source
Q: Describe and discuss the differences, if any, from step 4. Smaller max by a little bit
Q: Describe and discuss the differences, if any, from step 5: None, warning 'corect=yes is not used when making weighted
Extra Credit¶
Use the rmfimg tool to convert the RMF into an image and display in ds9.
In [145]:
rmfimg specextract_01.rmf specextract_01.rmf.img clob+
In [146]:
ds9 specextract_01.rmf.img -scale log -zoom to fit -view colorbar no -cmap heat \
-saveimage png ds9_exercise20.png -quit
display < ds9_exercise20.png
Sherpa: Modeling and Fitting¶
References
sherpa is the fitting and modeling package developed for CIAO. It allows users to obtain best fit model parameters based on a number of models, fitting methods, and statistics used for uncertainty.
The exercise shown here is an extremely simple sherpa use case. More advanced usage includes user supplied models, linked model parameters, simultaneously fitting multiple datasets, fitting images, etc.
Exercise 21¶
In this exercise Students will begin to use sherpa to model and fit the spectrum of the source.
- Load the ungrouped source spectrum obtained from Exercise 20 step 1 into sherpa
sherpa> load_data("source.pha") - Plot the spectrum using the plot_data command
- Group the data using group_counts(). Select a reasonable number of counts for the dataset. Replot the data.
- Restrict the energy range to a reasonable energy range using the notice() command. Replot the data.
- Set the source model to be an absorbed powerlaw using the set_source() command.
sherpa> set_source(xswabs.abs1*xspowerlaw.p1) - Set the absorption nH to the value listed in the header of the spectrum file.
- Use plot_model() to see the model folded through the instrument responses (ARF and RMF).
- Use show_all() to see all the data and model details.
- subtract() the background. Replot the data.
- fit() the data.
- Use plot_fit() to see the data with the model overplotted.
- Use plot_fit_delchi() to see the data with the error weighted residuals.
- Use conf() to obtain uncertainties on the fitted parameters.
- Use set_method() to try different minimization techniques and repeat the fit().
- Use set_stat() to try different fit statistics and repeat the fit().
- Try different grouping schemes.
- Try different spectral models in the set_source() command.
- Repeat this exercise using the spectrum, ARF and RMF obtain in Exercise 20 steps 2 and 3.
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
Q: Which fit method works best?
Q: Which spectrum fits best?
Q: Which model/param is best?
Q: How does grouping affect the results?
Q: What are the advantages/disadvantages of one method over another?
In [ ]:
Download the notebook.