! Program to fake up a response matrix. Assumes that the response is ! Gaussian and that the efficiency is given by a combination of data ! read in from effective area, detector efficiency, and filter transmission ! files. The last two of these can be ignored by setting the filename to ! 'none'. If the spectrometer is defined as dispersive then parameters are ! assumed to be given in Angstroms and if non-dispersive then in keV. The ! response bins and channels are linear with hidden parametes allowing a ! single break to a different linear relation. There are a number of options ! for the change in resolution with energy or wavelength : constant is ! constant in either energy or wavelength; sqroot is square root in energy; ! czt is the relation for CdZnTe detectors; and file reads the relation ! from a file. ! v1.00 kaa 1/22/96 ! v1.01 kaa 4/17/98 bug fix in stbins. was not setting the number of ! channels and energies correctly when there was no ! break in the channel/energy binsize. ! v1.10 kaa 2/16/98 Added the option to give an input RMF file instead ! of calculating it from the resolution information ! v1.11 ngan 9/21/00 Initialized the variable ierr to 0. ! v1.12 kaa 7/17/02 Fixed unsaved variables and incorrect initialization in getsig ! v1.13 kaa 7/30/04 Improved error if user gives an invalid value for ! the res_reln parameter. ! v1.14 kaa 2/02/05 Corrected bug when reading from a file resolution ! information for a dispersive instrument. This option ! has probably never been used. ! v1.15 kaa 5/17/05 Added the option to include escape peaks. ! v2.00 kaa 7/25/05 Rewrite to allow more general gain relations ! v2.01 kaa 7/12/06 Fixed bug (introduced in v2.00) in generating gaussian ! response shape ! v2.02 kaa 4/02/08 Added LINEAR option for resolution. ! v2.03 kaa 7/09/08 Fixed bug in generating gaussian response shape when channel ! boundaries are in decreasing energy order (viz gratings) ! v2.04 kaa 8/05/08 Write LO_THRES keyword as well as LO_THRESH since the ! former actually conforms to the FITS standard. ! v2.05 kaa 11/16/10 Improved diagnostics for user error. ! v2.06 kaa 4/26/11 Fixed a bug when reading energies from resp_file. ! v2.07 kaa 5/2/14 Ensure that E_MIN is always <= E_MAX. ! v2.08 kaa 5/15/15 Fixed bug when the response energies extend beyond the ! channel energies - had to fix bisearch in clcrsp.f ! v2.09 MFC 4/16/20 f90 version SUBROUTINE genrsp IMPLICIT NONE ! Pointers for dynamic arrays ! resp_matrix(max_elements) the response matrix elements ! resp_energies(0:n_energies) the response matrix energy ranges ! eff_area(2,n_a_energies) the energies and effective areas ! filt_eff(2,n_f_energies) the energies and filter transmissions ! det_eff(2,n_d_energies) the energies and detector efficiencies ! lineEn(n_r_energies) the energies for tabulated response data ! lineDt(3,n_lines,n_r_energies) the centroids (relative to energy), ! resolutions, and normalizations ! of lines ! resp_data(2,n_e_energies) definition of the response energies ! chan_data(2,n_c_energies) definition of the channels ! ch_bounds(n_channels,2) the channel boundary energies ! ngroup(n_energies) the number of response groups at ! each energy ! ichanb(max_tot_groups) start channel for response groups ! ichane(max_tot_groups) end channel for response groups ! INTEGER*8 resp_matrix, resp_energies, eff_area, filt_eff, det_eff ! INTEGER*8 resp_work, lineEn, lineDt, resp_data, chan_data ! INTEGER*8 ch_bounds, ngroup, ichanb, ichane character(6) version character(10) verdate parameter(version = '2.09', verdate='04/16/2020') INTEGER (kind=8) :: resp_matrix, resp_energies, eff_area, filt_eff, det_eff INTEGER (kind=8) :: resp_work, lineEn, lineDt, resp_data, chan_data INTEGER (kind=8) :: ch_bounds, ngroup, ichanb, ichane ! Other variables REAL fwhm, rsp_min, resp_low, resp_high, resp_break REAL chan_low, chan_high, chan_break INTEGER resp_number, resp_bnumber, chan_number, chan_bnumber INTEGER num_elements, max_elements, max_tot_groups INTEGER num_tot_groups, n_a_energies, n_channels, n_d_energies INTEGER n_energies, n_f_energies, n_r_energies, n_lines INTEGER n_e_energies, n_c_energies, ierr CHARACTER(20) resp_reln, chan_reln, resol_reln CHARACTER(72) resp_file, chan_file, resol_file CHARACTER(20) tlscpe, instrm CHARACTER(72) contxt, rmffil, efffil, detfil, filfil, inrfil LOGICAL disperse, clobber INTEGER lenact EXTERNAL lenact ! common block for dynamic memory INTEGER MEMI(100) REAL MEMR(100) EQUIVALENCE (MEMI, MEMR) COMMON /MEM/ MEMR CHARACTER(40) taskname COMMON /task/ taskname ierr = 0 CALL fcecho('GENRSP vers '//version//' '//verdate) taskname = 'genrsp '//version ! get the parameters CALL RSPPAR(inrfil, rmffil, disperse, tlscpe, instrm, & resp_reln, resp_file, resp_low, resp_high, & resp_number, resp_break, resp_bnumber, chan_reln, & chan_file, chan_low, chan_high, chan_number, & chan_break, chan_bnumber, efffil, detfil, filfil, & resol_reln, resol_file, fwhm, max_elements, rsp_min, & clobber, ierr) contxt = 'Failed to read parameters' IF (ierr.NE.0) GOTO 999 CALL fcecho(' ') ! get the sizes of the input data (effective area, detector efficiency, ! filter transmission, and if requested resolution, response energy, and ! channel data) CALL GTDTSZ(efffil, detfil, filfil, resol_file, resp_file, & chan_file, n_a_energies, n_d_energies, n_f_energies, & n_r_energies, n_e_energies, n_c_energies, n_lines, & ierr) contxt = 'Failed to get sizes of input files' IF (ierr.NE.0) GOTO 999 ! grab the memory for the input data eff_area = 0 filt_eff = 0 det_eff = 0 lineEn = 0 lineDt = 0 resp_data = 0 chan_data = 0 IF (n_a_energies .GT. 0) THEN CALL UDMGET(2 * n_a_energies, 6, eff_area, ierr) contxt = ' Not enough memory for eff_area array' IF (ierr.NE.0) GOTO 999 ENDIF IF (n_f_energies .GT. 0) THEN CALL UDMGET(2 * n_f_energies, 6, filt_eff, ierr) contxt = ' Not enough memory for filt_eff array' IF (ierr.NE.0) GOTO 999 ENDIF IF (n_d_energies .GT. 0) THEN CALL UDMGET(2 * n_d_energies, 6, det_eff, ierr) contxt = ' Not enough memory for det_eff array' IF (ierr.NE.0) GOTO 999 ENDIF IF (n_r_energies .GT. 0) THEN CALL UDMGET(n_r_energies, 6, lineEn, ierr) contxt = ' Not enough memory for lineEn array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(3 * n_lines * n_r_energies, 6, lineDt, ierr) contxt = ' Not enough memory for lineDt array' IF (ierr.NE.0) GOTO 999 ENDIF IF (n_e_energies .GT. 0) THEN CALL UDMGET(2 * n_e_energies, 6, resp_data, ierr) contxt = ' Not enough memory for resp_data array' IF (ierr.NE.0) GOTO 999 ENDIF IF (n_c_energies .GT. 0) THEN CALL UDMGET(2 * n_c_energies, 6, chan_data, ierr) contxt = ' Not enough memory for chan_data array' IF (ierr.NE.0) GOTO 999 ENDIF ! read in the effective areas, detector efficiencies, filter ! transmissions, lines (energies, resolutions, norms), and ! gain relation CALL RDINPD(n_a_energies, efffil, MEMR(eff_area), n_d_energies, & detfil, MEMR(det_eff), n_f_energies, filfil, & MEMR(filt_eff), n_r_energies, n_lines, resol_file, & MEMR(lineEn), MEMR(lineDt), n_e_energies, resp_file, & MEMR(resp_data), n_c_energies, chan_file, & MEMR(chan_data), ierr) contxt = 'Failed to read input data' IF (ierr.NE.0) GOTO 999 ! If no input RMF file was given then we need to make one so first ! calculate or estimate the sizes of the arrays that will be required. IF (inrfil .EQ. 'none' .OR. inrfil .EQ. 'NONE') THEN ! calculate the number of response energy bins and channels ! if the response energy definition file was given then use that ! number IF (n_e_energies .GT. 0) THEN n_energies = n_e_energies ! otherwise calculate from the other parameters ELSE n_energies = Resp_number IF (Resp_high .GT. Resp_break .AND.& Resp_break .GT. 0.) THEN n_energies = n_energies + Resp_bnumber ENDIF ENDIF ! if the response energy definition file was given then use that ! number IF (n_e_energies .GT. 0) THEN n_energies = n_e_energies ! otherwise calculate from the other parameters ELSE n_energies = Resp_number IF (Resp_high .GT. Resp_break .AND.& Resp_break .GT. 0.) THEN n_energies = n_energies + Resp_bnumber ENDIF ENDIF ! if the channel definition file was given then use that ! number IF (n_c_energies .GT. 0) THEN n_channels = n_c_energies ! otherwise calculate from the other parameters ELSE n_channels = Chan_number IF (Chan_high .GT. Chan_break .AND.& Chan_break .GT. 0.) THEN n_channels = n_channels + Chan_bnumber ENDIF ENDIF CALL xwrite(' ', 10) WRITE(contxt, '(a,i7,a)') '...', n_channels, & ' channels in spectrum' CALL xwrite(contxt, 10) WRITE(contxt, '(a,i7,a)') '...', n_energies, & ' energies in response' CALL xwrite(contxt, 10) ! estimate the total number of response groups max_tot_groups = n_energies * 10 ! in this case when an input RMF is given then we can read that ! to get the array sizes ELSE CALL GTRSSZ(inrfil, max_elements, n_energies, n_channels, & max_tot_groups, ierr) contxt = 'Failed to get sizes of input rmf data' IF (ierr.NE.0) GOTO 999 ENDIF ! now grab the memory for the response matrix arrays resp_matrix = 0 resp_energies = 0 resp_work = 0 ch_bounds = 0 ngroup = 0 ichanb = 0 ichane = 0 CALL UDMGET(max_elements, 6, resp_matrix, ierr) contxt = ' Not enough memory for resp_matrix array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(n_energies + 1, 6, resp_energies, ierr) contxt = ' Not enough memory for resp_energies array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(n_channels, 6, resp_work, ierr) contxt = ' Not enough memory for resp_work array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(2 * n_channels, 6, ch_bounds, ierr) contxt = ' Not enough memory for ch_bounds array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(n_energies, 4, ngroup, ierr) contxt = ' Not enough memory for ngroup array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(max_tot_groups, 4, ichanb, ierr) contxt = ' Not enough memory for ichanb array' IF (ierr.NE.0) GOTO 999 CALL UDMGET(max_tot_groups, 4, ichane, ierr) contxt = ' Not enough memory for ichane array' IF (ierr.NE.0) GOTO 999 ! Again if there was no input RMF then calculate it IF (inrfil .EQ. 'none' .OR. inrfil .EQ. 'NONE') THEN ! set up channel boundaries and response energies CALL STBINS(n_channels, chan_low, chan_high, chan_break, & chan_number, chan_bnumber, n_energies, resp_low, & resp_high, resp_break, resp_number, resp_bnumber, & disperse, resp_reln, n_e_energies, MEMR(resp_data), & chan_reln, n_c_energies, MEMR(chan_data), & MEMR(ch_bounds), MEMR(resp_energies)) ! calculate the response matrix CALL CLCRSP(n_energies, MEMR(resp_energies), n_channels, & MEMR(ch_bounds), rsp_min, MEMI(ngroup), & max_tot_groups, max_elements, disperse, resol_reln, & fwhm, n_r_energies, n_lines, MEMR(lineEn), & MEMR(lineDt), MEMI(ichanb), MEMI(ichane), & MEMR(resp_matrix), MEMR(resp_work), num_elements, & num_tot_groups, ierr) contxt = 'Failed to calculate response matrix' IF (ierr.NE.0) GOTO 999 IF (num_elements .EQ. 0) THEN contxt = & 'The calculated response has no elements greater than rsp_min' CALL fcecho(contxt) CALL exit(1) ENDIF ! or if there was an input RMF then read the information from that ELSE CALL INPRMF(inrfil, n_energies, MEMR(resp_energies), & n_channels, MEMR(ch_bounds), rsp_min, & MEMI(ngroup), max_tot_groups, max_elements, & MEMI(ichanb), MEMI(ichane), MEMR(resp_matrix), & num_elements, num_tot_groups, ierr) contxt = 'Failed to input response matrix' IF (ierr.NE.0) GOTO 999 ENDIF ! fold effective areas into response matrix IF (n_a_energies .GT. 0) THEN CALL FOLDIN(max_elements, n_energies, n_a_energies, & max_tot_groups, MEMR(eff_area), MEMR(resp_matrix), & MEMR(resp_energies), MEMI(ngroup), MEMI(ichanb), & MEMI(ichane), num_elements, num_tot_groups) ENDIF ! fold detector efficiencies into response matrix IF (n_d_energies .GT. 0) THEN CALL FOLDIN(max_elements, n_energies, n_d_energies, & max_tot_groups, MEMR(det_eff), MEMR(resp_matrix), & MEMR(resp_energies), MEMI(ngroup), MEMI(ichanb), & MEMI(ichane), num_elements, num_tot_groups) ENDIF ! fold filter transmissions into response matrix IF (n_f_energies .GT. 0) THEN CALL FOLDIN(max_elements, n_energies, n_f_energies, & max_tot_groups, MEMR(filt_eff), MEMR(resp_matrix), & MEMR(resp_energies), MEMI(ngroup), MEMI(ichanb), & MEMI(ichane), num_elements, num_tot_groups) ENDIF ! write response matrix CALL WRTRSP(max_elements, n_energies, n_channels, max_tot_groups, & MEMR(resp_matrix), MEMR(resp_energies), & MEMI(ngroup), MEMI(ichanb), MEMI(ichane), & MEMR(ch_bounds), rmffil, tlscpe, instrm, rsp_min, & clobber, ierr) contxt = ' Failed to write response matrix' IF (ierr .NE. 0) GOTO 999 999 CONTINUE IF (ierr.NE.0) THEN CALL xwrite(contxt, 10) WRITE(contxt, '(a,i4)') ' error = ', ierr CALL xwrite(contxt, 10) ENDIF CALL EXIT(0) END