;+ ; ; NAME: ; hxrbs_response ; ; PURPOSE: ; hxrbs response function, count rate per channel as a function of ; incident photon spectrum. ; ; CATEGORY: ; HXRBS, Spectral analysis, fitting ; ; CALLING SEQUENCE: ; ; hxrbs_response, flux=flux, em_flux=em_flux, ut_date=ut, gain=gain, ; edges=edges, newflux=newflux, e_loss=energies, ; omatrix=omatrix, con_factors=con_factors, un_con=un_con, error=error ; ; CALLED BY: ; drm_4_spex ; ; CALLS TO: ; none ; ; Inputs: ; FLUX - photons/cm2/s/keV defined at EM_FLUX, ; optional- only needed to obtain conversion factors ; ; EM_FLUX - energy in keV for FLUX ; optional- default is the energies of the loss matrix ; ; GAIN - HXRBS gain state ; mandatory on the first call, remembered until changed ; ; UT_DATE - Date of observation, any accepted anytim.pro format ; mandatory on the first call, remembered until changed ; RESOLUTION - HXRBS resolution coefficient to use, 0.75 at start of mission ; 1.3 at end ; ; Inputs/Outputs: ; EDGES(2,15) ; - Energy edges of HXRBS 15 channels on date=UT in keV ; These may be given as input values instead of the ; UT_DATE and GAIN ; NEWFLUX(I) ; - Photon flux defined midpoints of EDGES, ; - Interpolated to FLUX and EM_FLUX on midpoints of EDGES ; - if not input ; ; Outputs: ; OMATRIX(I,J) ; - Probability of an input flux in units of photon/cm^2/s in ; - bin J resulting in a rate of counts/cm2/keV/s in bin I ; E_LOSS(J) ; - Energy edges for rows and columns of energy-loss matrix ; CON_FACTORS(I) ; - Conversion factors. Ratio of the (Integrated flux/bin width) ; to NEWFLUX(I) ; UN_CON(I) ; - Uncertainties on the conversion factors. An expression ; of the uncertainty in the detector response ; ERROR - set if a parameter is entered out of range ; ; COMMON BLOCKS: ; hxrbs_response, e_matrix_sav, inmatrix, res_coef_com ; ; SIDE EFFECTS: ; none ; ; RESTRICTIONS: ; none ; ; PROCEDURE: ; ; CALCULATE THE HXRBS RESPONSE MATRIX AND/OR COUNT RATE ; FOR THE DETECTOR GAIN AND RESOLUTION ON THE DATE UT. ; ; Intended use: ; ; This procedure calculates the pulse-height response matrix for an ; input date and gain state. The matrix is computed from an energy-loss matrix ; obtained by running the HXRBS response program, BONOMO, for a set of 100 ; narrow line input functions. The energy loss matrix is broadened for the ; resolution of the detector to yield the pulse-height response matrix for the ; given date of observations. The pulse-height response matrix may be thought ; as the detector pulse-height spectrum in each of the 15 energy channels for an ; incident spectrum of a monoenergetic beam at one of the 100 incident energies. ; The procedure reads in the energy-loss matrix and its associated energies on ; its first call during an IDL session. This energy-loss matrix is resolution ; broadened only when the input date is changed. ; ; This procedure may be used in several ways: ; ; To obtain a pulse-height response matrix for a given date: ; ; HXRBS_RESPONSE, UT_DATE=UT_DATE, GAIN=GAIN, OMATRIX=OMATRIX, E_MATRIX=E_MATRIX, $ ; EDGES=EDGES ; ; This gives the pulse-height response matrix, OMATRIX, the edges of its energy ; input photon bins, E_MATRIX, the edges of the 15 HXRBS channels, EDGES. Using OMATRIX and ; an incident flux integrated over the E_MATRIX bins, one can compute the count ; rate in the 15 HXRBS channels, i.e. ; ; COUNTS (counts/keV/cm2/sec) = OMATRIX#(FLUX (Ph/cm2/s/keV) * (E_MATRIX(1:*)-E_MATRIX)) ; ; This calculation of the count rate is extremely quick after OMATRIX is ; obtained. ; ; CONVERSION FACTOR calculation ; ; CONVERSION FACTORS allow a simple transformation from the photon flux at a ; given energy to the count rate in a given energy channel for a specified ; model spectrum. In general, they are defined as the ratio ; of the count rate spectrum in counts/s/keV/cm^2 to the photon spectrum in ; photons/s/keV/cm^2 provided an incident spectrum is given. Specifically, ; the CONVERSION FACTORS used in the DECONVOLVE section of DCPFIT are the ratio ; of two quantities for a given channel I: ; ; Count rate(counts/sec) in channel I ; ----------------------------------- ; Model Photon Flux at EMID(I) * Area * Width(I) ; ; where EMID(I) is the mid-point energy in channel I, Area is the HXRBS area ; of 68.6 cm^2, and Width(I) is the channel width in keV. ; ; This procedure will also return the CONVERSION FACTORS as follows: ; ; On the first pass: ; ; HXRBS_RESPONSE, UT_DATE=UT_DATE, GAIN=GAIN, E_MATRIX=E_MATRIX, EDGES=EDGES ; ; Then, compute the photon flux, FLUX, on the midpoints of E_MATRIX and ; the photon flux, NEWFLUX, on the midpoints of EDGES outside of this procedure. ; The units of the photon flux must be photons/s/keV/cm^2. ; ; you can get the conversion factors by calling: ; ; HXRBS_RESPONSE, FLUX=FLUX, NEWFLUX=NEWFLUX, $ ; CON_FACTORS=CON_FACTORS, UN_CON=UN_CON ; ; where the 15 CONVERSION FACTORS and their fractional uncertainties are given by ; CON_FACTORS and UN_CON, respectively. ; ; If the energy edges are already known, then these arguments FLUX, NEWFLUX, ; CON_FACTORS, and UN_CON can be entered on the first pass throught the program. ; ; As with the calculation of the response matrix , OMATRIX, detailed above, ; this calculation is quick after the first pass. ; ; The gain state is specified by the function IGAIN(SARS) where SARS is a ; data word found in the DCFILE header. ; ; NON-STANDARD LIBRARY ROUTINES NEEDED: ; ATIME ; CHECKVAR ; CH2KVT ; CH2MV ; EDGE_PRODUCTS ; X_EOUT_DRM ; UTIME ; YMD2SEC ; SSWDB_SMM:HXRBS_MATRIX.DAT - data file with loss matrix ; ; MODIFICATION HISTORY: ; ras, 26-aug-94 ; ras, 20-sep-94 closed the matrix data file after reading ; Kim Tolbert, 19-Feb-2013. Moved from $SSW/smm/hxrbs/idl to ospex dir. Changed location for ; hxrbs_matrix.dat file from $SSWDB_SMM to $SSW_OSPEX/smm_hxrbsresp ;- pro hxrbs_response, flux=flux, em_flux=em_flux, ut_date=ut, gain=gain, $;inputs edges=edges, resolution=res_coeff_in, newflux=newflux, $ e_matrix=e_matrix, omatrix=omatrix, con_factors=con_factors, $ un_con=un_con, area = hxrbs_area, error=error error = 1 ;assume an error common hxrbs_response, e_matrix_sav, inmatrix, res_coef_com res_coeff = fcheck(res_coeff_in, fcheck(res_coef_com, 0.75) ) ;NOMINAL FRACTIONAL UNCERTAINTY IN CONVERSION FACTORS un_con =[14.6, 9.1, 8.1, 7.0, 6.5, 5.6, 5.0, 5.5, 4.9,$ 4.2, 3.2, 3.6, 2.0, 2.2, 2.6]/100. ;THESE VARIABLES ARE INITIALIZED TO MAKE THE PROCEDURE OBTAIN THE ;HXRBS CHANNEL EDGES AND RESOLUTION FROM A SET OF PRIORITIZED CHOICES AND TO ;REMEMBER THEM UNTIL THE CONTROLLING LOGIC CHANGES edges_input = keyword_set(edges) if n_elements(edges) eq 0 then begin ;If edges aren't set then date and gain are needed! ut = fcheck(anytim(ut, /sec),utime('80/6/21')) gain = fcheck(gain,4) ;CHECK FOR A VALID GAIN STATE if total(gain eq (indgen(8)+1)) eq 0 then begin ; gain is not from 1 to 8 error=1 print,'GAIN IS OUT OF RANGE' print,'Input GAIN = ',gain print,'Error: GAIN state should be from 1 to 8' goto, ERROR_RETURN endif ;!!!!!!!!!!!!!!!!!!!!!! ;CHECK FOR A DATE WITHIN THE MISSION LIFETIME ;Mission started 80/02/14 = 35337600. seconds from 79/01/01 ;Mission over by 89/12/01 = 3.44768e+08 seconds from 79/01/01 if (ut lt 35337600.) or (ut gt 3.44768e8) then begin error = 1 print,'DATE IS OUT OF RANGE print,'Input UT = ', atime(ut,/date) print,'Error: HXRBS lasted from 14-Feb-80 - 15-Nov-89' goto, ERROR_RETURN endif RESDAYS= [410, 2418, 2522, 3605, 3624, 3982] RESVALUES= [.77, .88, .95, 1.08, 1.30, 1.30 ] res_coeff = (interpol(resvalues,resdays,ut/86400.))(0) print, 'Using gain state of ',gain,' and time of ',atime(ut,/yoh) print, 'Resolution coefficient = ', res_coeff ;!!!!!!!!!!!!!!!!!!!!!! ;OUTPUT ENERGY BINS ;Get the energy edges for the date UT for the input gainstate from IGAIN(SARS) ; ;IF EDGES ARE NEWLY INPUT, THE DATE OR GAIN HAS CHANGED, ;THEN CHANGE THE 15 HXRBS CHANNEL EDGES ch2kvt, lld=1, gain=gain, ut=ut, e=edges ;CH2KVT checks pars. ;!!!!!!!!!!!!!!!!!!!!!! endif edge_products, edges, mean=emedges, width=wedges, edges_2=edges,edges_1=edges_1 hxrbs_area = 68.6 ;cm2 ;!!!!!!!!!!!!!!!!!!!!!! ;CALCULATE THE RESPONSE MATRIX if n_elements(inmatrix) eq 0 then begin ;read in the basic energy-loss matrix (without resolution broadening) ;ARE WE IN UNIX OR VMS? -> File name depends on operating system energy_loss_file= concat_dir(concat_dir('$SSW_OSPEX', 'smm_hxrbsresp'),'hxrbs_matrix.dat') ;!!!!!!!!!!!!!!!!!!!!!! ;READ IN THE ENERGY LOSS RESPONSE MATRIX IN XDR FORMAT nbins=0 openr,lu1,/get,energy_loss_file,/XDR readu,lu1,nbins e_matrix_sav = fltarr(nbins+1) ;edges of energy-loss matrix inmatrix = fltarr(nbins,nbins) ;energy-loss matrix readu,lu1, e_matrix_sav, inmatrix free_lun, lu1 ;CONVERT MATRIX FROM COUNTS/PHOTON TO COUNTS/(PHOTON/CM2) endif ;MIDPOINT ENERGIES OF ENERGY-LOSS MATRIX BINS edge_products, e_matrix_sav, mean=em_matrix, width=ww_matrix, edges_2=e_matrix ;SET THE DEFAULT INPUT FLUX ENERGIES ;the input matrix energies are the defaults for the input flux ;if em_flux isn't defined it is assumed to be on inmatrix energy mid-points checkvar,em_flux, em_matrix ;!!!!!!!!!!!!!!!!!!!!!! ;!!!!!!!!!!!!!!!!!!!!!! ;BROADEN THE ENERGY-LOSS MATRIX INTO THE MEASURED PULSE-HEIGHT SPECTRUM ;BY CONVOLVING WITH THE RESOLUTION MATRIX. ;COMPUTE THE RESOLUTION BROADENED MATRIX - GMATRIX ;CALCULATE THE PULSE-HEIGHT SPECTRUM FOR THIS RES_COEFF ON THIS DATE sigma = res_coeff * em_matrix ^.75 ;RESOLUTION OF ENERGY LOSS AT EM_MATRIX nmean = n_elements(em_matrix) pulse_shape = fltarr(nmean,nmean) for i=0, nmean-1 do $ pulse_shape(*,i) = gaussint( (e_matrix(1,*)-em_matrix(i))/sigma(i)) - $ gaussint( (e_matrix(0,*)-em_matrix(i))/sigma(i)) gmatrix = pulse_shape#inmatrix ;GMATRIX IS THE MEASURED PULSE-HEIGHT MATRIX IN COUNTS/PHOTON ;DEFINED ON THE BINS SPECIFIED BY E_MATRIX, ;Define gmatrix in counts/kev/incident photon gmatrix = temporary( gmatrix/rebin(reform(ww_matrix,nmean,1),nmean,nmean)) ;!!!!!!!!!!!!!!!!!!!!!! ;INTERPOLATE THE PULSE-HEIGHT MATRIX TO THE FINAL BINS omatrix = x_eout_drm( in_drm=gmatrix, in_e_out=e_matrix, out_e_out=edges, nsub=10) ;!!!!!!!!!!!!!!!!!!!!!! ;IF A PHOTON FLUX HAS BEEN INPUT, APPLY THE RESPONSE MATRIX AND COMPUTE ;CONVERSION FACTORS. IF NOT, THEN JUST RETURN THE RESPONSE MATRIX. if n_elements(flux) ne 0 then begin ;IF THE EM_FLUX ARE DIFFERENT FROM THE EM_MATRIX, INTERPOLATE FLUX TO EM_MATRIX wdiff = where(em_matrix - em_flux ne 0.0, ndiff) if total(abs(em_matrix-em_flux)) ne 0. then flux2 = interpol( flux, em_flux, em_matrix) else $ flux2 = flux ;FROM PHOTON/CM2/SEC/KEV TO PHOTON/CM2/SEC phot_per_bin = flux2 * ww_matrix ;CALCULATE THE COUNT RATE IN THE 15 HXRBS CHANNELS cnts_per_kevcm2 = omatrix#phot_per_bin ;EXPRESS IN UNITS OF COUNTS/CM2/SEC/KEV TO GET CONVERSION FACTORS newflux = interpol(flux,em_flux,emedges) con_factors = cnts_per_kevcm2 / newflux endif error = 0 ;success! ERROR_RETURN: return end