Coverage for larch/io/xrf_netcdf.py: 12%
121 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1#!/usr/bin/python
2"""
3support for netcdf file output files containing MCA spectra
4from Epics Mapping Mode with XIA xMXAP electronics
5"""
6import numpy as np
7import time
8import sys
9import os
11from scipy.io import netcdf_file
13def aslong(d):
14 """unravels and converts array of int16 (int) to int32 (long)"""
15 # need to unravel the array!!!
16 d = d.astype(np.int16).ravel()
17 d.dtype = np.int32
18 return d
20class xMAPBufferHeader(object):
21 def __init__(self,buff):
22 self.tag0 = buff[0] # Tag word 0
23 self.tag1 = buff[1] # Tag word 1
24 self.headerSize = buff[2] # Buffer header size
25 # Mapping mode (1=Full spectrum, 2=Multiple ROI, 3=List mode)
26 self.mappingMode = buff[3]
27 self.runNumber = buff[4] # Run number
28 # Sequential buffer number, low word first
29 self.bufferNumber = aslong(buff[5:7])[0]
30 self.bufferID = buff[7] # 0=A, 1=B
31 self.numPixels = buff[8] # Number of pixels in buffer
32 # Starting pixel number, low word first
33 self.startingPixel = aslong(buff[9:11])[0]
34 self.moduleNumber = buff[11]
35 self.channelID = np.array(buff[12:20]).reshape((4,2))
36 self.channelSize = buff[20:24]
37 self.bufferErrors = buff[24]
38 self.userDefined = buff[32:64]
40class xMAPMCAPixelHeader(object):
41 def __init__(self,buff):
42 self.tag0 = buff[0]
43 self.tag1 = buff[1]
44 self.headerSize = buff[2]
45 # Mapping mode (1=Full spectrum, 2=Multiple ROI, 3=List mode)
46 self.mappingMode = buff[3]
47 self.pixelNumber = aslong(buff[4:6])[0]
48 self.blockSize = aslong(buff[6:8])[0]
50class xMAPData(object):
51 def __init__(self,npix,nmod,nchan):
52 ndet = 4 * nmod
53 self.firstPixel = 0
54 self.numPixels = 0
55 self.counts = np.zeros((npix, ndet, nchan), dtype='i2')
56 self.realTime = np.zeros((npix, ndet), dtype='i8')
57 self.liveTime = np.zeros((npix, ndet), dtype='i8')
58 self.inputCounts = np.zeros((npix, ndet), dtype='i4')
59 self.outputCounts = np.zeros((npix, ndet), dtype='i4')
61CLOCKTICK = 0.320 # xmap clocktick = 320 ns
63def read_xrf_netcdf(fname, npixels=None, verbose=False):
64 # Reads a netCDF file created with the DXP xMAP driver
65 # with the netCDF plugin buffers
66 if verbose:
67 print( ' reading ', fname)
68 t0 = time.time()
69 # read data from array_data variable of netcdf file
70 read_ok = False
71 fh = None
72 try:
73 fh = netcdf_file(fname,'r')
74 read_ok = True
75 except:
76 time.sleep(0.010)
77 try:
78 fh = netcdf_file(fname,'r')
79 read_ok = True
80 except:
81 pass
83 if not read_ok or fh is None:
84 if fh is not None:
85 fh.close()
86 return None
88 array_data = fh.variables['array_data']
89 t1 = time.time()
91 # array_data will normally be 3d:
92 # shape = (narrays, nmodules, buffersize)
93 # but nmodules and narrays could be 1, so that
94 # array_data could be 1d or 2d.
95 #
96 # here we force the data to be 3d
97 shape = array_data.shape
98 if len(shape) == 1:
99 array_data.shape = (1, 1, shape[0])
100 elif len(shape) == 2:
101 array_data.shape = (1, shape[0], shape[1])
103 narrays, nmodules, buffersize = array_data.shape
104 modpixs = int(max(124, array_data[0, 0, 8]))
105 npix_total = 0
106 # real / live times are returned in microseconds.
107 for array in range(narrays):
108 for module in range(nmodules):
109 d = array_data[array,module, :]
110 bh = xMAPBufferHeader(d)
111 #if verbose and array==0:
112 dat = d[256:].reshape(modpixs, int((d.size-256)/modpixs ))
114 npix = bh.numPixels
115 if module == 0:
116 npix_total += npix
117 if array == 0:
118 # first time through, (array,module)=(0,0) we
119 # read mapping mode, set up how to slice the
120 # data, and build data arrays in xmapdat
121 mapmode = dat[0, 3]
122 if mapmode == 1: # mapping, full spectra
123 nchans = d[20]
124 data_slice = slice(256, 8448)
125 elif mapmode == 2: # ROI mode
126 # Note: nchans = number of ROIS !!
127 nchans = max(d[264:268])
128 data_slice = slice(64, 64+8*nchans)
129 xmapdat = xMAPData(narrays*modpixs, nmodules, nchans)
130 xmapdat.firstPixel = bh.startingPixel
132 # acquistion times and i/o counts data are stored
133 # as longs in locations 32:64
134 t_times = aslong(dat[:npix, 32:64]).reshape(npix, 4, 4)
135 p1 = npix_total - npix
136 p2 = npix_total
137 xmapdat.realTime[p1:p2, :] = t_times[:, :, 0]
138 xmapdat.liveTime[p1:p2, :] = t_times[:, :, 1]
139 xmapdat.inputCounts[p1:p2, :] = t_times[:, :, 2]
140 xmapdat.outputCounts[p1:p2, :] = t_times[:, :, 3]
142 # the data, extracted as per data_slice and mapmode
143 t_data = dat[:npix, data_slice]
144 if mapmode == 2:
145 t_data = aslong(t_data)
146 xmapdat.counts[p1:p2, :, :] = t_data.reshape(npix, 4, nchans)
148 t2 = time.time()
149 xmapdat.numPixels = npix_total
150 xmapdat.counts = xmapdat.counts[:npix_total]
151 xmapdat.realTime = CLOCKTICK * xmapdat.realTime[:npix_total]
152 xmapdat.liveTime = CLOCKTICK * xmapdat.liveTime[:npix_total]
153 xmapdat.inputCounts = xmapdat.inputCounts[:npix_total]
154 xmapdat.outputCounts = xmapdat.outputCounts[:npix_total]
155 if verbose:
156 print(' time to read file = %5.1f ms' % ((t1-t0)*1000))
157 print(' time to extract data = %5.1f ms' % ((t2-t1)*1000))
158 print(' read %i pixels ' % npix_total)
159 print(' data shape: ' , xmapdat.counts.shape)
160 fh.close()
161 return xmapdat
163def test_read(fname):
164 print( fname, os.stat(fname))
165 fd = read_xrf_netcdf(fname, verbose=True)
166 print(fd.counts.shape)