Simulation for XMM spectra (AO-4)

Last modified: February 26, 2007

This thread uses PINTofALE to generate simulated XMM EPIC-PN, EPIC-MOS, and RGS spectra of a coronal source normalized to a ROSAT/PSPC count rate. The XMM ARF, RMF and RSP files used here can be obtained from the PROPOSAL & PLANNING TOOLS page.

  1. Initialize variables
  2. Define a DEM
  3. Rescale DEM to match ROSAT count rate
  4. Construct spectra
  5. Simulate and plot
  6. Comparison with PIMMS
All the steps described here are repeated in the .par file: XMM_epic_rgs_ao4.par.

0. Initialize the settings

#	start IDL

;	set up the PoA environment
.run initale

;	if INITALE fails, run the script
;	@PoA_constructor

1. Initialize some variables

;	Because some variables will be used repeatedly in the course of this
;	thread, it might be useful to initialize them now. Nailing down these
;	variables now will also allow for something of a checklist for some
;	important quantities/files needed for the task at hand:

;	Local Environment We will set up the pathnames specific to
;	the local installation here

 !LDBDIR = '$CHIANTI'		; Atomic Line Database
		; choose from the predefined '$CHIANTI', '$SPEX', '$APED',
		; or specify the full path name to the line database

 pimmsdir='/soft/pimms/data'	; the full path name to the PIMMS database

;	Source We will characterize our model source by the following:

 !NH           = 3e20			 ; H column density [cm^-2]
 !EDENS        = 1.0e9   		 ; electron number density [cm^-3]
 !ABUND        = getabund('grevesse et al.') ; element abundances
					     ;(SEE: getabund())
 T_components  = [6.1, 6.8, 7.2]		; log(T[K]) components in EM
 EM_components = [6.1d11, 6.1d11, 7.1e11]	; Emission Measure [cm^-3]

;	Observation The observation can be described by:

 EXPTIME       =  50.			; nominal exposure time [ks]
 obs_rate      = 0.1                    ; ROSAT/PSPC count rate [ct/s]
					; (set to 0 or less to be ignored)

 pn_ARF        = 'pn-medium.arf'   	; EPIC-PN ARF filename or 'none'
 mos_ARF       = 'mos1-medium.arf' 	; EPIC-MOS ARF filename or 'none'
 rgs1_ARF      = 'none'			; RGS1 ARF filename or 'none'
 rgs2_ARF      = 'none'			; RGS2 ARF filename or 'none'
 pn_RMF        = 'pn.rmf'	     	; EPIC-PN RMF filename
 mos_RMF       = 'mos1.rmf'        	; EPIC-MOS RMF filename
 rgs1_RMF      = 'RGS1ORDER1.RSP'  	; RGS-1 RMF filename
 rgs2_RMF      = 'RGS2ORDER1.RSP'  	; RGS-2 RMF filename

;NOTE: the standard RGS response matrices include the
;effective area, and hence an ARF should not be specified for them.

;	Analysis Parameters 
;       We may restrict the analysis to the wavelength ranges of interest.
;       Here we set the ROSAT/PSPC passband used to measure obs_rate,
;       the desired XMM output passband, and the wavelength grid
;       for the idealized model spectrum. The output spectra will be        
;       defined on the default response energy grids restricted to 
;       the ranges defined here.  
 emin_rosat = 0.12	; minimum energy in ROSAT/PSPC passband [keV]
 emax_rosat = 2.48	; maximum energy in ROSAT/PSPC passband [keV] 
 emin_xmm = 0.2         ; minimum energy in XMM passband [keV]
 emax_xmm = 8.0 	; maximum energy in XMM passband [keV]
 wmin_rosat = !fundae.kevang/emax_rosat ; convert to [Ang]
 wmax_rosat = !fundae.kevang/emin_rosat ; convert to [Ang]
 wmin_xmm = !fundae.kevang/emax_xmm ; convert to [Ang]
 wmax_xmm = !fundae.kevang/emin_xmm ; convert to [Ang]
 !WMIN = wmin_rosat < wmin_xmm ; minimum wavelength for model [Ang]
 !WMAX = wmax_rosat > wmax_xmm ; maximum wavelength for model [Ang]
 nwbin   =  10000L ; number of bins in model spectrum 

2. Define a DEM

;	A Differential Emission Measure (DEM) is required to estimate the
;	amount of emission at various temperatures. Typically, a 2-temperature
;	model is used. Here we will use PINTofALE's mk_dem(), which constructs
;	a DEM array given a temperature grid and emission measure components.
;	We use as the temperature grid !LOGT. The emission measure components
;	are T_components and EM_components as defined above. We choose 'delta'
;	to treat the EM components as delta functions (see mk_dem() or the
;	Chandra/ACIS example for more options)

!DEM=mk_dem('delta', logT = !LOGT, pardem=T_components, indem=EM_components)

3. Rescale the DEM to match ROSAT count rate

;	We will assume that a ROSAT/PSPC count rate is available,
;	and that the simulation will match this rate.  We will assume
;	a count rate of 0.1 cts/s .  Data from other missions such
;	as ASCA, BeppoSAX, etc. can be dealt with in the same manner.

;	First find and read in the ROSAT/PSPC effective area.
;	You will need to know where your local PIMMS installation
;	is to do this.

rd_pimms_file, rosat_pspc_open, pspc_effar, pspc_wvlar, /wave

;	Make sure that the wavelengths are sorted in increasing order

ae=sort(pspc_wvlar) & pspc_wvlar=pspc_wvlar[ae] & pspc_effar=pspc_effar[ae]

;	The following is similar to the process described in the detailed
;	example thread (see Section 1 and Section 2)

A] Read in line cooling emissivities and calculate line intensities ; Read line cooling emissivities of all possible lines in the ; ROSAT/PSPC wavelength range from the atomic data base. ; NOTE: To avoid multiple reads of the line emissivity database, we ; shall read in the emissivities over the entire range of interest lconf=rd_line(atom,n_e=!EDENS,$ wrange=[MIN(pspc_wvlar)<!WMIN,MAX(pspc_wvlar)>!WMAX],$ dbdir=!LDBDIR,verbose=!VERBOSE,wvl=LWVL,logT=LLOGT,Z=Z,$ ion=ION,jon=JON,fstr=lstr) ; The output of rd_line() will only include level population, ; and not ion balances. We will use fold_ioneq() to fold ion balances. ; NOTE: This step should not be performed if !LDBDIR is set to APED, ; which already includes ion balances and abundances. if strpos(strlowcase(!LDBDIR),'aped',0) lt 0 then lconf=$ fold_ioneq(lconf,Z,JON,chidir=!CHIDIR,$ logT=LLOGT,eqfile=!IONEQF,verbose=!VERBOSE) ; And now calculate line intensities using lineflx(). v_ABUND = !ABUND ; NOTE: If !LDBDIR is set to APED, Anders & Grevesse abundances ; are already included in the emissivities. In such cases, either ; leave out the atomic numbers (Z) in the call to LINEFLX() below, ; or redfine the abundance array to be relative to the APED values, ; e.g., if strpos(strlowcase(!LDBDIR),'aped',0) ge 0 then $ v_ABUND = !ABUND/getabund('anders & grevesse') linint=lineflx(lconf,!LOGT,LWVL,Z,DEM=!DEM,abund=v_ABUND) ;[ph/s]
B] Read in continuum emissivities and calculate continuum intensities ; We can read in continuum emissivities using rd_cont(). ; It is important to note that the output emissivities of rd_cont() ; are in [1e-23 erg cm^3/s/Ang] and not [1e-23 erg cm^3/s] as with rd_line() ; NOTE: To avoid multiple reads of the continuum emissivity database, ; we shall read in the emissivities over the entire range of interest cconf=rd_cont(!CEROOT,n_e=!EDENS,$ wrange=[min(pspc_wvlar)<!WMIN,max(pspc_wvlar)>!WMAX],$ dbdir=!CDBDIR,abund=!ABUND,verbose=!VERBOSE,$ wvl=CWW,logT=ClogT) ; The continuum intensities per angstrom can be calculated again using ; lineflx(). Note that CWW contains the wavelength bin boundaries for ; the emissivity array. CWVL=0.5*(CWW[1:*]+CWW) conint=lineflx(cconf,!LOGT,CWVL,DEM=!DEM) ;[ph/s/Ang] ; Now to get just continuum intensity, we must multiply by an array ; containing the bin widths. If we define this array simply ; with: CDW=CWW[1:*]-CWW, we will get an ugly 'saw-toothed' figure. ; (a side-effect of the way the data-base is constructed) To work ; around this, we can use CWVL, the mid-bin values, and mid2bound(), ; which gives intelligent bin-boundary values given mid-bin values: CWB=mid2bound(CWVL) & CDW=CWB[1:*]-CWB conint=conint*CDW ;[ph/s/Ang]*[Ang]
C] Correct for inter-stellar absorption ; Derive ISM absorptions using ismtau() ltau=ismtau(LWVL,NH=!NH,fH2=!fH2,He1=!He1,HeII=!HeII,$ Fano=Fano,wam=wam,/bam,abund=!ABUND,verbose=!VERBOSE) ctau=ismtau(CWVL,NH=!NH,fH2=!fH2,He1=!He1,HeII=!HeII,$ Fano=Fano,wam=wam,/bam,abund=!ABUND,verbose=!VERBOSE) ltrans=exp(-ltau) & ctrans=exp(-ctau) ; Derive theoretical line fluxes linflx = linint * ltrans ;[ph/s/cm^2] ; Derive theoretical continuum fluxes conflx = conint * ctrans ;[ph/s/cm^2]
D] Bin spectra and fold in effective area ; make input theoretical spectrum grid nwbin_pspc = n_elements(pspc_effar) dwvl=float((max(pspc_wvlar)-min(pspc_wvlar))/nwbin_pspc) wgrid=findgen(nwbin+1L)*dwvl+min(pspc_wvlar) ; Rebin to form theoretical line spectrum using hastrogram() linspc = hastogram(abs(LWVL),wgrid,wts=linflx) ;[ph/s/cm^2/bin] ; Rebin to form theoretical continuum spectrum using rebinw() conspc = rebinw(conflx,CWVL,wgrid,/perbin) ;[ph/s/cm^2/bin] ; Derive predicted flux spectrum. WVLS=0.5*(WGRID[1:*]+WGRID) newEffAr=(interpol(pspc_effar,pspc_wvlar,WVLS) > 0) < (max(pspc_effar)) flxspc = (linspc + conspc) * newEffAr ; Derive predicted counts spectrum. flxspc=flxspc*EXPTIME*1e3 ;[ct/bin] ; Restrict the counts spectrum to the specified ACIS-I range oo = where(wgrid ge wmin_rosat and wgrid le wmax_rosat) flxspc = flxspc(oo)
; Now get the total count rate and renormalize the DEM to the ; observed rate of 0.1 ct/s. pred_rate = total(flxspc/EXPTIME/1e3) ;[ct/s] print,'' if obs_rate gt 0 then $ print,'Rescaling input DEM by a factor '+strtrim(obs_rate/pred_rate,2) print,'' rescale_factor=1.0 if obs_rate gt 0 then rescale_factor = obs_rate/pred_rate !DEM = !DEM * rescale_factor linint = linint * rescale_factor conint = conint * rescale_factor linflx = linflx * rescale_factor conflx = conflx * rescale_factor

4. Construct XMM Spectra

;       To construct the XMM spectra, we use the line and continuum
;	emissivities read in above (LCONF and CCONF), the line and
;	continuum intensities (LININT and CONINT) and fluxes (LINFLX
;	and CONFLX) computed above, and recompute the predicted fluxes
;	using XMM ARFs, and finally compute the observed spectra by
;	convolving with the RMFs.

A] Read in line cooling emissivities and calculate line intensities help,lconf,linint
B] Read in line cooling emissivities and calculate line intensities help,cconf,conint
C] Correct for inter-stellar absorption help,ltau,ltrans,linflx,ctau,ctrans,conflx
D] Bin spectra and fold in effective area EMAX = !fundae.kevang/!WMIN & EMIN = !fundae.kevang/!WMAX dnrg = float((EMAX-EMIN)/nwbin) ; bin size egrid = findgen(nwbin+1L)*dnrg+EMIN ; bin boundaries [keV] emid = 0.5*(egrid[1:*]+egrid) ; mid-bin values [keV] wvls = !fundae.kevang/emid ; [ang] linspc = hastogram(!fundae.kevang/abs(LWVL),egrid,wts=linflx) ; [ph/s/cm^2/bin] conspc = rebinw(conflx,!fundae.kevang/CWB,egrid,/perbin) ; [ph/s/cm^2/bin] ;Read in the effective areas using rdarf() if strlowcase(pn_ARF) eq 'none' then effar_pn=0*emid+1. else $ effar_pn = rdarf(pn_ARF,pn_ARF_str) if strlowcase(mos_ARF) eq 'none' then effar_mos=0*emid+1. else $ effar_mos = rdarf(mos_ARF,mos_ARF_str) if strlowcase(rgs1_ARF) eq 'none' then effar_rgs1=0*emid+1. else $ effar_rgs1 = rdarf(rgs1_ARF,rgs1_ARF_str) if strlowcase(rgs2_ARF) eq 'none' then effar_rgs2=0*emid+1. else $ effar_rgs2 = rdarf(rgs2_ARF,rgs2_ARF_str) ;figure out the wavelength grid for effective areas if n_tags(pn_ARF_str) eq 0 then nrgar_pn=emid else $ nrgar_pn = (0.5*(pn_ARF_str.ELO +pn_ARF_str.EHI)) if n_tags(mos_ARF_str) eq 0 then nrgar_mos=emid else $ nrgar_mos = (0.5*(mos_ARF_str.ELO+mos_ARF_str.EHI)) if n_tags(rgs1_ARF_str) eq 0 then nrgar_rgs1=emid else $ nrgar_rgs1 = (0.5*(rgs1_ARF_str.ELO+rgs1_ARF_str.EHI)) if n_tags(rgs2_ARF_str) eq 0 then nrgar_rgs2=emid else $ nrgar_rgs2 = (0.5*(rgs2_ARF_str.ELO+rgs2_ARF_str.EHI)) ;interpolate to put effective area on binned spectra grids new_effar_pn = (interpol(effar_pn ,nrgar_pn ,EMID) > 0) < (max(effar_pn)) new_effar_mos = (interpol(effar_mos ,nrgar_mos ,EMID) > 0) < (max(effar_mos)) new_effar_rgs1 = (interpol(effar_rgs1,nrgar_rgs1,EMID) > 0) < (max(effar_rgs1)) new_effar_rgs2 = (interpol(effar_rgs2,nrgar_rgs2,EMID) > 0) < (max(effar_rgs2)) ;[ct/s/bin] (if DEM is [cm-5]: [ct/s/cm2/bin]) flxspc_pn = (linspc + conspc) * new_effar_pn flxspc_mos = (linspc + conspc) * new_effar_mos flxspc_rgs1 = (linspc + conspc) * new_effar_rgs1 flxspc_rgs2 = (linspc + conspc) * new_effar_rgs2 ;Derive predicted counts spectrum flxspc_pn = flxspc_pn *EXPTIME*1e3 ;[ct/bin] flxspc_mos = flxspc_mos *EXPTIME*1e3 ;[ct/bin] flxspc_rgs1 = flxspc_rgs1 *EXPTIME*1e3 ;[ct/bin] flxspc_rgs2 = flxspc_rgs2 *EXPTIME*1e3 ;[ct/bin]
E] Convolve with RMF using ; read in RMFs pn_RMF_str=rd_ogip_rmf(pn_RMF) mos_RMF_str=rd_ogip_rmf(mos_RMF) rgs1_RMF_str=rd_ogip_rmf(rgs1_RMF) rgs2_RMF_str=rd_ogip_rmf(rgs2_RMF) conv_rmf, egrid, flxspc_pn, CHAN_pn, CTSPC_pn, pn_RMF_str conv_rmf, egrid, flxspc_mos, CHAN_mos, CTSPC_mos, mos_RMF_str conv_rmf, egrid, flxspc_rgs1, CHAN_rgs1, CTSPC_rgs1, rgs1_RMF_str conv_rmf, egrid, flxspc_rgs2, CHAN_rgs2, CTSPC_rgs2, rgs2_RMF_str
; Get co-added RGS spectrum. CHAN_rgs = CHAN_rgs1 ;[keV] CTSPC_rgs = CTSPC_rgs1 + CTSPC_rgs2 ;[ct/bin] valid grids same flxspc_rgs = flxspc_rgs1+flx_spc_rgs2 ;[ct/bin] (if sperate ARF) ; Restrict the spectra to the specified range oo_rgs = where(CHAN_rgs gt emin_xmm and CHAN_rgs lt emax_xmm) oo_mos = where(CHAN_mos gt emin_xmm and CHAN_mos lt emax_xmm) oo_pn = where(CHAN_pn gt emin_xmm and CHAN_pn lt emax_xmm) CTSPC_rgs1 = CTSPC_rgs1(oo_rgs) & CHAN_rgs1 = CHAN_rgs1(oo_rgs) CTSPC_rgs2 = CTSPC_rgs2(oo_rgs) & CHAN_rgs2 = CHAN_rgs2(oo_rgs) CTSPC_rgs = CTSPC_rgs(oo_rgs) & CHAN_rgs = CHAN_rgs(oo_rgs) CTSPC_mos = CTSPC_mos(oo_mos) & CHAN_mos = CHAN_mos(oo_mos) CTSPC_pn = CTSPC_mos(oo_pn) & CHAN_pn = CHAN_pn(oo_pn) ; Note that the output energy grid of the spectra will be by ; default the energy grid defined by the RMF. The spectrum ; however will only show lines and continuum between the ; selected wavelength ranges.

5. Obtain Poisson deviates and plot

;	The final step is a simulation of counts based on the spectrum
;	predicted above

nbin_pn  = n_elements(CTSPC_pn)  &  CTSIM_pn  = intarr(nbin_pn) 
nbin_mos = n_elements(CTSPC_mos) &  CTSIM_mos  = intarr(nbin_mos) 
nbin_rgs = n_elements(CTSPC_rgs) &  CTSIM_rgs = intarr(nbin_rgs) 

for i=0L,nbin_pn-1L  do if CTSPC_pn[i]  gt 0 then $
	CTSIM_pn[i] =randomu(seed,poisson=CTSPC_pn[i])
for i=0L,nbin_mos-1L do if CTSPC_mos[i] gt 0 then $
for i=0L,nbin_rgs-1L do if CTSPC_rgs[i] gt 0 then $

;	The results of the calculations are plotted below

;	The EPIC-pn model (in white) and simulated counts (in red)
window, 0
plot,   CHAN_pn, CTSPC_pn, title='XMM EPIC-PN MODEL SPECTRUM',/xl,/yl,$
        xtitle='[keV]', ytitle='[ct]', ystyle=1, xstyle=1, $
        xrange=[emin_xmm,emax_xmm], yrange=[0.5,2*max(CTSPC_pn)] & stample
oplot,  CHAN_pn, CTSIM_pn ,color = 2, psym=10,  thick=2   

;	The EPIC-mos model (in white) and simulated counts (in red)
window, 1
plot,   CHAN_mos, CTSPC_mos, title='XMM EPIC-MOS MODEL SPECTRUM',/xl,/yl,$
        xtitle='[keV]', ytitle='[ct]', ystyle=1, xstyle=1, $
        xrange=[emin_xmm,emax_xmm], yrange=[0.5,2*max(CTSPC_mos)] & stample
oplot,  CHAN_mos, CTSIM_mos, color = 2, psym=10, thick=2

;	The RGS model (in white) and simulated counts (in red).  Also shown
;	are the model counts from RGS1 (yellow) and RGS2 (green)
window, 2
plot,   !fundae.kevang/CHAN_rgs, CTSIM_rgs, title='XMM RGS MODEL SPECTRUM',$
        xtitle='['+!AA+']', ytitle='[ct]', ystyle=1, xstyle=1, /nodata,$
oplot,  !fundae.kevang/CHAN_rgs, CTSIM_rgs, color=2, psym=10, thick=2
oplot,  !fundae.kevang/CHAN_rgs1, CTSPC_rgs1, color=3, psym=10, thick=2
oplot,  !fundae.kevang/CHAN_rgs2, CTSPC_rgs2, color=4, psym=10, thick=2

;	As above, but zoomed in to show the low counts detail.
window, 3
plot,   !fundae.kevang/CHAN_rgs, CTSIM_rgs, title='XMM RGS MODEL SPECTRUM',$
        xtitle='['+!AA+']', ytitle='[ct]', ystyle=1, xstyle=1, /nodata,$
oplot,  !fundae.kevang/CHAN_rgs, CTSIM_rgs, color=2, psym=10
oplot,  !fundae.kevang/CHAN_rgs1, CTSPC_rgs1, color=3, psym=10
oplot,  !fundae.kevang/CHAN_rgs2, CTSPC_rgs2, color=4, psym=10

;       We may summarize results as follows:
; Simulated counts spectra are in:
;    CTSIM_pn(chan_pn),CTSIM_mos(CHAN_mos),CTSIM_rgs(CHAN_rgs) [ct/bin]  
; Redistributed counts spectra are in: 
;    CTSPC_pn(chan_pn),CTSPC_mos(CHAN_mos),CTSPC_rgs(CHAN_rgs) [ct/bin]
; Predicted counts spectra are in: 
;    FLXSPC_pn(wvls),FLXSPC_mos(wvls), FLXSPC_rgs(wvls) [ct/bin]
;    Note: Standard RGS response matrices include the effective area, so, unless
;    a seperate rmf and effective area files are used, FLXSPC_rgs will
;    NOT contain predicted counts spectra.  
; Theoretical line fluxes are in: 
;    linspc(lwvl) [ph/s[/cm^2]]
; Theoretical continuum fluxes are in: 
;    conspc(cwvl) [ph/s[/cm^2]]

6. Comparison with PIMMS

;       Results using this thread  are roughly consistent with 
;       PIMMS results. PIMMS returned the following XMM count rates
;       for a ROSAT/PSPC count rate of 0.1 cts/s and a Raymond-Smith 
;       model of temperature 0.5437 keV: 
;             XMM PN MED  :   7.792E-01 cts/s
;             XMM MOS MED :   1.958E-01 cts/s
;             XMM RGS1 O1 :   2.629E-02 cts/s
;             XMM RGS2 O1 :   3.627E-02 cts/s
;       Running this thread with  !ABUND = getabund('Allen'), 
;       T_components = [6.8], and  obs_rate = 0.1 gives: 
;             XMM PN MED  :  0.74  cts/s
;             XMM MOS MED :  0.19  cts/s
;             XMM RGS1 O1 :  0.027 cts/s
;             XMM RGS2 O1 :  0.035 cts/s

pintofale() (LL)