minor updates in ATHENA frame processing

Ccs/packet_config_ATHENA_DE.py

@@ -179,8 +179,8 @@ class HkPkt:
 
 # Science interface
 
-FRAME_SIZE_X = 512
-FRAME_SIZE_Y = 512
+NPIX_LD = 512
+NPIX_FD = 64
 
 EVT_PKT_ELEMENT_LEN = 6  # length of event packet element (timestamp, header, pixel, footer) in bytes
 STRUCT_EVT_PIX = {  # name, bitsize, bitoffset (LSB)
@@ -368,10 +368,11 @@ class EventFrame:
     FP frame containing event data
     """
 
-    def __init__(self, raw, process=True, notime=False):
+    def __init__(self, raw, framesize=NPIX_LD, process=True, notime=False):
         self._raw = raw
         self._len = len(raw)
         self.evt_list_proc = None
+        self.framesize = framesize
 
         if notime:
             hoff = 1
@@ -414,10 +415,10 @@ class EventFrame:
 
         # initialise empty array
         if nans:
-            arr = np.zeros((FRAME_SIZE_Y, FRAME_SIZE_X))
-            arr[:] = np.NAN
+            arr = np.zeros((self.framesize, self.framesize))
+            arr[:] = np.nan
         else:
-            arr = np.zeros((FRAME_SIZE_Y, FRAME_SIZE_X), dtype=np.uint16)
+            arr = np.zeros((self.framesize, self.framesize), dtype=np.uint16)
 
         for p in self.evt_list_proc:
             arr[p.line_n, p.pixel_n] = p.energy
@@ -476,6 +477,8 @@ class FpmPktParser:
         self.echobytes = echobytes
         self.defaultbytes = defaultbytes
 
+        self.nbytes_tot = 0
+
         self.lastpkt = None
 
     def __call__(self, sock):
@@ -512,6 +515,7 @@ class FpmPktParser:
         msg += self.read_socket(sock, mlen, strict=strict)
 
         self.lastpkt = msg
+        self.nbytes_tot += len(msg)
 
         return msg
 
@@ -544,7 +548,8 @@ class FpmProcessor:
     Processes FPM packets and assembles event frames.
     """
 
-    def __init__(self, ashex=True, process=True, notime=False, tofile=None, queue_to=10):
+    def __init__(self, framesize=NPIX_LD, ashex=True, process=True, notime=False, tofile=None, queue_to=10):
+        self.framesize = framesize
         self.ashex = ashex
         self.process = process
         self.notime = notime
@@ -579,9 +584,9 @@ class FpmProcessor:
 
             # process non-SCI data
             if pkt[0] != IfAddr.SCI:
-                return pkt.hex(sep=' ') if self.ashex else pkt
+                return [pkt.hex(sep=' ')] if self.ashex else [pkt]
 
-            # process data from SCI itnerface
+            # process data from SCI interface
             nevts = (len(pkt) - 1) // EVT_PKT_ELEMENT_LEN
             buf = io.BytesIO(pkt[1:])
             for i in range(nevts):
@@ -603,15 +608,15 @@ class FpmProcessor:
         else:
 
             if self.curfrm != -1:
-                self.frames.append(self.mk_evt_frame(verbose=verbose))
+                self.frames.append(self.mk_evt_frame(framesize=self.framesize, verbose=verbose))
 
             self.cfdata = ed
             self.framecnt += 1
             self.curfrm = ed[0]
 
-    def mk_evt_frame(self, verbose=False):
+    def mk_evt_frame(self, framesize=NPIX_LD, verbose=False):
         try:
-            frame = EventFrame(self.cfdata)
+            frame = EventFrame(self.cfdata, framesize=framesize)
             if verbose:
                 print(frame)
         except Exception as err:
@@ -623,7 +628,7 @@ class FpmProcessor:
     def close_file(self):
         self.tofile.close()
 
-    def proc_file(self, fname, npkts_sci=16, binary=False):
+    def proc_file(self, fname, npkts_sci=16, binary=False, verbose=True):
 
         assert os.path.isfile(fname)
 
@@ -643,7 +648,7 @@ class FpmProcessor:
                 print((pt + msg).hex(sep=' '))
             elif pt[0] == IfAddr.SCI:
                 for i in range(npkts_sci):
-                    self.proc(buf, verbose=True)
+                    self.proc(buf, verbose=verbose)
             else:
                 print('ERROR', pt, buf.tell())
 
@@ -657,6 +662,16 @@ class FpmProcessor:
         self.frames.clear()
 
 
+class FrameList(list):
+
+    def get_frame_ids(self):
+        return [x.frameid for x in self if isinstance(x, EventFrame)]
+
+
+def get_frame_ids(framelist):
+    return [x.frameid for x in framelist if isinstance(x, EventFrame)]
+
+
 def filter_frames(objlist, empty_frames=True):
 
     if empty_frames:

Ccs/tools/dataprocessing/athena_de_frame_proc.py

@@ -28,23 +28,39 @@ def run(fname, binary=False):
     return results
 
 
-def save_frames(frames, outfile, nans=False):
+def save_frames(frames, outfile, nans=False, ascube=False):
 
     from astropy.io import fits
 
     hdl = fits.HDUList()
     empty_frames = 0
 
+    cubedata = []
+
     for frame in frames:
 
         if len(frame.evt_list_proc) == 0:
             empty_frames += 1
             continue
 
-        hdu = fits.ImageHDU(data=frame.as_array(nans=nans))
-        hdu.header['FRAMENUM'] = frame.header.frame_n
-        hdu.header['NEVTS'] = len(frame.evt_list_proc)
+        if ascube:
+            cubedata.append(frame.as_array(nans=nans))
+
+        else:
+            hdu = fits.ImageHDU(data=frame.as_array(nans=nans))
+            hdu.header['FRAMENUM'] = frame.header.frame_n
+            hdu.header['NEVTS'] = len(frame.evt_list_proc)
 
+            hdl.append(hdu)
+
+    if ascube:
+        cube = np.stack(cubedata)
+        hdu = fits.ImageHDU(data=cube)
+        if not nans:
+            nevts = (cube != 0).sum()
+        else:
+            nevts = (~np.isnan(cube)).sum()
+        hdu.header['NEVTSTOT'] = nevts
         hdl.append(hdu)
 
     print('Processed {} frames.'.format(len(frames)))
@@ -69,11 +85,12 @@ class FrameViewer:
         else:
             self.frames = [frame for frame in framelist if frame.nevts > 0]
 
-        self.data = np.zeros((len(self.frames), de.FRAME_SIZE_Y, de.FRAME_SIZE_X))
+        # self.data = np.zeros((len(self.frames), de.FRAME_SIZE_Y, de.FRAME_SIZE_X))
         self.nans = nans
 
-        for i, frame in enumerate(self.frames):
-            self.data[i, :, :] = frame.as_array(nans=self.nans)
+        # for i, frame in enumerate(self.frames):
+            # self.data[i, :, :] = frame.as_array(nans=self.nans)
+        self.data = np.stack([frame.as_array(nans=self.nans) for frame in self.frames])
 
     def show(self, idx=0, cmap='inferno', interpolation='none'):
 
@@ -101,7 +118,10 @@ class FrameViewer:
 
 if __name__ == '__main__':
     if len(sys.argv) < 3:
-        print('USAGE: ./athena_de_frame_proc.py <EVT_FILENAME> <FITS_FILE>')
+        print('USAGE: ./athena_de_frame_proc.py <EVT_FILENAME> <FITS_FILE>\n'
+              '\t-b   read file in binary mode\n'
+              '\t-n   set non-event pixel values to nan\n'
+              '\t-c   collect frames in image cube instead of individual HDUs')
         sys.exit()
 
     if '-b' in sys.argv:
@@ -110,16 +130,30 @@ if __name__ == '__main__':
     else:
         binary = False
 
+    if '-n' in sys.argv:
+        nans = True
+        sys.argv.remove('-n')
+    else:
+        nans = False
+
+    if '-c' in sys.argv:
+        ascube = True
+        sys.argv.remove('-c')
+    else:
+        ascube = False
+
+    # fname = '/home/marko/space/athena/DEdata/dataproc/test_1.txt'
+    # fitsfile = '/home/marko/space/athena/DEdata/dataproc/test_1.fits'
     fname, fitsfile = sys.argv[1:3]
     fitsfile = os.path.abspath(fitsfile)
 
     ow = ''
     if os.path.isfile(fitsfile):
-        while ow.lower() not in ['y', 'n']:
+        while ow.lower().strip() not in ['y', 'n']:
             ow = input('File {} already exists, overwrite? (y/n) '.format(fitsfile))
 
         if ow.lower() == 'n':
             sys.exit()
 
     frames = run(fname, binary=binary)
-    save_frames(frames, fitsfile, nans=False)
+    save_frames(frames, fitsfile, nans=nans, ascube=ascube)