wibethdecoder.py

Go to the documentation of this file.

#!/usr/bin/env python3
 
from hdf5libs import HDF5RawDataFile
 
import daqdataformats
import detdataformats
import fddetdataformats
from rawdatautils.unpack.wibeth import *
from rawdatautils.utilities.wibeth import *
import detchannelmaps
 
import click
import time
import numpy as np
 
 
@click.command()
@click.argument('filenames', nargs=-1, type=click.Path(exists=True))
@click.option('--nrecords', '-n', default=-1, help='How many Trigger Records to process (default: all)')
@click.option('--nskip', default=0, help='How many Trigger Records to skip (default: 0)')
@click.option('--channel-map', default=None, help="Channel map to load (default: None)")
@click.option('--print-headers', is_flag=True, help="Print WIB2Frame headers")
@click.option('--print-adc-stats', is_flag=True, help="Print ADC Pedestals/RMS")
@click.option('--check-timestamps', is_flag=True, help="Check WIB2 Frame Timestamps")
@click.option('--det', default='HD_TPC', help='Subdetector string (default: HD_TPC)')
@click.option('--quiet', is_flag=True, help="Print only final summary information")
 
def main(filenames, nrecords, nskip, channel_map, print_headers, print_adc_stats, check_timestamps, det, quiet):
 
    for filename in filenames:
    
        h5_file = HDF5RawDataFile(filename)
 
        records = h5_file.get_all_record_ids()
 
        if nskip > len(records):
            print(f'Requested records to skip {nskip} is greater than number of records {len(records)}. Exiting...')
            return
        if nrecords>0:
            if (nskip+nrecords)>len(records):
                nrecords=-1
            else:
                nrecords=nskip+nrecords
 
        records_to_process = []
        if nrecords==-1:
            records_to_process = records[nskip:]
        else:
            records_to_process = records[nskip:nrecords]
        print(f'Will process {len(records_to_process)} of {len(records)} records in {filename}.')
 
        #have channel numbers per geoid in here
        ch_map = None
        if channel_map is not None:
            ch_map = detchannelmaps.make_tpc_map(channel_map)
        offline_ch_num_dict = {}
        offline_ch_plane_dict = {}
 
        for r in records_to_process:
 
            if not quiet:
                print(f'Processing (Record Number,Sequence Number)=({r[0],r[1]})')
 
            wib_geo_ids = h5_file.get_geo_ids_for_subdetector(r,detdataformats.DetID.string_to_subdetector(det))
 
            for gid in wib_geo_ids:
                #geo_info = detchannelmaps.HardwareMapService.parse_geo_id(gid)
                det_link = 0xffff & (gid >> 48);
                det_slot = 0xffff & (gid >> 32);
                det_crate = 0xffff & (gid >> 16);
                det_id = 0xffff & gid;
                subdet = detdataformats.DetID.Subdetector(det_id)
                det_name = detdataformats.DetID.subdetector_to_string(subdet)
            
                if not quiet:
                    print(f'\tProcessing subdetector {det_name}, crate {det_crate}, slot {det_slot}, link {det_link}')
 
                frag = h5_file.get_frag(r,gid)
                frag_hdr = frag.get_header()
                frag_type = frag.get_fragment_type()
                frag_ts = frag.get_trigger_timestamp()
                frag_window_begin = frag.get_window_begin()
                frag_window_end = frag.get_window_end()
                if(frag_type!=daqdataformats.FragmentType.kWIBEth):
                    print('\tNot WIBEth fragment type {frag_type}. Continue.')
                    continue
 
                if not quiet:
                    print(f'\tTrigger timestamp for fragment is {frag_ts}')
 
                n_frames = get_n_frames(frag);
 
                if not quiet:
                    print(f'\tFound {n_frames} WIBEth Frames.')
 
                wf = fddetdataformats.WIBEthFrame(frag.get_data())
 
                #print header info
                if print_headers:
                    print('\n\t==== WIB HEADER (First Frame) ====')
                    print_header(wf,prefix='\t\t')
 
                #fill channel map info if needed
                if(offline_ch_num_dict.get(gid) is None):
                    if channel_map is None:
                        offline_ch_num_dict[gid] = np.arange(64)
                        offline_ch_plane_dict[gid] = np.full(64,9999)
                    else:
                        dh = wf.get_daqheader()
                        offline_ch_num_dict[gid] = np.array([ch_map.get_offline_channel_from_det_crate_slot_stream_chan(dh.det_id, dh.crate_id, dh.slot_id, dh.stream_id, c) for c in range(64)])
                        offline_ch_plane_dict[gid] = np.array([ ch_map.get_plane_from_offline_channel(uc) for uc in offline_ch_num_dict[gid] ])
 
 
                #unpack timestamps into numpy array of uin64
                if check_timestamps:
                    timestamps = np_array_timestamp(frag)
                    timestamps_diff = np.diff(timestamps)
                    timestamps_diff_vals, timestamps_diff_counts = np.unique(timestamps_diff, return_counts=True)
                    
                    if(n_frames>0):
                        if not quiet or len(timestamps_diff_counts)>1:
                            print('\n\t==== TIMESTAMP CHECK ====')
                            print(f'\t\tTimestamps (First, Last, Min, Max): ({timestamps[0]},{timestamps[-1]},{np.min(timestamps)},{np.max(timestamps)})')
                            print(f'\t\tTimestamp diffs: {timestamps_diff_vals}')
                            print(f'\t\tTimestamp diff counts: {timestamps_diff_counts}')
                            print(f'\t\tAverage diff: {np.mean(timestamps_diff)}')
                            data_extent_before_readout_window = frag_window_begin - timestamps[0]
                            data_extent_after_readout_window = timestamps[-1] - frag_window_end
                            print(f'\t\tData extent beyond readout window (before,after): {data_extent_before_readout_window},{data_extent_after_readout_window} (approximate)')
                            if data_extent_before_readout_window < 0.0 or data_extent_after_readout_window < -32.0:
                                print(f'\t\t\tWARNING: the WIBEth data does not fully cover the readout window ({frag_window_begin},{frag_window_end})')
 
 
                if print_adc_stats:
 
                    #unpack adcs into a n_frames x 256 numpy array of uint16
                    adcs = np_array_adc(frag)
                    adcs_rms = np.std(adcs,axis=0)
                    adcs_ped = np.mean(adcs,axis=0)
                
                    print('\n\t====WIB DATA====')
 
                    for ch,rms in enumerate(adcs_rms):
                        print(f'\t\tch {offline_ch_num_dict[gid][ch]} (plane {offline_ch_plane_dict[gid][ch]}): ped = {adcs_ped[ch]:.2f}, rms = {adcs_rms[ch]:.4f}')
 
                #print("\n")
            #end gid loop
 
            if check_timestamps:
                timestamps_frame0 = np.array([ fddetdataformats.WIBEthFrame(h5_file.get_frag(r,gid).get_data()).get_timestamp() for gid in wib_geo_ids ])
                timestamps_frame0_diff = timestamps_frame0 - timestamps_frame0[0]
 
                if not quiet or np.any(timestamps_frame0_diff):
                    print('\n\t==== TIMESTAMP ACROSS WIBS CHECK ====')
                    print(f'\t\tTimestamp diff relative to first WIB',timestamps_frame0_diff)
        
        #end record loop
 
    #end file loop
    print(f'Processed all requested records')
 
    
if __name__ == '__main__':
    main()