Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/PINTofALE/pro/chanarf.pro Дата изменения: Wed Sep 27 21:21:37 2006 Дата индексирования: Tue Oct 2 00:34:10 2012 Кодировка:

function chanarf,rmfstr,effar,spec=spec,chnrg=chnrg,verbose=verbose, _extra=e
;+
;function chanarf
; compute the average effective area over each channel, as a
; weighted average of the effective area from all the energies
; that contribute to that channel.
;
; makes a difference mainly for low and medium spectral
; resolution data. for high-res grating data, not so much.
;
;syntax
; carf=chanarf(rmfstr,effar,spec=spec)
;
;parameters
; rmfstr [INPUT; required] response matrix structure, in the
; same format as returned from RD_OGIP_RMF()
; effar [INPUT] effective area as a function of energy [cm^2]
; * assumed to be on the same grid as the input energies
; of RMFSTR, i.e., RMFSTR.ELO and RMFSTR.EHI
; * ignored if size doesn't match RMFSTR
;
;keywords
; spec [INPUT] a spectrum to further weight the contribution
; of EFFAR
; * assumed to be on the same grid as RMFSTR.ELO and RMFSTR.EHI
; * ignored if size doesn't match EFFAR
; * ideally should be in units of [ph/keV/...]
; chnrg [OUTPUT] the average energy of a photon in this channel
; * it should be similar to 0.5*(RMFSTR.EMN+RMFSTR.EMX), and
; serves as a check on the gain calibration
; verbose [INPUT] controls chatter
; _extra [JUNK] here only to prevent crashing
;
;history
; vinay kashyap (MMVI.IX)
;-

; usage
ok='ok' & np=n_params() & nr=n_elements(rmfstr) & nrt=n_tags(rmfstr)
if np eq 0 then ok='Insufficient parameters' else $
if nr eq 0 then ok='RMFSTR is undefined' else $
if nrt eq 0 then ok='RMFSTR must be a structure; see RD_OGIP_RMF()'
if ok ne 'ok' then begin
print,'Usage: carf=chanarf(rmfstr,effar,spec=spec,chnrg=chnrg,$'
print,' verbose=verbose)'
print,' compute the average effective area over each channel'
if np ne 0 then message,ok,/informational
return,-1L
endif

; inputs
vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1
;
ok='ok' & rtags=tag_names(rmfstr)
iNNRG=where(strpos(rtags,'NNRG') ge 0)
if iNNRG[0] ne -1 then NNRG=rmfstr.NNRG else ok='RMFSTR.NNRG not found'
iELO=where(strpos(rtags,'ELO') ge 0)
if iELO[0] ne -1 then ELO=rmfstr.ELO else ok='RMFSTR.ELO not found'
iEHI=where(strpos(rtags,'EHI') ge 0)
if iEHI[0] ne -1 then EHI=rmfstr.EHI else ok='RMFSTR.EHI not found'
iNCHAN=where(strpos(rtags,'NCHAN') ge 0)
if iNCHAN[0] ne -1 then NCHAN=rmfstr.NCHAN else ok='RMFSTR.NCHAN not found'
iEMN=where(strpos(rtags,'EMN') ge 0)
if iEMN[0] ne -1 then EMN=rmfstr.EMN else ok='RMFSTR.EMN not found'
iEMX=where(strpos(rtags,'EMX') ge 0)
if iEMX[0] ne -1 then EMX=rmfstr.EMX else ok='RMFSTR.EMX not found'
iN_GRP=where(strpos(rtags,'N_GRP') ge 0)
if iN_GRP[0] ne -1 then N_GRP=rmfstr.N_GRP else ok='RMFSTR.N_GRP not found'
iF_CHAN=where(strpos(rtags,'F_CHAN') ge 0)
if iF_CHAN[0] ne -1 then F_CHAN=rmfstr.F_CHAN else ok='RMFSTR.F_CHAN not found'
iN_CHAN=where(strpos(rtags,'N_CHAN') ge 0)
if iN_CHAN[0] ne -1 then N_CHAN=rmfstr.N_CHAN else ok='RMFSTR.N_CHAN not found'
iMATRIX=where(strpos(rtags,'MATRIX') ge 0)
if iMATRIX[0] ne -1 then MATRIX=rmfstr.MATRIX else ok='RMFSTR.MATRIX not found'
iFIRSTCHAN=where(strpos(rtags,'FIRSTCHAN') ge 0)
if iFIRSTCHAN[0] ne -1 then FIRSTCHAN=rmfstr.FIRSTCHAN else ok='RMFSTR.FIRSTCHAN not found'
if ok ne 'ok' then begin
message,ok,/informational
if vv gt 0 then message,'Returning without further ado',/informational
return,-1L
endif
;
areff=total(matrix,1)
nea=n_elements(effar) & if nea eq nnrg then areff=areff*effar
;
wts=ehi-elo
nsp=n_elements(spec) & if nsp eq nnrg then wts=wts*spec

; make a full 2D image of the RMF
rmfimg=fltarr(nchan,nnrg) & szn=size(n_chan)
for inrg=0L,nnrg-1L do begin
ngrp=n_grp[inrg] & jbeg=0
for ig=0,ngrp-1 do begin
if szn[0] gt 1 then begin
ibeg=f_chan[ig,inrg] & iw=n_chan[ig,inrg]
endif else begin
ibeg=f_chan[inrg] & iw=n_chan[inrg]
endelse
if keyword_set(firstchan) then ibeg=ibeg-1 ;IDL index correction
if iw gt 0 then rmfimg[ibeg:ibeg+iw-1,inrg]=matrix[jbeg:jbeg+iw-1,inrg]
jbeg=jbeg+iw
endfor
endfor

; now compute the average weighted ARF in each detector channel
carf=fltarr(nchan)
for ic=0L,nchan-1L do carf[ic]=total(rmfimg[ic,*]*wts*areff)/total(rmfimg[ic,*]*wts)
chnrg=fltarr(nchan) & nrg=0.5*(elo+ehi)
for ic=0L,nchan-1L do chnrg[ic]=total(rmfimg[ic,*]*wts*nrg)/total(rmfimg[ic,*]*wts)

; show
if vv gt 10 then begin
plot,elo,effar,thick=3,line=1 & oplot,emn,carf,thick=3,line=2
if vv gt 1000 then stop
endif

return,carf
end