;--------------------------------------------------------------------------- ; Document name: hsi_visibility_file__define.pro ; Time-stamp: ;--------------------------------------------------------------------------- ; NAME: ; VISIBILITY FILE CLASS DEFINITION ; ; PURPOSE: ; This class provides a processot to hsi_visibility__define ti read ; visibility files. This calss cannot be used alone and needs to be ; accessed through hsi_visibility__define ; ; CATEGORY: ; Utilities hessi/imaging ; ; CONSTRUCTION: ; obj = Obj_New( 'hsi_visibility_file' ) ; (but this is not recommended, use hsi_visibility) ; ; SEE ALSO: ; hsi_visibility__define ; hsi_visibility_raw__define ; ; HISTORY: ; 18-Sep-2019, RAS, In process method, use reverse indices in histogram to organize the pointer array of visibilities after ; reading in the complete bag from the FITS file instead of where,where_arr multiple times. ; For large vis file (117tx15e), old way took ~20 min, new way ~.6 sec. ; 14-Sep-2018, Kim. In process method, added eventlist_strategy and pmtras_diagnostic to tags to remove from ; control. Also added /silent in call to reader->getdata. ; 09-Jul-2018, Kim. In process method, added ev_filename, eb_index, tb_index to tags to remove from control ; structure and cleaned up by putting all tags to remove in one array, rem_tags ; 23-Jun-2017, Kim. In process, if norm_ph_factor tag is not in structure, add it, and set to 0. ; 23-Mar-2017, Kim. In process, remove all tags from control that don't require reprocess ; (previously just removed det_index_mask, now vis_edit, etc.) ; 12-Sep-2012, ras, obtain sp_energy_binning from caller when vis_type is 'reg_electron' ; 17-Dec-2007, Kim. In process, check status after read, and abort if bad ; 13-Dec-2007 -- in process, remove det_index_mask from control struct before setting; ; sept-2005 created, acs ; ;------------------------------------------------------------------- PRO hsi_visibility_file_test o = hsi_visibility() o->set, vis_filename = 'aa.fits' im = o->getdata() obj_destroy, o o = hsi_image() o->set, vis_filename = 'gaga.fits' ;o->set, obs_time= '2005-01-17 09:43:' + ['23', '27'] o->set, im_energy = [25,100] o->plot end ;-------------------------------------------------------------------- FUNCTION HSI_Visibility_File::INIT, caller self.caller = caller RETURN, 1 END ;---------------------------------------------------------- function hsi_visibility_file::need_update ; dont check further down, all visibilities loaded. return, self.caller->get( /need_update, /this ) or $ self.caller->getreload() or $ self.caller->yes_new_eb() or $ self.caller->yes_new_tb() end ;-------------------------------------------------------------------- Function HSI_Visibility_File::Get_Control_Info, reader=reader, control=control, info=info, $ hdr=hdr filename = self.caller->get( /vis_input_fits ) default, control, 1 default, info, 0 control = 1 - info type = (['info','control'])[control] param = hsi_fits_param( filename, reader=reader, type=type, hdr=hdr) return, param end ;-------------------------------------------------------------------- PRO HSI_Visibility_File::Process, $ read=read, $ _EXTRA=_extra cbe = self.caller->get( class = 'hsi_calib_eventlist', /obj ) ; that means it's just a change of detector mask if ~keyword_set(READ) && self.caller->get( /need_update ) eq 0 then return control = Self->Get_control_info(/control, reader=reader, hdr=hdr) ;reader = fitsread() ; ;reader->set, filename = self.caller->get( /vis_input_fits ) ; ;hdr = reader->getheader( /struct, status=status ) ;if not status then begin ; msg = 'Error reading Visibility FITS file.' ; @hsi_message ;endif ; if is_number(control) then return obs_time = [ anytim( hdr.date_obs ), $ anytim( hdr.date_END ) ] self.caller->Set, OBS_TIME = obs_time, /FORCE ;eband = [ hdr.energy_l, hdr.energy_h] ;self.caller->set, energy_band = eband, /FORCE current_vis_type = Self.Caller->Get(/VIS_TYPE) ;control = reader->getdata( extname = 'CONTROL PARAMETERS' ) control = str_top2sub( control ) ;control = rem_tag(control, 'det_index_mask') rem_tags = ['det_index_mask', 'vis_edit', 'vis_chi2lim', 'vis_conjugate', 'vis_normalize', 'vis_plotfit', $ 'vis_type', 'image_algorithm', 'im_energy_binning', 'sp_energy_binning', 'ev_filename', 'eventlist_strategy', $ 'pmtras_diagnostic', 'eb_index', 'tb_index'] control = rem_tag(control, rem_tags) ; this is to fix the old files (pre oct 2008) because we got ; vis_time_intervals and not im_time_intervals if tag_exist( control, 'vis_time_intervals' ) then $ self.caller->set, im_time_interval = control.vis_time_intervals, /FORCE info = Self->Get_control_info(/info, reader=reader, hdr=hdr) ;info = reader->getdata( extname = 'INFO PARAMETERS' ) info = str_top2sub( info ) ;For now we set the vis_type parameter and then set it back to the input value. ;Later, this will be checked for consistency with the bins in visibility::getdata() ;If it can be used for the specified purpose this will work, if not hopefully ;the other software will deal with it gracefully self.caller->set, _extra = info, /FORCE self.caller->set, _extra = control, /FORCE ext = reader->nextextension() vis = reader->getdata(/silent) obj_destroy, reader ;if norm_ph_factor tag is not in structure, add it, and set to 0. if ~tag_exist(vis,'norm_ph_factor') then vis = add_tag(vis,0.,'norm_ph_factor') ;If vis were made for regularized electrons then sp_energy_binning should be used as ;photon vis are computed on these and that's what we're bringing back for now ;set sp_energy_binning ;look for unique first values of erange ;ord = sort(vis.erange[0]) ;vis = vis[ord] ;ix = uniq(vis.erange[0]) ;erange = vis[ix].erange erange = hsi_vis_erange( vis ) obe = self.caller->Get(/obj, class='hsi_binned_eventlist') ;If it's a vis fits file then the binned eventlist must be set to sp_energy_binning for ;the im_energy_binning reconciliation to work obe->set, sp_energy_binning = erange, /force im_energy_binning = erange self.caller->set, im_energy_binning = im_energy_binning, /force ;will pass through reconcile binning and set im_energy as needed ;based on vis_type, i.e. for regularized electrons twice as many im bins as sp bins of the same width an the same energy base ; trange = get_edge_products( self.caller->get( /im_time_interval ), /edges_2 ) ; The commented out code below is the pre 18-sep-2019 version - very slow for large vis files ; ; n_e = n_elements( erange[0,*] ) ; n_t = n_elements( trange[0,*] ) ; nvs = n_elements( vis ) ; vis_arr = ptrarr( n_t, n_e ) ; ; ; for i=0, n_e-1 do begin ; elist = where( total( abs(vis.erange- reproduce(erange[*,i],nvs)),1)/avg(erange[*,i]) lt 1e-5) ; ; ;elist = where_within( vis.erange, erange[*,i] ) ; ; for j=0, n_t-1 do begin ; ; ;tlist = where_within( vis.trange, trange[*,j] ) ; tlist = where( total( abs(vis.trange- reproduce(trange[*,j],nvs)),1) lt 1e-3) ; ; limits = minmax( where_arr( elist, tlist ) ) ; ; ; print, elist[limits[0]] ; ; vis_arr[j,i] = ptr_new( vis[elist[limits[0]:limits[1]]] ) ; ; ; vis_map->setvalues, j,i, elist[limits[0]], limits[1]-limits[0]+1 ; ; endfor ; ; endfor ; Modern version using reverse indices to place all the ti and ei in the correct pointer array bin ; while only sorting once. Much faster than above method. erange0 = get_uniq( erange[0,*]) trange0 = get_uniq( trange[0,*]) n_e = n_elements( erange0 ) n_t = n_elements( trange0 ) vis_arr = ptrarr( n_t, n_e ) ixt = value_locate( trange0, vis.trange[0]) ixe = value_locate( erange0, vis.erange[0]) ixte = (ixt * n_e) + ixe ;every vis has an ixte index, and histogram will find it and set the reverse indices accordingly h = histogram( ixte, rev = rx ) ;t and e are in the wrong order! e should precede t as that's how the data are accumulated vis_arr = ptrarr( n_t, n_e ) for i = 0, n_e-1 do for j = 0, n_t-1 do if h[ j * n_e + i] ge 1 then vis_arr[j,i] = ptr_new( vis[ reverseindices( rx, j * n_e +i)] ) ;self.caller->set, vis_map = vis_map ; x_idx = reader->getpar( 'x_pos' ) ; y_idx = reader->getpar( 'y_pos' ) ; print, x_idx, y_idx ; cbe[x_idx, y_idx ] = ptr_new( data ) ; i = i+1 ;endwhile self.Caller->setdata, vis_arr self.Caller->Set, vis_type = current_vis_type, /force self.caller->set, need_update = 0, /FORCE END ;--------------------------------------------------------------------------- PRO HSI_Visibility_File__Define self = { HSI_Visibility_File, $ caller: obj_new() } END ;--------------------------------------------------------------------------- ; End of 'hsi_visibility_packet__define.pro'. ;---------------------------------------------------------------------------