Simulation for Chandra ACIS-I (AO-9)

Last modified: February 26, 2007

This thread uses PINTofALE to generate a simulated Chandra spectrum. The Chandra ARF and RMF used were obtained from the Chandra Proposal Information page and placed in !ARDB. Here, we will generate an ACIS-I spectrum for a star that has a nominal ROSAT/PSPC count-rate of 0.1 ct/s (you must have PIMMS installed locally to do this). The steps involved are:

1. define a DEM
2. rescale to match previous observation (e.g., ROSAT/PSPC)
3. make a simulated ACIS-I spectrum

0. Initialize the settings

#	start IDL
idl

;	set up the PoA environment ;
.run initale

;NOTE:
;	if INITALE fails, run the script
;	@PoA_constructor

;	Change some of the preset defaults to values reasonable
;	for this exercize
!LDBDIR='$SPEX'
	;Warning: APED emissivities include ion balance from
	;Mazzotta et al., and abundances from Anders & Grevesse.
	;Care must be taken to not include them twice:
	if strpos(strlowcase(!LDBDIR),'aped',0) ge 0 then !IONEQF='none'
	if strpos(strlowcase(!LDBDIR),'aped',0) ge 0 then $
		!ABUND = !ABUND/getabund('anders & grevesse')
!NH=3e20

1. Define a DEM

;	A Differential Emission Measure (DEM) is required to estimate
;	the amount of emission at any given temperature.  Typically
;	a 2-temperature model is used, as in the example below:
!DEM = mk_dem('delta', logT=!LOGT, pardem=alog10([7.55e6, 23.4e6]),$
	indem=[5.1d12,6.1d12]*0.05)

;	NOTE: the above is the only option available in the
;	chandra_acis_ao7.par file.

;	In general, stellar coronal plasmas are characterized by a
;	multi-temperature emission measure.  It is possible to
;	define a DEM from scratch interactively, using MK_DEM
;	(left click to define points along a DEM curve, and
;	right click to exit)
;
!DEM = mk_dem('cursor', logT=!LOGT, pardem=[1e10,1e14], $
  xcur=mk_dem_xcur, ycur=mk_dem_ycur)

;	A DEM thus derived may be edited in detail with DEMACS
;
!DEM = demacs(!LOGT, DEM0=!DEM, logT0=!LOGT, group=group, igroup=igroup)

;	more complicated DEMs, (e.g., splines) may also be defined:
;
in_pardem=[3.00000,4.01421,5.35229,5.49464,6.46261,7.06759,7.40923]
in_indem=[500.050,188.068,2.71397,2.16901,159.778,24.0164,0.133006]*1e11
!DEM = mk_dem('spline', logT=!LOGT, indem=in_indem, pardem=in_pardem) > 0

;	for the sake of definiteness, we will use a 2-temperature
;	delta-function emission measure distribution in the following
;	example.

T_components  = [6.9,    7.4]		;log(T[K]) components in the EM
EM_components = [5.1d12, 6.1d12]	;Emission Measure [cm^-3]
!DEM=mk_dem('delta',logT=!LOGT,pardem=T_components,indem=EM_components)

plot,!logT,!DEM,xtitle='log!d10!n(T [K])',ytitle='DEM [cm!u-5!n/logK]'

;NOTE:
;	this is a completely arbitrary DEM(T), and its only saving
;	grace is that it looks sensible.  It definitely would not
;	generate the correct number of ROSAT counts and must be
;	renormalized.

2. Rescale to match, e.g., ROSAT/PSPC

;	We will assume that a ROSAT/PSPC count rate is available,
;	and that the simulation will match this rate.  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.
;
pimmsdir = '/soft/pimms/data'	;<<<change to your local installation
rosat_pspc_open = get_pimms_file('ROSAT','PSPC','OPEN',pdir=pimmsdir)
;
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]

;	some initializations
;
src_DIST=10.	;place the source (say) 10 pc away

;	the DEM needed for the script MAKE_SPECTRUM (which will
;	be called below) must be the DEM defined at the source
;	(distance will be corrected for within the script).  Hence,
;	rescale the currently defined DEM (which has reasonable
;	values for a DEM "at Earth").
;
src_DCOR = 10.D^(alog10(4.*!DPI)+2.*alog10(src_DIST*3.1d18))
!DEM = !DEM * src_DCOR

;	now make a ROSAT/PSPC spectrum using MAKE_SPECTRUM,
;	which is a script and requires that its inputs be
;	in specially named variables and have the correct
;	units.  A more flexible approach is outlined in the
;	example sessions.
;
v_LOGT=!LOGT			;temperature grid
v_DEM=!DEM			;DEM in [cm^-3]
v_DIST=src_DIST			;distance in [pc]
v_EFFAR=pspc_effar > 0		;effective area in [cm^2]
v_WVLAR=pspc_wvlar		;effective area wavelength grid in [Ang]
v_WGRID=mid2bound(pspc_wvlar)	;sorted wavelength bin boundaries in [Ang]
v_EXPTIME=50.			;nominal exposure time [ks] (say)
v_RMF='none'			;name of OGIP-compatible RMF
;
;NOTE:
;	note that convolving with an RMF is unnecessary at this stage
;	because all we need is a total count rate estimate
;
.run make_spectrum

;	get the total count-rate and normalize the DEM to the
;	"observed" rate of 0.1 ct/s
;	
obs_rate = 0.1					;[ct/s]
pred_rate = total(v_flxspc/v_EXPTIME/1e3)	;[ct/s]
print,'normalization factor: ',obs_rate/pred_rate
;normalization factor:     0.0030849559
!DEM = !DEM * obs_rate/pred_rate

3. Make ACIS spectrum

;	first read in the Chandra ACIS-I on-axis ARF
;
acisi_effar=rdarf(!ARDB+'/acisi_aimpt_cy07.arf',arstr)

;	figure out the spectral binning required
;
elo=arstr.ELO & ehi=arstr.EHI & wlo=12.3985/ehi & whi=12.3985/elo

;	now reverse such that wavelengths are in increasing order
;
wlo=reverse(wlo) & whi=reverse(whi) & acisi_effar=reverse(acisi_effar)
wgrid=[wlo,max(whi)] & midwvl=0.5*(wgrid+wgrid[1:*])

;	set up to call MAKE_SPECTRUM again
;
v_LOGT=!LOGT			;temperature grid
v_DEM=!DEM			;DEM in [cm^-3]
v_EFFAR=acisi_effar		;effective area in [cm^2]
v_WVLAR=midwvl			;effective area wavelength grid in [Ang]
v_WGRID=wgrid			;sorted wavelength bin boundaries in [Ang]
v_EXPTIME=50.			;nominal exposure time [ks]
v_RMF=!ARDB+'/acisi_aimpt_cy07.rmf'	;name of OGIP-compatible RMF

;	and make the observed counts spectrum
;
.run make_spectrum

;	simulate an actual observation by finding Poisson deviates
;
nbin=n_elements(v_CTSPC)
sim_spec=intarr(nbin)
for i=0L,nbin-1L do $
  if v_CTSPC[i] gt 0 then sim_spec[i]=randomu(seed,poisson=v_CTSPC[i])

;	plot
;
plot,v_CHAN,v_CTSPC,xtitle='[keV]',ytitle='[ct]',/xl,xr=[0.5,7],$
  /yl,yr=[0.1,300],thick=2
oplot,v_CHAN,sim_spec,psym=10,thick=2,color=3

4.Compare with MARX simulation

;	First note that the above results are consistent with rough
;	PIMMS results.  For a ROSAT/PSPC count rate of 0.1 ct/s,
;	a Raymond-Smith spectrum of temperature 0.54 keV produces
;	approx 0.14 ct/s with Chandra/ACIS-I.  This is within a
;	factor of 2 of the value predicted here:
;
print,total(v_CTSPC)/v_EXPTIME/1e3
;        0.174697957

;	MARX requires a spectrum in an ASCII file of 2 columns,
;	with column 1 containing the energy grid in [keV] and
;	column 2 containing the flux density in [ph/s/cm^2/keV]

;	Now, v_LINSPC[v_WVLS] and v_CONSPC[v_WVLS] are already
;	almost in this form.  They must be converted from
;	[ph/s/cm^2/bin] to [ph/s/cm^2/keV].  A point to note is
;	that v_WVLS was derived within MAKE_SPECTRUM as the
;	mid-bin value of the bin boundaries in v_WGRID.

;	convert from [Ang] to [keV]
;
v_ENERGY=12.3985/v_WVLS
Egrid=12.3985/v_WGRID

;	and get the bin width
;
dE=abs(Egrid[1:*]-Egrid)

;	derive MARX compatible flux density
;
spec_MARX = (v_LINSPC + v_CONSPC) / dE

;	reverse the arrays (so that they are in increasing order of energy)
;
spec_MARX = reverse(spec_MARX)
v_ENERGY  = reverse(v_ENERGY)

;	and print out to file
;
openw,umarx,'spec_marx.dat',/get_lun
for i=0L,n_elements(v_ENERGY)-1L do printf,umarx,v_ENERGY[i],spec_MARX[i]
close,umarx & free_lun,umarx

;	see the MARX manual at http://space.mit.edu/ASC/MARX/
;	for details on how to obtain and run MARX.  Note that updating
;       to MARX 4.2.0 will allow for for Cycle 7 simulations including ACIS
;	contamination (this will be necessary to match the ARF used in this
;	example). 
;        
;	We have run MARX independently on the spectrum written out
;	above, using the recipe prescribed in section 5 of the manual. 
;
;	cp /soft/marx/marx.par .
;	pset marx.par ExposureTime=50000 OutputDir='sim_marx' \
;	  GratingType='NONE' DetectorType='ACIS-I' SourceFlux=-1 \
;	  SpectrumType='FILE' SpectrumFile='spec_marx.dat' \
;	  SourceType='POINT' DitherModel='INTERNAL'
;	marx TSTART=2006.46
;	marx2fits sim_marx/ spec_marx.evt
;	dmextract infile="spec_marx.evt[EVENTS][(x,y)=circle(4018,4141,10)]
;         [bin pi=1:1024:1]" outfile="spec_marx.pi" clobber=yes verbose=2

;	now read in the binned spectrum from the MARX simulation
M_CTSPC=mrdfits('spec_marx.pi',1,hdr)

;	and overplot (in red) on the PINTofALE simulation (green)
;	(note that we do not expect perfect agreement because we have
;	used a very rough ARF in section 3 above.  nevertheless, the
;	agreement is good.)
;
plot,v_CHAN,v_CTSPC,/xl,xr=[0.2,7],/xs,/yl,yr=[0.1,400],/ys,thick=2,$
  xtitle='[keV]',ytitle='[ct]',title='Comparing MARX and PINTofALE'
oplot,v_CHAN,sim_spec,psym=10,color=3
oplot,v_CHAN,M_CTSPC.COUNTS,psym=10,thick=2,color=2