;+
;Suzaku_ao2.par
;	generate a simulated Suzaku XIS spectrum for a coronal source
;	using a specified multi-temperature DEM, normalized to a
;	Chandra/ACIS-I count rate, and also report the HXD rate. 
;
;description
;	http://hea-www.harvard.edu/PINTofALE/doc/EXAMPLE_Suzaku.html
;
;calling sequence
;	.run initale
;	@Suzaku_ao2.par
;
;(VK Dec2005)
;(bug corrections; VK 4Jan2006)
;-

;------------------------------------------------------------------------------
;	control variables
;------------------------------------------------------------------------------
;	Local Environment: set up the pathnames specific to installation here

pimmsdir	= '/soft/pimms/data'	;change to your local installation

!LDBDIR		= '$SPEX'			; Atomic Line Database
	        ; choose from the predefined '$CHIANTI', '$SPEX', '$APED',
		; or specify the full path name to the line database
!IONEQF		= 'ioneq/mazzotta_etal.ioneq'	; ion balance
!CEROOT		= 'cie'			; 'cie' or 'chianti'

;	Source: characterize model source by the following:

!NH		= 3e20			 	; H column density [cm^-2]
!EDENS		= 1.0e9   		 	; electron number density [cm^-3]
!ABUND		= getabund('asplund et al.') 	; element abundances
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]
T_components	= [6.8]		; log(T[K]) components in EM
EM_components	= [7.1e11]	; Emission Measure [cm^-3]

;	Observation parameters

EXPTIME		= 10.			   	; nominal exposure time [ks]
acisi_rate	= 1.0	; Chandra/ACIS-I rate [ct/s] - 0 or less to ignore
acisi_emin	= 0.3	; minimum ACIS-I energy to consider [keV]
acisi_emax	= 8.0	; maximum ACIS-I energy to consider [keV]
cornorm		= 0.0				; background correction norm

;	ARFs, RMFs, and background PHAs
;	NOTE: 'none' is always an option

;	XIS files can be downloaded from
;	http://heasarc.gsfc.nasa.gov/docs/astroe/prop_tools/xis_mat.html

xisFI_RMF	= 'ae_xi2_20051219.rmf'		; XIS FI RMF filename
xisFI_ARF	= 'ae_xi2_onaxis_20050916.arf'	; XIS FI ARF filename
xisFI_BKG	= 'ae_xi2_back.pha'		; XIS FI bkg PHA filename
;
xisBI_RMF	= 'ae_xi1_20051210d.rmf'	; XIS BI RMF filename
xisBI_ARF	= 'ae_xi1_onaxis_20050916.arf'	; XIS BI ARF filename
xisBI_BKG	= 'ae_xi1_back.pha'		; XIS BI bkg PHA filename

;	HXD files can be downloaded from
;	http://heasarc.gsfc.nasa.gov/docs/astroe/prop_tools/hxd_mat.html

hxdPIN_RSP	= 'ae_hxd_pinxinom_20051104.rsp'	;HXD PIN RSP filename
hxdPIN_BKG	= 'ae_hxd_pinbkg_20051105.pha'		;HXD PIN bkg PHA filename
;
hxdGSO_RSP	= 'ae_hxd_gso_20051019.rsp'		;HXD GSO RSP filename
hxdGSO_BKG	= 'ae_hxd_gsobkg_20051105.pha'		;HXD GSO bkg PHA filename

;	flow control parameters
;these variables will be filled in downstream, and will not be recomputed in
;subsequent runs of this script unless the appropriate lines are uncommented

;  ae2_lstr	= 0	;the line emissivity database structure
;  ae2_cstr	= 0	;the XIS continuum emissivity database structure
;  ae2_hcstr	= 0	;the HXD continuum emissivity database structure
;  rm_xisFI	= 0	;the XIS/FI response matrix
;  ea_xisFI	= 0	;the XIS/FI effective area
;  bg_xisFI	= 0	;the XIS/FI background spectrum
;  rm_xisBI	= 0	;the XIS/BI response matrix
;  ea_xisBI	= 0	;the XIS/BI effective area
;  bg_xisBI	= 0	;the XIS/BI background spectrum
;  rm_hxdPIN	= 0	;the HXD/PIN response matrix
;  bg_hxdPIN	= 0	;the HXD/PIN background spectrum
;  rm_hxdGSO	= 0	;the HXD/GSO response matrix
;  bg_hxdGSO	= 0	;the HXD/GSO background spectrum

;------------------------------------------------------------------------------
;	compute the 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

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

;------------------------------------------------------------------------------
;	read in the Suzaku calibration data
;------------------------------------------------------------------------------

;	XIS/FI
dwvl=0.000434 & xis_FIwgrid=10.^(findgen(4097)*dwvl)
if strlowcase(xisFI_RMF) ne 'none' and n_tags(rm_xisFI) eq 0 then $
  rm_xisFI=rd_ogip_rmf(xisFI_RMF,effar=qe_xisFI)
if n_tags(rm_xisFI) ne 0 then xis_FIwgrid=12.3985/[min(rm_xisFI.ELO),rm_xisFI.EHI]
if strlowcase(xisFI_RMF) eq 'none' and n_tags(rm_xisFI) eq 0 then qe_xisFI=0.*xis_FIwgrid[1:*]+1.
if strlowcase(xisFI_ARF) ne 'none' and not keyword_set(ea_xisFI) then ea_xisFI=rdarf(xisFI_ARF,ea_xisFIstr)
if strlowcase(xisFI_ARF) eq 'none' and not keyword_set(ea_xisFI) then ea_xisFI=0.*xis_FIwgrid[1:*]+1.
if strlowcase(xisFI_BKG) ne 'none' and n_tags(bg_xisFI) eq 0 then bg_xisFI=mrdfits(xisFI_BKG,1,hbg_xisFI)
if strlowcase(xisFI_BKG) eq 'none' and n_tags(bg_xisFI) eq 0 then $
  bg_xisFI=create_struct('CHANNEL',lindgen(n_elements(xis_FIwgrid)-1),'COUNTS',lonarr(0*xis_FIwgrid[1:*]))
bg_xisFI_exp=1. & if keyword_set(hbg_xisFI) then bg_xisFI_exp=float(sxpar(hbg_xisFI,'EXPOSURE'))>1.

;	XIS/BI
dwvl=0.000434 & xis_BIwgrid=10.^(findgen(4097)*dwvl)
if strlowcase(xisBI_RMF) ne 'none' and n_tags(rm_xisBI) eq 0 then $
  rm_xisBI=rd_ogip_rmf(xisBI_RMF,effar=qe_xisBI)
if n_tags(rm_xisBI) ne 0 then xis_BIwgrid=12.3985/[min(rm_xisBI.ELO),rm_xisBI.EHI]
if strlowcase(xisBI_RMF) eq 'none' and n_tags(rm_xisBI) eq 0 then qe_xisBI=0.*xis_BIwgrid[1:*]+1.
if strlowcase(xisBI_ARF) ne 'none' and not keyword_set(ea_xisBI) then ea_xisBI=rdarf(xisBI_ARF,ea_xisBIstr)
if strlowcase(xisBI_ARF) eq 'none' and not keyword_set(ea_xisBI) then ea_xisBI=0.*xis_BIwgrid[1:*]+1.
if strlowcase(xisBI_BKG) ne 'none' and n_tags(bg_xisBI) eq 0 then bg_xisBI=mrdfits(xisBI_BKG,1,hbg_xisBI)
if strlowcase(xisBI_BKG) eq 'none' and n_tags(bg_xisBI) eq 0 then $
  bg_xisBI=create_struct('CHANNEL',lindgen(n_elements(xis_BIwgrid)-1),'COUNTS',lonarr(0*xis_BIwgrid[1:*]))
bg_xisBI_exp=1. & if keyword_set(hbg_xisBI) then bg_xisBI_exp=float(sxpar(hbg_xisBI,'EXPOSURE'))>1.

;	HXD/PIN
dwvl=0.0094 & hxd_PINwgrid=10.^(findgen(257)*dwvl-0.89)
if strlowcase(hxdPIN_RSP) ne 'none' and n_tags(rm_hxdPIN) eq 0 then $
  rm_hxdPIN=rd_ogip_rmf(hxdPIN_RSP,effar=ea_hxdPIN)
if n_tags(rm_hxdPIN) ne 0 then hxd_PINwgrid=12.3985/[min(rm_hxdPIN.ELO),rm_hxdPIN.EHI]
if n_tags(rm_hxdPIN) ne 0 then ea_hxdPIN=total(rm_hxdPIN.MATRIX,1) else ea_hxdPIN=0.*hxd_PINwgrid[1:*]+1.
if strlowcase(hxdPIN_BKG) ne 'none' and n_tags(bg_hxdPIN) eq 0 then bg_hxdPIN=mrdfits(hxdPIN_BKG,1,hbg_hxdPIN)
if strlowcase(hxdPIN_BKG) eq 'none' and n_tags(bg_hxdPIN) eq 0 then $
  bg_hxdPIN=create_struct('CHANNEL',lindgen(n_elements(hxd_PINwgrid)-1),'COUNTS',lonarr(0*hxd_PINwgrid[1:*]))
bg_hxdPIN_exp=1. & if keyword_set(hbg_hxdPIN) then bg_hxdPIN_exp=float(sxpar(hbg_hxdPIN,'EXPOSURE'))>1.

;	HXD/GSO
dwvl=0.0059 & hxd_GSOwgrid=10.^(findgen(513)*dwvl-1.92)
if strlowcase(hxdGSO_RSP) ne 'none' and n_tags(rm_hxdGSO) eq 0 then $
  rm_hxdGSO=rd_ogip_rmf(hxdGSO_RSP,effar=ea_hxdGSO)
if n_tags(rm_hxdGSO) ne 0 then hxd_GSOwgrid=12.3985/[min(rm_hxdGSO.ELO),rm_hxdGSO.EHI]
if n_tags(rm_hxdGSO) ne 0 then ea_hxdGSO=total(rm_hxdGSO.MATRIX,1) else ea_hxdGSO=0.*hxd_GSOwgrid[1:*]+1.
if strlowcase(hxdGSO_BKG) ne 'none' and n_tags(bg_hxdGSO) eq 0 then bg_hxdGSO=mrdfits(hxdGSO_BKG,1,hbg_hxdGSO)
if strlowcase(hxdGSO_BKG) eq 'none' and n_tags(bg_hxdGSO) eq 0 then $
  bg_hxdGSO=create_struct('CHANNEL',lindgen(n_elements(hxd_GSOwgrid)-1),'COUNTS',lonarr(0*hxd_GSOwgrid[1:*]))
bg_hxdGSO_exp=1. & if keyword_set(hbg_hxdGSO) then bg_hxdGSO_exp=float(sxpar(hbg_hxdGSO,'EXPOSURE'))>1.

;------------------------------------------------------------------------------
;	read in line emissivity database
;------------------------------------------------------------------------------

wmin=min(xis_BIwgrid) < min(xis_FIwgrid)
wmax=max(xis_BIwgrid) < max(xis_FIwgrid)
incieq=1
if strpos(strlowcase(!LDBDIR),'apec',0) ge 0 then incieq=0 else $
 if strpos(strlowcase(!LDBDIR),'aped',0) ge 0 then incieq=0 else $
  if strpos(strlowcase(!LDBDIR),'gu',0) ge 0 then incieq=0
if n_tags(ae2_lstr) eq 0 then $
  ae2_lstr=rd_list('all|'+strtrim(wmin,2)+'-'+strtrim(wmax,2)+'|'+!LDBDIR,sep='|',$
  incieq=incieq,n_e=!EDENS,eqfile=!IONEQF,verbose=!VERBOSE)

;	If !LDBDIR is set to APED, Anders & Grevesse abundances
;	are already included in the emissivities.  In which case,
;	redefine the abundance array to be relative to APED.

abnorm=0*!ABUND+1.
if strpos(strlowcase(!LDBDIR),'apec',0) ge 0 or $
  strpos(strlowcase(!LDBDIR),'aped',0) ge 0 then $
  abnorm=getabund('anders & grevesse')

;------------------------------------------------------------------------------
;	read in continuum emissivity database
;------------------------------------------------------------------------------

if n_tags(ae2_cstr) eq 0 then $
  fc=rd_cont(!CEROOT,n_e=!EDENS,wrange=[wmin,wmax],dbdir=!CDBDIR,$
  abund=!ABUND,fcstr=ae2_cstr,twoph=twoph)

if n_tags(ae2_hcstr) eq 0 then $
  fc=rd_cont(!CEROOT,n_e=!EDENS,wrange=12.3985/[60.,1.],dbdir=!CDBDIR,$
  abund=!ABUND,fcstr=ae2_hcstr,twoph=twoph)

;------------------------------------------------------------------------------
;	do not edit anything below unless you know what you are doing
;------------------------------------------------------------------------------

;------------------------------------------------------------------------------
;	compute nominal fluxes and ISM opacities
;------------------------------------------------------------------------------

;	line intensities
linint=	lineflx(ae2_lstr.LINE_INT,ae2_lstr.LOGT,abs(ae2_lstr.WVL),$
	ae2_lstr.Z,DEM=!DEM,abund=!ABUND/abnorm)		;[ph/s/cm^2]

;	continuum intensities
conint=	lineflx(ae2_cstr.CONT_INT,ae2_cstr.LOGT,abs(ae2_cstr.midWVL),$
	DEM=!DEM)						;[ph/s/cm^2/Ang]
cwb=	mid2bound(ae2_cstr.midWVL)
cdw=	cwb[1:*]-cwb
conint=	conint*cdw						;[ph/s/cm^2]

;	continuum intensities for HXD
hconint=lineflx(ae2_hcstr.CONT_INT,ae2_hcstr.LOGT,abs(ae2_hcstr.midWVL),$
	DEM=!DEM)						;[ph/s/cm^2/Ang]
cwb=	mid2bound(ae2_hcstr.midWVL)
cdw=	cwb[1:*]-cwb
hconint=hconint*cdw						;[ph/s/cm^2]

;	ISM absorption
ltau=	ismtau(abs(ae2_lstr.WVL),NH=!NH,fH2=!fH2,He1=!He1,HeII=!HeII,$
	Fano=Fano,wam=wam,/bam,abund=!ABUND,verbose=!VERBOSE)
ctau=	ismtau(abs(ae2_cstr.midWVL),NH=!NH,fH2=!fH2,He1=!He1,HeII=!HeII,$
	Fano=Fano,wam=wam,/bam,abund=!ABUND,verbose=!VERBOSE)
hctau=	ismtau(abs(ae2_hcstr.midWVL),NH=!NH,fH2=!fH2,He1=!He1,HeII=!HeII,$
	Fano=Fano,wam=wam,/bam,abund=!ABUND,verbose=!VERBOSE)
ltrans=	exp(-ltau)
ctrans=	exp(-ctau)
hctrans=exp(-hctau)

;	line fluxes
linflx=	linint * ltrans		; [ph/s/cm^2]
conflx=	conint * ctrans		; [ph/s/cm^2]
hconflx=hconint * hctrans	; [ph/s/cm^2]

;------------------------------------------------------------------------------
;	renormalize DEM by comparing to predicted ACIS rate
;------------------------------------------------------------------------------

;	We will assume that a Chandra/ACIS-I count rate is available,
;	and require that the simulation match this rate.  Data from
;	other missions such as ROSAT, ASCA, BeppoSAX, etc. can be dealt
;	with in the same manner.

;	First find and read in the Chandra/ACIS-I effective area
;	you will need to know where your local PIMMS installation
;	is to do this.

rd_pimms_file,get_pimms_file('chandra','acis-i',pdir=pimmsdir),$
	acisi_effar,acisi_wvlar,/wave

;	fold in effective area
acisi_wmin=12.3985/acisi_emax & acisi_wmax=12.3985/acisi_emin
os=sort(acisi_wvlar) & acisi_wvlar=acisi_wvlar[os] & acisi_effar=acisi_effar[os]
;
lwvl=abs(ae2_lstr.WVL) & ow_lacisi=where(lwvl ge acisi_wmin and lwvl le acisi_wmax)
effar=(interpol(acisi_effar,acisi_wvlar,lwvl[ow_lacisi])>0) < max(acisi_effar)
flxlin_acisi=	linflx[ow_lacisi]*effar		;[ct/s]
;
cwvl=abs(ae2_cstr.midWVL) & ow_cacisi=where(cwvl ge acisi_wmin and cwvl le acisi_wmax)
effar=(interpol(acisi_effar,acisi_wvlar,cwvl)>0) < max(acisi_effar)
flxcon_acisi=	conflx[ow_cacisi]*effar		;[ct/s]

;	total predicted rate in Chandra/ACIS-I
pred_rate=total(flxlin_acisi)+total(flxcon_acisi)

print,''
if acisi_rate gt 0 then $
  print,'Rescaling input DEM by a factor '+strtrim(acisi_rate/pred_rate,2)
  print,'to match required ACIS-I rate of '+strtrim(acisi_rate,2)+' ct/s'
print,''
rescale_factor=1.0
if acisi_rate gt 0 then rescale_factor = acisi_rate/pred_rate
!DEM = !DEM * rescale_factor
linflx	= linflx * rescale_factor
conflx	= conflx * rescale_factor
hconflx	= hconflx * rescale_factor

;------------------------------------------------------------------------------
;	compute HXD rate
;------------------------------------------------------------------------------

;	HXD/PIN
wv_hxdPIN=0.5*(hxd_PINwgrid[1:*]+hxd_PINwgrid)
;
lwvl=abs(ae2_lstr.WVL)
effar=(interpol(ea_hxdPIN,wv_hxdPIN,lwvl)>0) < max(ea_hxdPIN)
flxlin=	linflx*effar		;[ct/s]
;
cwvl=abs(ae2_hcstr.midWVL)
effar=(interpol(ea_hxdPIN,wv_hxdPIN,cwvl)>0) < max(ea_hxdPIN)
flxcon=	hconflx*effar		;[ct/s]
;
flx_hxdPIN = total(flxlin)+total(flxcon)

;	HXD/GSO
wv_hxdGSO=0.5*(hxd_GSOwgrid[1:*]+hxd_GSOwgrid)
;
lwvl=abs(ae2_lstr.WVL)
effar=(interpol(ea_hxdGSO,wv_hxdGSO,lwvl)>0) < max(ea_hxdGSO)
flxlin=	linflx*effar		;[ct/s]
;
cwvl=abs(ae2_hcstr.midWVL)
effar=(interpol(ea_hxdGSO,wv_hxdGSO,cwvl)>0) < max(ea_hxdGSO)
flxcon=	hconflx*effar		;[ct/s]
;
flx_hxdGSO = total(flxlin)+total(flxcon)

;------------------------------------------------------------------------------
;	compute XIS rate
;------------------------------------------------------------------------------

;	XIS/FI
wv_XISFI=0.5*(XIS_FIwgrid[1:*]+XIS_FIwgrid)
;
lwvl=abs(ae2_lstr.WVL)
effar=(interpol(qe_XISFI*ea_XISFI,wv_XISFI,lwvl)>0) < max(qe_XISFI*ea_XISFI)
flxlin=	linflx*effar		;[ct/s]
;
cwvl=abs(ae2_cstr.midWVL)
effar=(interpol(qe_XISFI*ea_XISFI,wv_XISFI,cwvl)>0) < max(qe_XISFI*ea_XISFI)
flxcon=	conflx*effar		;[ct/s]
;
flx_XISFI = total(flxlin)+total(flxcon)

;	XIS/BI
wv_XISBI=0.5*(XIS_BIwgrid[1:*]+XIS_BIwgrid)
;
lwvl=abs(ae2_lstr.WVL)
effar=(interpol(qe_XISBI*ea_XISBI,wv_XISBI,lwvl)>0) < max(qe_XISBI*ea_XISBI)
flxlin=	linflx*effar		;[ct/s]
;
cwvl=abs(ae2_cstr.midWVL)
effar=(interpol(qe_XISBI*ea_XISBI,wv_XISBI,cwvl)>0) < max(qe_XISBI*ea_XISBI)
flxcon=	conflx*effar		;[ct/s]
;
flx_XISBI = total(flxlin)+total(flxcon)

;------------------------------------------------------------------------------
;	compute XIS spectra
;------------------------------------------------------------------------------

;       To construct the Suzaku/XIS spectra, we use the line and continuum
;	fluxes calculated above (LINFLX and CONFLX), rebin them into a
;	consistent energy grid, and convolve with the response matrix

;	XIS/FI

;	bin the spectra

egrid=12.3985/xis_FIwgrid & os=sort(egrid) & egrid=egrid[os] & tmp=egrid
lnrg=12.3985/abs(ae2_lstr.WVL) & os=sort(lnrg)
xisFI_lspc=	hastogram(lnrg[os],egrid,wts=linflx[os])	;[ph/s/cm^2/bin]
;
egrid=12.3985/xis_FIwgrid & os=sort(egrid) & egrid=egrid[os]
cnrg=12.3985/ae2_cstr.WVL & os=sort(cnrg)
xisFI_cspc=	rebinw(conflx[os],cnrg[os],egrid,/perbin)	;[ph/s/cm^2/bin]
;
xisFI_spc=	xisFI_lspc + xisFI_cspc

;	convolve with RMF
tmp=12.3985/xis_FIwgrid & ee_xisFI=0.5*(tmp[1:*]+tmp)
if strlowcase(xisFI_RMF) ne 'none' then $
  conv_rmf,egrid,xisFI_spc,xisFI_channel,xisFI_predct,rm_xisFI,$
  effar=ea_xisFI,nrgar=0.5*(ea_xisFIstr.ELO+ea_xisFIstr.EHI),$
  verbose=!VERBOSE

;	units are now in ct/s/bin, so multiply by EXPTIME
xisFI_predct = xisFI_predct * EXPTIME * 1e3	;[ct/bin]

;	XIS/BI

;	bin the spectra
egrid=12.3985/xis_BIwgrid & os=sort(egrid) & egrid=egrid[os] & tmp=egrid
lnrg=12.3985/abs(ae2_lstr.WVL) & os=sort(lnrg)
xisBI_lspc=	hastogram(lnrg[os],egrid,wts=linflx[os])	;[ph/s/cm^2/bin]
;
egrid=12.3985/xis_BIwgrid & os=sort(egrid) & egrid=egrid[os]
cnrg=12.3985/ae2_cstr.WVL & os=sort(cnrg)
xisBI_cspc=	rebinw(conflx[os],cnrg[os],egrid,/perbin)	;[ph/s/cm^2/bin]
;
xisBI_spc=	xisBI_lspc + xisBI_cspc

;	convolve with RMF
tmp=12.3985/xis_BIwgrid & ee_xisBI=0.5*(tmp[1:*]+tmp)
if strlowcase(xisBI_RMF) ne 'none' then $
  conv_rmf,egrid,xisBI_spc,xisBI_channel,xisBI_predct,rm_xisBI,$
  effar=ea_xisBI,nrgar=0.5*(ea_xisBIstr.ELO+ea_xisBIstr.EHI),$
  verbose=!VERBOSE

;	units are now in ct/s/bin, so multiply by EXPTIME
xisBI_predct = xisBI_predct * EXPTIME * 1e3	;[ct/bin]

;------------------------------------------------------------------------------
;	obtain Poisson deviates of source and background spectra
;------------------------------------------------------------------------------

print,'simualting XIS/FI source spectrum'
nFI=n_elements(xisFI_predct)
xisFI_counts=lonarr(nFI) & xisFI_countsbg=lonarr(nFI) & xisFI_countstot=lonarr(nFI)
tmp=xisFI_predct & for i=0L,nFI-1L do if tmp[i] gt 1e-20 then $
	xisFI_counts[i]=randomu(seed,poisson=tmp[i])
;
print,'simualting XIS/FI background spectrum'
tmp=bg_xisFI.counts*(EXPTIME*1e3)*(1.+cornorm)/BG_XISFI_EXP
for i=0L,nFI-1L do if tmp[i] gt 1e-20 then $
	xisFI_countsbg[i]=randomu(seed,poisson=tmp[i])
;
xisFI_countstot = xisFI_counts + xisFI_countsbg

print,'simualting XIS/BI source spectrum'
nBI=n_elements(xisBI_predct)
xisBI_counts=lonarr(nBI) & xisBI_countsbg=lonarr(nBI) & xisBI_countstot=lonarr(nBI)
tmp=xisBI_predct & for i=0L,nBI-1L do if tmp[i] gt 1e-20 then $
	xisBI_counts[i]=randomu(seed,poisson=tmp[i])
;
print,'simualting XIS/FI background spectrum'
tmp=bg_xisBI.counts*(EXPTIME*1e3)*(1.+cornorm)/BG_XISBI_EXP
for i=0L,nBI-1L do if tmp[i] gt 1e-20 then $
	xisBI_countsbg[i]=randomu(seed,poisson=tmp[i])
;
xisBI_countstot = xisBI_counts + xisBI_countsbg

;------------------------------------------------------------------------------
;	make plots
;------------------------------------------------------------------------------

peasecolr & !p.multi=[0,1,2]

plot,xisFI_channel,xisFI_countstot,psym=10,thick=2,$
	xtitle='E [keV]',ytitle='counts',$
	title='XIS/FI '+strtrim(string(EXPTIME,'(f5.2)'),2)+' ks',$
	xrange=[0.2,8],/xs,/xlog,yrange=max(xisFI_countstot)*[1e-4,1],/ylog
oplot,xisFI_channel,xisFI_counts,psym=10,thick=2,col=1
oplot,xisFI_channel,xisFI_countsbg,psym=10,thick=2,col=2
oplot,xisFI_channel,xisFI_predct,col=3

plot,xisBI_channel,xisBI_countstot,psym=10,thick=2,$
	xtitle='E [keV]',ytitle='counts',$
	title='XIS/BI '+strtrim(string(EXPTIME,'(f5.2)'),2)+' ks',$
	xrange=[0.2,8],/xs,/xlog,yrange=max(xisBI_countstot)*[1e-4,1],/ylog
oplot,xisBI_channel,xisBI_counts,psym=10,thick=2,col=1
oplot,xisBI_channel,xisBI_countsbg,psym=10,thick=2,col=2
oplot,xisBI_channel,xisBI_predct,col=3
stample,stacol=4

!p.multi=0

;------------------------------------------------------------------------------
;	report
;------------------------------------------------------------------------------

print,''
print,'HXD/PIN bremss count rate = ',flx_hxdPIN,'	total counts = ',flx_hxdPIN*EXPTIME*1e3
print,'HXD/PIN background counts = ',(1.+cornorm)*total(bg_HXDPIN.COUNTS)*EXPTIME*1e3/bg_hxdPIN_exp
print,''
print,'HXD/GSO bremss count rate = ',flx_hxdGSO,'	total counts = ',flx_hxdGSO*EXPTIME*1e3
print,'HXD/GSO background counts = ',(1.+cornorm)*total(bg_HXDGSO.COUNTS)*EXPTIME*1e3/bg_hxdGSO_exp
print,''

print,''
print,'XIS/FI count rate = ',flx_XISFI,'	total counts = ',flx_XISFI*EXPTIME*1e3
print,'XIS/FI background counts = ',(1.+cornorm)*total(bg_XISFI.COUNTS)*EXPTIME*1e3/bg_XISFI_exp
print,''
print,'XIS/BI count rate = ',flx_XISBI,'	total counts = ',flx_XISBI*EXPTIME*1e3
print,'XIS/BI background counts = ',(1.+cornorm)*total(bg_XISBI.COUNTS)*EXPTIME*1e3/bg_XISBI_exp
print,''

print,''
print,'------------------------------------------------------------------------------'
print,'	Outputs are in the following variables:'
print,'------------------------------------------------------------------------------'
print,''
print,'------------------------------------------------------------------------------'
print,'XIS/FI --'
print,'model spectrum, pre-RMF :	xisFI_spc [ph/s/cm^2/bin] (egrid,xis_FIwgrid)'
print,'predicted counts intensity :	xisFI_predct [ct/channel] (xisFI_channel)'
print,'simulated source counts :	xisFI_counts [ct/channel] (xisFI_channel)'
print,'simulated background counts :	xisFI_countsbg [ct/channel] (xisFI_channel)'
print,'simulated total counts  :	xisFI_countstot [ct/channel] (xisFI_channel)'
print,''
print,'------------------------------------------------------------------------------'
print,'XIS/BI --'
print,'model spectrum, pre-RMF :	xisBI_spc [ph/s/cm^2/bin] (egrid,xis_BIwgrid)'
print,'predicted counts intensity :	xisBI_predct [ct/channel] (xisBI_channel)'
print,'simulated source counts :	xisBI_counts [ct/channel] (xisBI_channel)'
print,'simulated background counts :	xisBI_countsbg [ct/channel] (xisBI_channel)'
print,'simulated total counts  :	xisBI_countstot [ct/channel] (xisBI_channel)'
print,''
print,'------------------------------------------------------------------------------'
print,'HXD/PIN --'
print,'bremss count rate :	flx_hxdPIN [ph/s]'
print,''
print,'------------------------------------------------------------------------------'
print,'HXD/GSO --'
print,'bremss count rate :	flx_hxdGSO [ph/s]'
print,'------------------------------------------------------------------------------'
print,''