pro funct2, x, a, f, pder ;FUNCT procedure to run with CURVEFIT common funct, function_name f = call_function( function_name, x, a(0:*) ) ;return the values if n_params() eq 4 then begin; CALCULATE PARTIAL DERIVATIVES npar = n_elements(a) pder = dblarr(n_elements(x), npar) for i=0,npar-1 do $ pder(*,i)=f_pder(funct=function_name, par=a(0:*), data=x, n=i) endif return end ;+ ; Name: ; DECOMPOSE_BATSE_DRM ; ; PURPOSE: ; Decompose a BATSE SPEC response matrix into photopeak line strengths ; and K escape peak line strengths as well as the residual matrix. ; For BATSE LAD, there is no decomposition, just the basic responses ; ; CALLING SEQUENCE: ; decompose_batse_drm, theta=theta, det_id=det_id, $ ; photo_eff=photo_eff, esc_eff=esc_eff, drm_resid=drm_resid, e_in, $ ; out_frac=out_frac, sigma = sigma ; ; INPUTS: ; theta - source zenith angle ; det_id - use input response file name, bsas/bdas database files, for this id ; files installed on SSWDB_BATSE ; det_type - 'spec' or 'lad' ; OPTIONAL INPUTS: ; only_resp - if set, then only generate the response for the angle ; don't extract any photo_efficiencies ; if only_resp set to -1 then even lad will have photoefficiencies extracted! ; OUTPUTS: ; e_in - photon energies, keV ; photo_eff- photoefficiency flux in counts/det/ (photon/cm2) ; esc_eff - escape efficiency ; drm_resid - residual drm after photo and esc peaks subtracted ; drm_resid is in units of counts/keV / (photon/cm2) ; out - energy grid of output energy loss in units of fraction of input photon energy ; grid points are centered in each output interval ; sigma - resolution sigma in keV, fwhm = 2.36 * sigma ; OPTIONAL OUTPUTS: ; gm_use - parameters of gaussians used to fit peaks and escape peaks ; resp - basic batse response matrix built from stored model coefficients, ; polynomial in cosine ; ; PROCEDURE: ; ; CALLS: ; ; COMMON BLOCKS: ; ; RESTRICTIONS: ; ; MODIFICATION HISTORY: ; ras, 5-apr-95, changed from jun93 to mar94 data files ; ras, 10-apr-95, changed back to mar94 to jun93 data files ; ras, 21-apr-95, changed to support lad drm too ; ras, 21-jul-95, start fitting at response row 3 skipping 0,1,2 ; Version 5, ras, 15-jun-1997, Allow only_resp to force photo_efficiencies! ; RAS, 7-Jul-1997, changed PERM_DATA to SSWDB_BATSE ;- ; pro decompose_batse_drm, theta=theta, det_id=det_id, det_type=det_type, 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=verbose, only_resp=only_resp if strlowcase(det_type(0)) eq 'lad' then begin drmname = 'inst_response_lad_dec92.*' expon = 2 ew10 = 9. endif else begin ew10 = 12. ;drmname = 'inst_detector_response_spec_jun93.*' drmname = 'inst_detector_response_spec_apr95.*' expon = 1 endelse fdrm=findfile(concat_dir(chklog('SSWDB_BATSE'), drmname)) data=0 num_in_energy=0 file_def=0 @inst_detector_response restore,fdrm(det_id),verb=fcheck(verbose,1) out_frac=file_def.OUT_ENERGY_FRAC edge_products, out_frac, mean=out theta=fcheck(theta,10.) f=cos( (theta/!radeg)^expon + 0.261799) resp=data.response(0,0,*) for i=1,3 do resp = resp+data.response(i,0,*)*f^i resp = reform(resp) ;output x input e_in=data.input_energy if (strlowcase(det_type(0)) eq 'lad' and (fcheck(only_resp) ne -1)) or $ fcheck(only_resp,0) eq 1 then return ;just get the response for LADs gm=fltarr(6,n_elements(e_in)) gm_line=fltarr(3,n_elements(e_in)) gm_esc =gm_line wesc = where( e_in ge 33.2 and e_in le 75.0, nesc ) common funct, function_name function_name = 'gaussian' ;in function_com f_model is used to convert the function for curvefit into a procedure! for i=3,num_in_energy-1 do begin if e_in(i) lt 33. then $ gfit = gaussfit( [out,max(out)+findgen(10)*.1+.1], $ [resp(*,i),fltarr(10)],a) else begin wout = where(out gt 0.9) gfit = gaussfit( [out(wout),max(wout)+findgen(10)*.1+.1], $ [resp(wout,i),fltarr(10)],a) endelse wout = indgen(20)+44 apar = a(0:2) yfit = curvefit(out(wout),resp(wout,i),wout*0+1.,apar,func='funct2') a = apar gm_line(0,i) = apar(0:2) gm(0,i)=apar if i ge min(wesc) and i le max(wesc) then begin ;extract a second peak test = resp(*,i)-gaussian(out,gm(0:2,i)) testfit = gaussfit(out,test>0.,b) gm_esc(0,i) = b(0:2) endif endfor fl=(gm_line(0,*)*gm_line(2,*))(*)*e_in*sqrt(2*!pi) ;line flux fl_esc =(gm_esc(0,*)*gm_esc(2,*))(*)*e_in*sqrt(2*!pi) w=indgen(60)+3 y = (gm(2,*))(*) xx=e_in sigma= [exp(poly( alog(e_in(0:30)),poly_fit(alog(xx(w)),alog(y(w)),3))),$ exp(poly( alog(e_in(31:*)),poly_fit(alog(xx(w)),alog(y(w)),2)))] gm_old = gm gm(1,*) = 1. gm(2,*) = sigma gm(0,*) = f_div(gm_old(0,*) * gm_old(2,*) ,sigma) photo_eff = fl photo_eff(18) = (smooth(fl,5))(18:*) esc_eff = fl_esc wesc_ext= [wesc(0)-2+indgen(2),wesc, max(wesc)+1+indgen(4)] ;interpolate escape peak to slightly higher energy esc_poly_coef= poly_fit(e_in(wesc),esc_eff(wesc), 2, escfit) while escfit(nesc-1) gt escfit(nesc-2) do begin nesc = nesc-1 wesc = wesc(0:nesc-1) esc_poly_coef= poly_fit(e_in(wesc),esc_eff(wesc), 2, escfit) endwhile esc_eff(wesc_ext) = poly(e_in(wesc_ext), esc_poly_coef) > 0.0 ;esc_eff(wesc_ext) = smooth( (spline(e_in(wesc),(fl_esc(wesc)),e_in(wesc_ext),1) >0.0 ) ,5) esc_eff(wesc_ext(0:1)) = 0.0 ;these are below escape energy! ;for response at 10 keV or less, will calculate the response from low_e_spec_eff ;above 10 keV, use photo_eff w10 = where(e_in gt ew10, nw10) drm_resid = resp if strlowcase(det_type(0)) eq 'spec' then $ photo_eff(0:w10(0)) = low_e_spec_eff(e_in(0:w10(0)), theta) ;create final output vectors and arrays! gm_use = gm_old gm_use(2,*) = sigma gm_use(0,*) = photo_eff/(e_in*sqrt(2*!pi))/sigma gm_use(3:5,*) = gm_esc(0:2,*) ;put ones in the escape centroid and width as defaults to avoid arithmetic errors sigma_esc = interpol(sigma, e_in, (e_in*gm_esc(1,*)) >1.) * gm_esc(1,*) gm_use(5,*) = sigma_esc wzero = where(gm_use(5,*) eq 0,nzero) if nzero ge 1 then gm_use(5,wzero)=1.0 gm_use(3,*) = f_div(esc_eff/(e_in*sqrt(2*!pi)),sigma_esc) for i=w10(0), num_in_energy-1 do $ drm_resid(0, i) = resp(*,i) - gaussian(out, gm_use(0:2,i)) - $ gaussian(out, gm_use(3:5,i)) end