A. Markowitz April 2022 Extracting spectra and light curves for EPIC pn & MOS Imaging mode data ======================================================================= ---------------------------------- Step 1: Initialize XMM-SAS: > setuppython ## my own personal aliased command; you set up python as per your needs > setupxmm ## my own personal alias command; you set up you PATH to point to your ## perl and to your installations of WCStools, XPA, ds9 as necessary python3 -m venv /work/agm/software/pyve; setenv PYTHON /work/agm/software/pyve/bin/python ; setenv PATH ${PATH}:/work/agm/software/pyve/bin setenv PATH ${PATH}:/work/agm/software/wcstools/wcstools-3.9.6:/work/agm/software/xmm/xpa:/work/agm/software/xmm/ds9_82; setenv SAS_PERL /usr/bin/perl; source /work/agm/software/xmm/xmmsas_20210317_1624/setsas.csh ------------------------------------------------ Step 2: Create subdirectories where you can work > cd /work/agm/xmm/0413/0862770701 > mkdir ODF ##### where the raw data files are placed > mkdir proc ##### where we do our pipleine re-processing and extraction > tar xvzf 0862770701.tar.gz This creates a second tar file to be unpacked: > tar xvf 3893_0862770701.TAR Reminder: 3893 = XMM revolution (orbit) number ------------------------------------------------ Step 3: Set up environment variables: setenv SAS_CCFPATH /work/agm/caldb/data/xmm/ccf # <-- Yes, include the "ccf" at the end! setenv SAS_ODF /work/agm/xmm/0413/0862770701/ODF # <-- Currently just points to the ODF directory. It will be # changed later, after the odfingest step, to point to the SUM.SAS file setenv SAS_CCF $SAS_ODF/ccf.cif # <-- This file doesn't exist yet. But it will soon. ------------------------------------------------ Step 3.9: Following section 5 of the XMM ABC guide, everything in the ODF directory must be unzipped: > gunzip *.gz ------------------------------------------------ Steps 4-5: Following Section 5 of the ABC guide "cifbuild" and "odfingest" while in the ODF subdirectory: Step 4: CIFBUILD to create the CCF Index File (ccf.cif): In directory /work/agm/xmm/0413/0862770701/ODF: > cifbuild Step 5a: ODFINGEST (still in the ODF directory) > odfingest # odfingest creates a "....SUM.SAS" file. It is only necessary to run odfongest once on any dataset. If for some reason odfingest must be rerun, you must first delete the earlier SUM.SAS file it produced, or else errors will result. Step 5b: setenv SAS_ODF /full/path/to/ODF/full_name_of_*SUM.SAS, e.g., > setenv SAS_ODF /work/agm/xmm/0413/0862770701/ODF/3893_0862770701_SCX00000SUM.SAS ------------------------------------------------ Step 6a: Following Section 6.1 of the ABC guide: Go to the /proc directory, and run "emproc" & "epproc" to produce calibrated event files for all EPIC cameras: > cd ../proc > emproc ## no arguments needed > epproc ## no arguments needed Step 6b: Rename the events files to something more convenient: mv 3893_0862770701_EMOS1_S001_ImagingEvts.ds M1evli.fits mv 3893_0862770701_EMOS2_S002_ImagingEvts.ds M2evli.fits mv 3893_0862770701_EPN_S003_ImagingEvts.ds PNevli.fits ------------------------------------------------ (Step 7: Extract RGS spectra using rgsproc > rgsproc ## no arguments needed. ------------------------------------------------ Step 8: From the events files, list basic information: observing modes, optical blocking filters, start/stop times: fkeyprint {eventsfile}+1 {keyword} e.g. fkepyrint PNevli.fits+1 DATE-OBS keyword pn MOS1 MOS2 DATE-OBS 2021-03-12T09:53:06 2021-03-12T09:27:25 2021-03-12T09:27:46 DATE-END 2021-03-13T02:29:08 2021-03-13T02:34:31 2021-03-13T02:34:31 TSTART 7.31930054563052E+08 7.31928508736652E+08 7.31928529867299E+08 TSTOP 7.31989816388424E+08 7.31990134539171E+08 7.31990135245893E+08 SUBMODE PrimeFullWindow PrimeFullWindow PrimeFullWindow FILTER Medium Medium Medium TSTART & TSTOP: mission elapsed time = number of seconds since 1998.0 It is not uncommon for the pn to start later and/or end earlier (by 5-30 minutes) than the MOS cameras. -------------------------------------------- Step 9: Following section 6.2 of the ABC guide, we use "evselect" to create first-look images for pn/MOS1/MOS2, with no energy or time filtering. These images will help us identify the central source's location on the detector, and help us plan where to choose regions from which to extract background (avoiding other point sources in the field of view, avoiding hot/bad pixels/rows, etc.) > evselect table=M1evli.fits withimageset=yes imageset=M1firstlookimage.fits xcolumn=X ycolumn=Y imagebinning=imageSize ximagesize=600 yimagesize=600 > evselect table=M2evli.fits withimageset=yes imageset=M2firstlookimage.fits xcolumn=X ycolumn=Y imagebinning=imageSize ximagesize=600 yimagesize=600 > evselect table=PNevli.fits withimageset=yes imageset=PNfirstlookimage.fits xcolumn=X ycolumn=Y imagebinning=imageSize ximagesize=600 yimagesize=600 # withimageset=yes --> we're making an image # imageset = name of output image > ds9 PNfirstlookimage.fits & Reminder: According to NED, the source coordinates are WISEA J041338.08-215204.7 04h13m38.1s -21d52m05s Note the position in detector coordinates of the central source. Draw a 40" region circle around it. Consider off-source circular positions for the background (move the 40" circle around a few arcmins, staying on the same CCD chip, and find a location devoid of background sources, chip edges, etc.). Note the positions in detector coordinates For all three cameras: detector coords are: SRC: 26927, 27711 BGD: 26274, 25405 (There's a nice source-free area a couple arcmins to the SSE of the CLAGN) Note: For pn: sometimes the source ends up a bit too close to the chip edge, and the 40" = 800 pixel radius circle may overlap with the last 1-2 rows near the chip edge. In such cases, it's okay to use e.g., 35" = 700 pixel or 37"=740 pixel radius circles instead -------------------------------------------- Step 10: Create "filtered" event files (which are much smaller in size and faster to work with), by filtering out bad events, and filtering out events marked as photons < 0.2 keV (pn or MOS) or >12 keV (MOS) or > 15 keV (pn), following section 6.3 of the ABC guide: Additionally, filter on pattern, also called event grade (see sect 3.3.11 of the XMM user's handbook, https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/uhb/epic_evgrades.html): Pattern0 = single events (least sensitive to pile-up) Pattern1-4 = double events Pattern5-12 = triple + quadruple events (most sensitive to pile-up) However, for pn, in certain observing modes, Pattern1-4 data yields a calibration error ~0.3-0.4 keV, though pattern 0 is okay there. So we will extract (and fit) pattern0 and pattern1-4 for pn separately. MOS: Because MOS has a lower count rate than pn for a given flux, it's less susceptible to pile-up. Usually it is okay to extract pattern0-12 together unless your source is very bright. Note: I have added "((FLAG & 0xfb0000) == 0)" to the expression here (additional compared to the ABC guide) evselect table=M1evli.fits withfilteredset=yes expression='(PATTERN <= 12)&&(PI in [200:12000])&&#XMMEA_EM&&((FLAG & 0xfb0000) == 0)' \ filteredset=M1evli_filt_p012.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes evselect table=M2evli.fits withfilteredset=yes expression='(PATTERN <= 12)&&(PI in [200:12000])&&#XMMEA_EM&&((FLAG & 0xfb0000) == 0)' \ filteredset=M2evli_filt_p012.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes evselect table=PNevli.fits withfilteredset=yes expression='(PATTERN == 0)&&(PI in [200:15000])&&#XMMEA_EP&&((FLAG & 0xfb0000) == 0)' \ filteredset=PNevli_filt_p0.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes evselect table=PNevli.fits withfilteredset=yes expression='((PATTERN >= 1)&&(PATTERN <= 4))&&(PI in [200:15000])&&#XMMEA_EP&&((FLAG & 0xfb0000) == 0)' \ filteredset=PNevli_filt_p14.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes # withfilteredset=yes ---> we're creating a filtered event file -------------------------------------------- Step 11: Create and display a background light Curve to check for space radiation-induced background flares (Following section 6.4 of the ABC guide) Relatively harder bands are more sensitive to such soft proton background flares. Traditionally, I examined the 10-13 keV light curve, but I've been examining 5-13 keV in recent times as well. Use a 40" (800-pixel) extraction radius centered on the bgd detector coords chosen above. AND/OR extract the light curve for the entire detector (no region filtering) --- I actually do this out of laziness, and this is only advisable for VERY faint sources; you do not want the source intrinsic variability with the flaring background! If there are any bright point sources in the field of view, you can extract over the whole detector while excluding those sources, e.g., using !((X,Y) IN circle(26274,25405,800)) Step 11a: evselect table=PNevli_filt_p0.fits withrateset=yes rateset=PN0_1013bgd_ltcrv.fits filtertype=expression expression='(PI in [10000:13000])&&((X,Y) IN circle(26274,25405,800))' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes evselect table=PNevli_filt_p0.fits withrateset=yes rateset=PN0_0513bgdwhole_ltcrv.fits filtertype=expression expression='(PI in [5000:13000])' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes evselect table=PNevli_filt_p0.fits withrateset=yes rateset=PN0_0513bgd_ltcrv.fits filtertype=expression expression='(PI in [5000:13000])&&((X,Y) IN circle(26274,25405,800))' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes evselect table=PNevli_filt_p0.fits withrateset=yes rateset=PN0_1013bgdwhole_ltcrv.fits filtertype=expression expression='(PI in [10000:13000])' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes evselect table=PNevli_filt_p14.fits withrateset=yes rateset=PN14_1013bgd_ltcrv.fits filtertype=expression expression='(PI in [10000:13000])&&((X,Y) IN circle(26274,25405,800))' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes evselect table=M2evli_filt_p012.fits withrateset=yes rateset=M2_1013bgd_ltcrv.fits filtertype=expression expression='(PI in [10000:13000])&&((X,Y) IN circle(26274,25405,800))' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes Step 11b: fv PN0_1013bgd_ltcrv.fits & Row "RATE" --> Click "Plot" --> Click "TIME" then "X", and click "RATE" then "Y" --> Plot ## 10-13 keV, whole chip, pn0: The bkgd spikes are evident. Clip at 0.7 ct/s ## 10-13 keV, bkgd circle only: PN0: There is a small, narrow spike about a third of the way into the observation. It peaks at ~0.08 - 0.10 ct/s. A couple of small spikes ~0.06 ct/s, very minor ## 10-13 keV, bkgd circle only: PN14: There is a small, narrow spike about a third of the way into the observation. It peaks at ~0.08 - 0.10 ct/s. A couple of small spikes ~0.06 ct/s, very minor ## 10-13 keV, bkgd circle only: MOS2: The spike is barely evident. Peaks at 0.06 ct/s in one 50-sec bin These spikes are very minor compared to other, much larger background spikes I have seen in other observations. 0.05 ct/s seems to be a reasonable threshold for excluding data. Note: FV is not good at displaying all the decimals in the X-axis (time units are MET = mission elapsed time since 1998.0), so it can be difficult to judge the exact start/stop times of flares. fdump PN0_1013bgd_ltcrv.fits+1 outfile=PN0_1013bgd_ltcurv.txt col="TIME,RATE" cp PN0_1013bgd_ltcurv.txt jun emacs jun & ## edit out header cat jun | awk '{print $3}' | sort -g > jun2 Obs Start is MET 731 930 079 Flares occur MET 731 951 079 -- MET 731 951 929 Obs end is MET 731 989 829 Duration of observation is 59750 sec. Number of bins with background flares at 0.08 ct/s or greater: 4 bins / 1196 total bins Number of bins with background flares at 0.06 ct/s or greater: 12 bins / 1196 total bins So I think we can use the pn0 bkgd lightcurve to filter out 1.0 % of the data without worries. -------------------------- Step 12: Applying Time Filters and creating a second filtered set of events files (section 6.5 of the ABC guide) Step 12a: Let's create a new Good Time Interval (GTI) file with the task tabgtigen, filtering on the RATE in the 10-13 keV pn pattern0 whole-detector bkgd light curve. Using 10-13 keV, whole detector: tabgtigen table=PN0_1013bgdwhole_ltcrv.fits gtiset=goodgtisetwhole.fits timecolumn=TIME expression='(RATE<=0.70)' Step 12b: Use evselect to apply goodgtisetwhole.fits to all once-filtered event files to create twice-filtered event files: evselect table=M1evli_filt_p012.fits withfilteredset=yes expression='GTI(goodgtisetwhole.fits,TIME)' filteredset=M1evli_filt2_p012.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes evselect table=M2evli_filt_p012.fits withfilteredset=yes expression='GTI(goodgtisetwhole.fits,TIME)' filteredset=M2evli_filt2_p012.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes evselect table=PNevli_filt_p0.fits withfilteredset=yes expression='GTI(goodgtisetwhole.fits,TIME)' filteredset=PNevli_filt2_p0.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes evselect table=PNevli_filt_p14.fits withfilteredset=yes expression='GTI(goodgtisetwhole.fits,TIME)' filteredset=PNevli_filt2_p14.fits filtertype=expression keepfilteroutput=yes updateexposure=yes filterexposure=yes Step 12c: Extract a new 10-13 keV (or 5-13 keV) pn background light curve to verify that these steps worked: (same as step 11a, except the table is now "filt2" not "filt", and the output file name is different: evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_1013bgdwhole_ltcrv2.fits filtertype=expression expression='(PI in [10000:13000])' maketimecolumn=yes timecolumn=TIME timebinsize=50 makeratecolumn=yes fv PN0_1013bgdwhole_ltcrv2.fits & ----------------------- Step 12.5: Use edetect_chain (Sect 6.6 of the ABC guide) to look for point sources and verify that your background region is devoid of point sources. (One needs an events file that has been filtered to be free of background flares, bad pixels, etc., before doing this) Usually, you can verify with ds9 that the region you choose for background extraction is consistent with the rest of the chip. So I usually skip this step. But if you want something more autoamted, then you can do this: 12.5a) First, make the attitude file: atthkgen atthkset=attitude.fits timestep=1 12.5b) make the soft, total, and hard X-ray images with evselect; extract from PNevli_filt2_p0.fits, M1evli_filt2_p012.fits, M2evli_filt2_p012.fits evselect table=M1evli_filt2_p012.fits withimageset=yes imageset=mos1-s.fits imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 filtertype=expression expression='(FLAG == 0)&&(PI in [300:2000])' evselect table=M1evli_filt2_p012.fits withimageset=yes imageset=mos1-h.fits imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 filtertype=expression expression='(FLAG == 0)&&(PI in [2000:10000])' evselect table=M1evli_filt2_p012.fits withimageset=yes imageset=mos1-t.fits imagebinning=binSize xcolumn=X ximagebinsize=22 ycolumn=Y yimagebinsize=22 filtertype=expression expression='(FLAG == 0)&&(PI in [300:10000])' 12.5c) Run edetect_chain: edetect_chain imagesets='mos1-s.fits mos1-h.fits' \ eventsets='M1evli_filt2_p012.fits' attitudeset=attitude.fits \ pimin='300 2000' pimax='2000 10000' likemin=10 witheexpmap=yes \ ecf='1.85 0.47' eboxl_list=eboxlist_l.fits \ eboxm_list=eboxlist_m.fits eml_list=emllist.fits esp_withootset=no Energy conversion factors (ECFs) convert the source count rates into fluxes. They have units of 10^11 cts cm2/erg Gamma=1.9, NH=1e10 --> 0.3-2: ECF = 1.85 ct/s/(1e-11) 2-10: ECF = 0.472 ct/s/(1e-11) 12.5d) display the results with srcdisplay: srcdisplay boxlistset=emllist.fits imageset=mos1-t.fits regionfile=regionfile.txt sourceradius=0.01 withregionfile=yes This launches ds9; save the sources into a regionfile, and overplot that with the pn image (ds9 PNfirstlookimage.fits & .....) ---------------------------- Step 13: Extract the time-averaged source & background spectra with evselect (sect. 6.7 of the ABC Guide) Note: the parameter "specchannelmax" should be set to 11999 for the MOS, or 20479 for the PN Step 13a: Source extraction (point source): 40" = 800 pixel radius, centered on the source (detector coords = 26927, 27711): evselect table='PNevli_filt2_p0.fits' energycolumn='PI' withfilteredset=yes filteredset='PNevli_filt3_p0.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26927,27711,700))' withspectrumset=yes spectrumset='PN0_SRC.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 evselect table='PNevli_filt2_p14.fits' energycolumn='PI' withfilteredset=yes filteredset='PNevli_filt3_p14.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26927,27711,700))' withspectrumset=yes spectrumset='PN14_SRC.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 evselect table='M1evli_filt2_p012.fits' energycolumn='PI' withfilteredset=yes filteredset='M1evli_filt3_p012.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26927,27711,800))' withspectrumset=yes spectrumset='M1012_SRC.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 evselect table='M2evli_filt2_p012.fits' energycolumn='PI' withfilteredset=yes filteredset='M2evli_filt3_p012.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26927,27711,800))' withspectrumset=yes spectrumset='M2012_SRC.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 filteredset = an even-more filtered output events file that you will need to check for pile-up later <--- need to write both SRC and BKGD filtered-events files spectralbinsize = size of bins, in eV ## In some of my old notes, I see spectralbinsize=5 for pn and spectralbinsize=15 for MOS, but the ABC Guide's examples have ## spectralbinsize=5 for mos. So I guess we should use spectralbinsize=5 for both instruments. Step 13b: Extraction of background spectra: Same as Step 13b, but now centered on the source-free region (detector coords = 26274,25405): evselect table='PNevli_filt2_p0.fits' energycolumn='PI' withfilteredset=yes filteredset='PNevli_bgd3_p0.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26274,25405,700))' withspectrumset=yes spectrumset='PN0_BGD.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 evselect table='PNevli_filt2_p14.fits' energycolumn='PI' withfilteredset=yes filteredset='PNevli_bgd3_p14.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26274,25405,700))' withspectrumset=yes spectrumset='PN14_BGD.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=20479 evselect table='M1evli_filt2_p012.fits' energycolumn='PI' withfilteredset=yes filteredset='M1evli_bgd3_p012.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26274,25405,800))' withspectrumset=yes spectrumset='M1012_BGD.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 evselect table='M2evli_filt2_p012.fits' energycolumn='PI' withfilteredset=yes filteredset='M2evli_bgd3_p012.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE(26274,25405,800))' withspectrumset=yes spectrumset='M2012_BGD.ds' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 (Steps 13a & 3b: You can apply TIME or RATE filters here for time- or flux-resolved spectral extraction!) ---------------------------- Step 14: Check for Pile Up using "epatplot" (section 6.8 of ABC guide): epatplot set=PNevli_filt3_p0.fits plotfile=PN0_epat.ps useplotfile=yes withbackgroundset=yes backgroundset=PNevli_bgd3_p0.fits epatplot set=PNevli_filt3_p14.fits plotfile=PN14_epat.ps useplotfile=yes withbackgroundset=yes backgroundset=PNevli_bgd3_p14.fits epatplot set=M1evli_filt3_p012.fits plotfile=M1012_epat.ps useplotfile=yes withbackgroundset=yes backgroundset=M1evli_bgd3_p012.fits epatplot set=M2evli_filt3_p012.fits plotfile=M2012_epat.ps useplotfile=yes withbackgroundset=yes backgroundset=M2evli_bgd3_p012.fits Look for differences between the predicted and observed pattern distributions, especially towards harder energies (e.g. reduction in double events and an excess in single events). In the case of J0413-21: this source is so faint that there are insufficient stats for epatplot in the pn .... so the source is clearly not piled up. In MOS: no deviation between predicted and observed pattern distributions. If there is pileup: see Section 6.9 of the ABC guide for tips on how to mitigate: re-extract source but ignore central peak, where pileup is most likely to occur; use only single events, etc. --------------------- Step 14.5: Add the background file to the source spectrum via the BACKFILE keyword: fparkey PN0_BGD.ds PN0_SRC.ds+1 BACKFILE fparkey PN14_BGD.ds PN14_SRC.ds+1 BACKFILE fparkey M1012_BGD.ds M1012_SRC.ds+1 BACKFILE fparkey M2012_BGD.ds M2012_SRC.ds+1 BACKFILE Note: Do this step _before_ grppha; grppha will copy these keypwords from the ungrouped file to the grouped file. If you do grppha first, then you muts manually add the BACKFILE keywords to the grouped file, too. --------------------- Step 15: Create the Photon Redistribution Matrix (RMF) and Ancillary File (ARF) (setting setbackscale=yes) Step 15a: rmfgen: rmfgen rmfset=PN0.rmf spectrumset=PN0_SRC.ds rmfgen rmfset=PN14.rmf spectrumset=PN14_SRC.ds rmfgen rmfset=M1012.rmf spectrumset=M1012_SRC.ds rmfgen rmfset=M2012.rmf spectrumset=M2012_SRC.ds (rmfgen only takes a couple minutes) Step 15b: arfgen: arfgen arfset=PN0.arf spectrumset=PN0_SRC.ds withrmfset=yes rmfset=PN0.rmf withbadpixcorr=yes badpixlocation=PNevli_filt2_p0.fits setbackscale=yes arfgen arfset=PN14.arf spectrumset=PN14_SRC.ds withrmfset=yes rmfset=PN14.rmf withbadpixcorr=yes badpixlocation=PNevli_filt2_p14.fits setbackscale=yes arfgen arfset=M1012.arf spectrumset=M1012_SRC.ds withrmfset=yes rmfset=M1012.rmf withbadpixcorr=yes badpixlocation=M1evli_filt2_p012.fits setbackscale=yes arfgen arfset=M2012.arf spectrumset=M2012_SRC.ds withrmfset=yes rmfset=M2012.rmf withbadpixcorr=yes badpixlocation=M2evli_filt2_p012.fits setbackscale=yes (arfgen only takes about a minute to run) badpixlocation = tells arfgen "see this event file for bad pix info"; should be the same as the event file from which the source spectrum was extracted in Step 13a. Step 15c: Add the ARF and RMF files to the source spectra header keywords: fparkey PN0.arf PN0_SRC.ds+1 ANCRFILE fparkey PN14.arf PN14_SRC.ds+1 ANCRFILE fparkey M1012.arf M1012_SRC.ds+1 ANCRFILE fparkey M2012.arf M2012_SRC.ds+1 ANCRFILE fparkey PN0.rmf PN0_SRC.ds+1 RESPFILE fparkey PN14.rmf PN14_SRC.ds+1 RESPFILE fparkey M1012.rmf M1012_SRC.ds+1 RESPFILE fparkey M2012.rmf M2012_SRC.ds+1 RESPFILE Note, Spring 2022: We were seeing some source spectra that, after loading them into Xspec, yielding negative count rates when the background was subtracted. There was an error traced to this step. The solution was to run the "backscale" command manually: backscale spectrumset=PN0_SRC.ds badpixlocation=PNevli_filt2_p0.fits backscale spectrumset=PN0_BGD.ds badpixlocation=PNevli_filt2_p0.fits backscale spectrumset=PN14_SRC.ds badpixlocation=PNevli_filt2_p14.fits backscale spectrumset=PN14_BGD.ds badpixlocation=PNevli_filt2_p14.fits backscale spectrumset=M1012_SRC.ds badpixlocation=M1evli_filt2_p012.fits backscale spectrumset=M1012_BGD.ds badpixlocation=M1evli_filt2_p012.fits backscale spectrumset=M2012_SRC.ds badpixlocation=M2evli_filt2_p012.fits backscale spectrumset=M2012_BGD.ds badpixlocation=M2evli_filt2_p012.fits -------------- Step 16a: Group spectra to a minimum of 15-20 cts/bin to ensure adequate Chi^2 statistics: grppha infile=PN0_SRC.ds outfile=PN0_SRCg.ds comm="group min 20" tempc="exit" clobber=yes grppha infile=PN14_SRC.ds outfile=PN14_SRCg.ds comm="group min 20" tempc="exit" clobber=yes grppha infile=M1012_SRC.ds outfile=M1012_SRCg.ds comm="group min 20" tempc="exit" clobber=yes grppha infile=M2012_SRC.ds outfile=M2012_SRCg.ds comm="group min 20" tempc="exit" clobber=yes Step 16b: Add the ARF and RMF files to the background spectra header keywords, in case we want to load the backgrounds into xspec, too: fparkey PN0.arf PN0_BGD.ds+1 ANCRFILE fparkey PN14.arf PN14_BGD.ds+1 ANCRFILE fparkey M1012.arf M1012_BGD.ds+1 ANCRFILE fparkey M2012.arf M2012_BGD.ds+1 ANCRFILE fparkey PN0.rmf PN0_BGD.ds+1 RESPFILE fparkey PN14.rmf PN14_BGD.ds+1 RESPFILE fparkey M1012.rmf M1012_BGD.ds+1 RESPFILE fparkey M2012.rmf M2012_BGD.ds+1 RESPFILE Step 16c: Make a note of the "good exposure times after screening": foreach fn (PN0_SRCg.ds PN14_SRCg.ds M1012_SRCg.ds M2012_SRCg.ds) fkeyprint $fn+1 EXPOSURE end pn: 40.8 ksec MOS1 & 2: both 46.5 ksec -------------------------------------------------- Step 17: Load all four source spectral into xspec to verify that files got made okay! xspec setplot ene cpd /xs data 1:1 PN0_SRCg.ds ; ign 1: **-0.22 10.0-** data 2:2 PN14_SRCg.ds ; ign 2: **-0.50 10.0-** ## PN14: Ignore only below 0.20 keV, then plot: you see the calibration issue data 3:3 M1012_SRCg.ds ; ign 3: **-0.20 10.0-** data 4:4 M2012_SRCg.ds ; ign 4: **-0.20 10.0-** plot ld plot eff (verify that arf generation has yielded reasonable values) fit a very quick simple model (e.g. plain power-law) and verify that flux is reasonable (also to verify arf generation) -------------------------- Step 18: Preliminary light curve extraction with evselect (energy range = full band for starters) Let's get a first look at the broadband (0.5-10 keV) continuum variability properties; later, we'll re-extract sub-band light curves corresponding to different spectral components dominating. Note: PN0 has roughly 3x the count rate as PN14 or either MOS, so typically I extract light curves from PN0 only. # reminder: Detector coords for SRC: 26927, 27711 BGD: 26274, 25405 # 1-8 keV Src & Bkgd: evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_HXsrc_ltcrv.fits filtertype=expression expression='(PI in [1000:8000])&&((X,Y) IN circle(26927,27711,700))' maketimecolumn=yes timecolumn=TIME timebinsize=500 makeratecolumn=yes evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_HXbgd_ltcrv.fits filtertype=expression expression='(PI in [1000:8000])&&((X,Y) IN circle(26274,25405,700))' maketimecolumn=yes timecolumn=TIME timebinsize=500 makeratecolumn=yes # 0.3-1.0 keV src & bkgd: evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_SXsrc_ltcrv.fits filtertype=expression expression='(PI in [300:1000])&&((X,Y) IN circle(26927,27711,700))' maketimecolumn=yes timecolumn=TIME timebinsize=500 makeratecolumn=yes evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_SXbgd_ltcrv.fits filtertype=expression expression='(PI in [300:1000])&&((X,Y) IN circle(26274,25405,700))' maketimecolumn=yes timecolumn=TIME timebinsize=500 makeratecolumn=yes # 0.3-8.0 keV src & bkgd: evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_TXsrc_ltcrv.fits filtertype=expression expression='(PI in [300:8000])&&((X,Y) IN circle(26927,27711,700))' maketimecolumn=yes timecolumn=TIME timebinsize=500 makeratecolumn=yes evselect table=PNevli_filt2_p0.fits withrateset=yes rateset=PN0_TXbgd_ltcrv.fits filtertype=expression expression='(PI in [300:8000])&&((X,Y) IN circle(26274,25405,700))' maketimecolumn=yes timecolumn=TIME timebinsize=500 makeratecolumn=yes fv PN0_TXsrc_ltcrv.fits & fv PN0_SXsrc_ltcrv.fits & fv PN0_HXsrc_ltcrv.fits & Note: These are "first-look" light curves only. For "production-quality" light curves, we must apply some corrections via "epiclccorr" later! https://heasarc.gsfc.nasa.gov/docs/xmm/sas/help/epiclccorr/node3.html Also, evselect cannot do the background subtraction; we can write out the files to ASCII and use fortran (or python or C or whatever) to do the subtraction. Note: There are sometimes brief tens or hundreds of seconds telemetry drop-outs (count rate goes to zero) -- you must expunged those from the light curve data.