Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://hea-www.harvard.edu/PINTofALE/doc/EXAMPLE_Chandra.html
Дата изменения: Unknown Дата индексирования: Tue Oct 2 03:24:45 2012 Кодировка: Поисковые слова: comet |
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:
# 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
; 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.
; 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
; 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
; 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![]()