function get_be_add, x, SET=SET

common get_be_add, xsave

if keyword_set(set) and exist(x) then xsave = x 

xsave = fcheck( x, fcheck(xsave, 0.35))

return, xsave

end

;+
;  Name:
;       BATSE_SPEC_DRM
;	
; PURPOSE: build batse spec detector response matrix
;	including proper treatment of detector energyloss/light-output
;	relationship
;	
; CATEGORY: BATSE, SPEX
;
; CALLING SEQUENCE:
;	
; INPUTS:
;	edges_out - drm output channel edges (2xn) in keV
;	theta - zenith angle of source, theta is 0 degrees for a source on-axis
;	det_id- detector id, 0-7 for BATSE spec
;	high- if set include a higher density above an MeV on the input side
; OPTIONAL INPUTS:
;	
; OUTPUTS:
;	edges_in - photon intervals used for drm input
;	drm - detector response matrix 
;	pairfit - coefficients of fits to energies above pair-threshold energy
;	sigma - resolution in standard deviations (sigma) from fits to response
;		database
;	e_in - energies corresponding to sigma in keV, original database photon energies
;	phtn_e- use this set of incident photon energies, 2xn, in keV
;
; OPTIONAL OUTPUTS:
;
; PROCEDURE:
;
; CALLS:
;
; COMMON BLOCKS:
;
; RESTRICTIONS:
;
; MODIFICATION HISTORY:
;	ras, 10-mar-1995
;	ras, 18-mar-1995, replaced spline interpolation with linear interpolation
;	for photoefficiency
;	added normalization needed when changing from energy-loss to light-output gaussians
;	increased output bin sampling internally, then summed over to final resolution, ras, 27jul95
;	fixed final rebinning to ignore output channels where the width is 0.0, ras, 25-sep-1995
;	ras, 1-nov-95
;	because pair fits can be bad at high energy, had to check pairfit quality range,
;
;	version 4, ras 17 oct 1996, added phtn_e
;	Version 5, richard.schwartz@gsfc.nasa.gov, 18-feb-1998, protect against NAN universally using
;	idl's finite function, not the woeful ssw routine.
;-
pro batse_spec_drm, edges_out=edges_out, edges_in=edges_in2, theta=theta, $
	high=high, det_id=det_id, drm=drm, resid=drm_resid, pairfit=pairfit, $
	scale_sigma=scale_sigma, sigma=sigma, e_in=e_in, phtn_e=phtn_e

;get the database drm decomposed into photopeak, escape peak, and residual flux

decompose_batse_drm, theta=theta, det_id=det_id, det_type='spec', e_in=e_in,  $
	photo_eff=photo_eff, esc_eff=esc_eff, drm_resid=drm_resid, out=out, $
	sigma = sigma, resp=resp, gm_use=gm_use, verbose=0

nin = n_elements(e_in)
sigma = sigma*e_in*fcheck(scale_sigma,1.)
;do explicit fits starting above 1600 keV for pair-production escape and double escape

i1 = (where(e_in ge 1600., ni1))(0)
epair = e_in(i1) ;first pair energy

pairfit = fltarr(12,ni1 > 1)
if ni1 ge 0 then begin
  for i=i1,nin - 1 do begin 
   w=where(e_in(i)*out gt .6*(e_in(i)*gm_use(1,i)-1022.) , nw)
   x = e_in(i)*out(w)
   y = resp(w,i)
   afit = [0,0,0, gm_use(0:2,i)*[1,e_in(i),e_in(i)*1.414]]
   afit = [afit, afit(3:5), afit(3:5)]
   afit([7,10])=afit([7,10])-[511.,1022.] 
   afit([6,9])=afit([6,9])*[.2,.1]   
   afit(0)=avg(y-mgauss( x, afit))
   afit([5,8,11]) =interpol(sigma,e_in,afit([4,7,10]))*sqrt(2.)
   ;yfit = curvfit(x,y,y*0.0+1.,niter=5, /quiet, afit,funct='mgauss')
   yfit = funct_fit(x,y, afit,funct='mgauss')
   pairfit(0,i-i1)=afit
  endfor
;Evaluate the quality of the fits and replace bad fits with interpolated values
   wgood = where( (abs(pairfit(4,*)-pairfit(7,*)-511.) le 200.) and $
		 (abs(pairfit(4,*)-pairfit(10,*)-1022.) le 200.), ngood )
   if ngood ge 1 then begin
	pairfit = pairfit(*,wgood)
	epair = minmax(pairfit(4,*)) ;starting energy for pairfits
   endif
	
endif

;are the output channels sufficiently resolved to avoid edge effects in sampling?
;
w50=where(edges_out(0,*) le 50., nw50)
usenew = 0
if nw50 ge 1 then begin
	if (where( edges_out(1,w50)-edges_out(0,w50) gt 2.))(0) ne -1 then begin
		edge_products, edges_out, edges_1=edges_out1
		nedges = [edges_out1, 5.+findgen(50),findgen(10)*5.+55.,findgen(90)*10.+110.]
		ord = uniq( nedges, sort(nedges))
		iold = long( interpol(indgen(n_elements(ord)), nedges(ord), edges_out1))
		edge_products, nedges(ord), edges_2=edges
		usenew = 1
	endif
endif
if not usenew then edges = edges_out

edge_products, edges, edges_2=edges, mean=eout_mean, edges_1=evec, wid=wedges

nout = n_elements(eout_mean)

if n_elements(phtn_e) eq 0 then begin
  e_in2 = [e_in,findgen(38)+3.,33.+findgen(40)*0.25,6.25+findgen(3)*0.25]
  if fcheck(high,0) then e_in2 = [e_in, findgen(50)*4+2000.,510.,512.,1000.+findgen(60)*50.]
  e_in2 = e_in2(uniq(e_in2, sort(e_in2)))
endif else e_in2 = phtn_e
edge_products, e_in2, mean=ein2_mean, edges_2=edges_in2
;now at least one bin every keV to 33 keV
;want the input energies interpolated over the edges_in2

nin2  = n_elements(ein2_mean)

drm = fltarr(nout, nin2)

drm_resid = drm_resid / 127.
photo_eff = photo_eff / 127.
esc_eff = esc_eff / 127.
;above an MeV add the photo_eff back into the drm_resid
w1000 = where(e_in ge 1000.0, nw1000)
;if nw1000 ge 1 then for i=w1000(0),w1000(nw1000-1) do drm_resid(0,i) = drm_resid(*,i) + $
;	gaussian( out*e_in(i),[photo_eff(i)/(sigma(i)*sqrt(2.*!pi)),e_in(i),sigma(i)] )
photo_eff2 = photo_eff
;photo_eff2(w1000) = 0.0 

w10 = where(edges_in2(0,*) gt 16., nw10)

;The indexes into out for eout_mean for a peak energy of 100.0 keV
ifid = interpol( findgen(n_elements(out)),out, eout_mean/100.) >0.0

for i=w10(0),nin2-1 do begin
  if edges_in2(0,i) le 1000 then nsub=2 else nsub = 2 
  esub=edges_in2(0,i)+(findgen(nsub)+1./(nsub*2))*(edges_in2(1,i)-edges_in2(0,i))/nsub   
  wm = where(eout_mean le edges_in2(1,i)*1.2, nwm)
  isub = interpol(findgen(nin), e_in, esub)
  if nwm ge 1 then for j=0, nsub-1 do begin
	drm(wm,i) = drm(wm,i) + interpolate(drm_resid,ifid(wm)*100/esub(j),fltarr(nwm)+ isub(j))/nsub
  endfor


endfor

g_c, evec, el
el = el*evec
edge_products, el, edges_2=elight, mean=eml, width=wl


for i=w10(0),nin2-1 do begin   
  wgt = where(eout_mean gt ein2_mean(i)*1.1, ngt)   
  if ngt ge 1 then drm(wgt,i) = 0   
endfor

new =intarr(nout)
for i=0,nout-1 do new(i) = (where(eml(i) le edges(0,*)))(0)

wbad = where(new eq -1, nbad)
if nbad ge 1 then new(wbad) = nout-1
map = fltarr(nout,nout)

newsigma = interpol(sigma,e_in, eout_mean)
for i=0,nout-1 do map(*,i) = f_div(wedges(i),wedges) * $
    (gaussint( (edges(1,*)-eout_mean(new(i)))/newsigma(new(i))) - $
     gaussint( (edges(0,*)-eout_mean(new(i)))/newsigma(new(i))))

drm = temporary(map#drm)
drm_resid = drm

wesc = where(esc_eff gt 0, nesc)

esc_energy = e_in - 29.2 > 4.

nsub = 10
;need to add the extra attenuation from the millipore paper at all energies.
photo_eff2 = photo_eff2 * low_e_spec_eff( e_in < 999., theta, be_add=get_be_add()) / $
	low_e_spec_eff(e_in < 999., theta, be_add=0.0)
;stop
for i=0,nin2-1 do begin   
   esub=edges_in2(0,i)+(findgen(nsub)+1./(nsub*2))*(edges_in2(1,i)-edges_in2(0,i))/nsub   
   sigsub = interpol(sigma,e_in,esub) 
;   if esub(0) le 12. then ph_sub = low_e_spec_eff(esub, theta, be_add=0.35)/127. else $
   ph_sub = interpol(photo_eff2,e_in,esub)
   ph_sub = ph_sub/nsub/sigsub/sqrt(2*!pi)   
   g_c, esub, esubl_old, light2= esubl
   esubl = esubl * esub   
   ;correct the normalization for integrating over light units!
   ;the number of interactions must be preserved!
   ;should revisit and check for low energy edge effects
   ph_normal = f_div( total( f_gauss_intg(edges, [ph_sub(nsub/2),esub(nsub/2),sigsub(nsub/2)])) ,$
	 total( f_gauss_intg(elight, [ph_sub(nsub/2),esubl(nsub/2),sigsub(nsub/2)])))
   ph_sub =ph_sub * ph_normal

  if edges_in2(0,i) ge epair(0) and edges_in2(1,i) le epair(1) then begin
	;clear out drm to be used by pairfit
	    wm = where( eout_mean ge .6*(edges_in2(0,i)-1022.),nwm)
	    if nwm ge 1 then drm(wm,i) = 0.0
;interpolate energies in esub by the energies of the full-energy peak!
	    ipair = interpol(findgen(n_elements(pairfit(4,*))), (pairfit(4,*))(*), esub)
	    xpair = ipair mod 1.
	    ipair = fix(ipair)
   endif

    if edges_in2(0,i) lt epair(0) or edges_in2(1,i) gt epair(1) then $

     for j=0,nsub-1 do drm(0,i) = drm(*,i)+ f_gauss_intg(elight, [ph_sub(j),esubl(j),sigsub(j)])   $
   else begin  ;use pair fits!

     for j=0,nsub-1 do begin 
	afit = pairfit(*,ipair(j))*(1.-xpair(j)) + $
	       xpair(j) * pairfit(*,(ipair(j)+1)<(n_elements(pairfit(4,*))-1) )
        afit([7,10]) = afit(4) -[511.,1022.]
	if nwm ge 1 then drm(wm,i)= drm(wm,i) + mgauss( eout_mean(wm), afit )/(nsub*127.)
     endfor
     wmm = where( eout_mean ge (afit(4)+3*(afit(5)/1.414)), nwmm)
     if nwmm ge 1 then drm(wmm,i) = 0.0 ;make sure it goes to zero above peak!
   endelse
	

endfor


nsub = 10
for i=0,nin2-1 do begin   
   if edges_in2(0,i) ge e_in(wesc(0)) and edges_in2(0,i) le e_in(wesc(nesc-1)) then begin   
	esub=edges_in2(0,i)+(findgen(nsub)+1./(nsub*2))*(edges_in2(1,i)-edges_in2(0,i))/nsub   
	esub= esub - 29.2   
   	sigsub = interpol(sigma,e_in,esub) & $
   	ph_sub = interpol(esc_eff,esc_energy,esub)/nsub/sigsub/sqrt(2*!pi)   
   	g_c, esub, esubl   
   	esubl = esubl * esub   
   	for j=0,nsub-1 do $
      	drm(0,i) = drm(*,i)+ f_gauss_intg(elight, [ph_sub(j),esubl(j),sigsub(j)])   
   endif   
endfor

;HAVE WE ADDED EDGES TO THE OUTPUT CHANNELS, IF SO WE HAVE TO INTEGRATE BACK TO THE ORIGINAL
;EDGES.
if usenew then begin
	drm1 = drm
	drm  = fltarr( n_elements(edges_out(0,*)),n_elements(drm1(0,*)) )
	edge_products, edges_out, wid=wedges_out
	awedges = rebin( reform(wedges, n_elements(wedges),1), n_elements(wedges), n_elements(drm1(0,*)))
	for i=0,n_elements(iold)-2 do if wedges_out(i) gt 0 then drm(i,*) = $
		total(drm1(iold(i):iold(i+1)-1,*)* awedges(iold(i):iold(i+1)-1,*),1) / wedges_out(i)
	awedges = 0
	drm1 = 0
endif


;protect against NAN
wnan = where(finite(drm) eq 0, num_nan)
if num_nan ge 1 then drm( wnan ) = 0.0
printx,string(/print,'minmax(drm)',minmax(drm))



end