Python data and runfile modules

Recent changes:

• mce_script/trunk r960: MCEButterworth class can be used to filter/unfilter data arrays.
• mce_script/trunk r902: MCEFile has options to automatically "unwrap" data bit-windowing effects; also can compute and remove DC filter gain.

The interfaces described here are based on mce_script/trunk r408.

mce_data.py

The module mce_data.py provide roughly the same functionality that mas_data.pro and mas_runfile.pro provide for IDL.

System setup

Make sure python can find mce_data.py (and possibly mce_runfile.py). $MAS_PYTHON should be in your PYTHONPATH environment variable.

mce@mce-ubc-2:~$ export PYTHONPATH=$PYTHONPATH:$MAS_PYTHON

This should be done automatically if you "source $MAS_ROOT/template/mas_env.bash" in your .bashrc".

MCEFile / SmallMCEFile classes

The MCEFile class is an alias for the SmallMCEFile class. These can be used interchangeably. The "Small" is just meant to suggest that you shouldn't try to load more frame data than your machine's RAM can handle. The module will complain and prevent you from loading more than 100 MB of data in one go.

Basic usage

In python, import the 'mce_data' module and create a SmallMCEFile object:

>>> from mce_data import *
>>> fn = '/data/cryo/current_data/data002'
>>> f = SmallMCEFile(fn)

To suppress loading of the runfile, or override the runfile name, use one of these:

>>> f = SmallMCEFile( fn, runfile=False)
>>> f = SmallMCEFIle( fn, runfile='template_runfile.run')

Then load the data like this:

>>> d = f.Read()

The data member of d is a 2d array, [n_detectors, n_frames]:

>>> d.data.shape
(328, 100)

The list of rows and columns associated with the 328 detectors is available through col_list and row_list:

>>> d.col_list
[8, 9, 10, 11, 12, 13, 14, 15, 8, 9, 10, 11, 12, 13, 14, 15,
 8, 9, 10, 11, 12, 13, 14, 15, 8, 9, 10, 11, 12, 13, 14, 15,
 ... 
 8, 9, 10, 11, 12, 13, 14, 15]

To display a pixel's values in all the data frames of a file:

>>> d.data[col+row*8]

Or to see pixel index i:

>>> d.data[i,:]

To display all pixels in frame f:

>>> d.data[:,f]

The data from the first frame header is available in the header member:

>>> d.header
{'status': 2052, 'data_rate': 47, 'userfield': 0, 'num_rows_reported': 41,
 'runfile_id': 1231969044, 'row_len': 64, 'header_version': 6, 'rc_present':
 [False, True, False, False], 'address0_ctr': 23142826, 'num_rows': 41,
 'ramp_value': 0, 'frame_counter': 0, 'sync_box_num': 0, 'ramp_addr': 0}

Data processing

Note that data are automatically rescaled to single normalization associated with each signal type (e.g. unfiltered feedback data are always scaled into DAC units).

To suppress the rescaling:

>>> data_noscale = f.Read(rescale=False)

For filtered feedback data, you may want to automatically compute and remove the DC gain of the low pass filter:

>>> data_dac = f.Read(unfilter='DC')

(The "unfilter" flag has no effect on signals other than the filtered feedback. The unfilter=True option is meant to allow full deconvolution of the filter, but this is not yet implemented.)

Raw frame data

To get raw frames (including header and checksum), pass "raw_frames=True" to the Read call:

>>> raw = f.Read(raw_frames=True)
>>> raw.shape
(100, 1100)

Then you can look at, e.g., the checksum:

>>> raw[:,1099]
>>> dd[:,1099]
array([918158926, 918144283, 918156148, 918156615, 918153934, 918159205,
       ...
       918146379, 918145339, 918158232, 918146376], dtype=int32)

Force row/column ordering

To reformat the d.data array so that its indices are (row, column, frame), pass "row_col=True" to Read:

>>> d = f.Read(row_col=True)
>>> d.data.shape
(33, 32, 100)

Then to see frame k:

>>> d.data[:,:,k]

Extracting fields from mixed mode data

The MCE signal names, for mce_data purposes, are:

error    - co-added error (reference mode 0)
fb       - feedback in sq1 DAC units (reference mode 1)
fb_filt  - filtered feedback (reference mode 2)
fj       - flux jump counter

For mixed modes, the default is to extract the feedback or filtered feedback signal by default. To extract an alternate signal, pass "field=..." to the Read call:

>>> d = f.Read(field='fj')
>>> print d.data
 array([[0,0,0,0,1,0,0,0,0,
 ...

It is possible to extract a set of fields simultaneously, using the "fields=" option in the Read call. This will override the "field=" option. The extracted fields appear in a dictionary in d.data.

>>> d = f.Read(fields=['fj', 'fb_filt'])
>>> print d.data['fj']
    ...
>>> print d.data['fb_filt']
    ...

You can extract all fields by passing "fields='all'":

>>> d = f.Read(fields='all')
>>> print d.data['fj']
    ...
>>> print d.data['fb_filt']
    ...

To force the bitfield extraction code to assume a particular data_mode, pass "data_mode=..." to Read.

>>> d = f.Read(data_mode=10, field='fb_filt')

MCERunfile class

Note that in older versions of mce_script, the MCERunfile class is contained in the file mce_runfile.py. In that case, replace "from mce_data import ..." with "from mce_runfile import ..." in what follows.

Load a runfile

>>> from mce_runfile import *
>>> runfile_name = '/data/cryo/current_data/1220531790_dat.run'
>>> rf = MCERunfile(runfile_name)

Simple data

Recall that the structure of runfiles is such that a line of data has an address defined by its 'block' and 'key' (where the key is the tag + specifiers...). The contents of rf include a dictionary of dictionaries of all the block / key pairs. For example:

>>> print rf.data['HEADER']['RB rc1 data_mode']
 00000010
>>> print rf.data['SQUID']['SQ_tuning_dir']
 1220510497

However, the member function "Item" allows us to repackage the runfile data by specifying a data type ('string', 'int', 'float') and whether or not we expect an array or a single value. For example:

>>> print rf.Item('HEADER', 'RB rc1 data_mode')
['00000010']
>>> print rf.Item('HEADER', 'RB rc1 data_mode', type='int')
[10]
>>> print rf.Item('HEADER', 'RB rc1 data_mode', type='int', array=False)
10
>>> print rf.Item('HEADER', 'RB rc1 data_mode', type='float')
[10.0]

Two-dimensional data

Some runfile entries are really 2d arrays, entered row by row. For example, in the 'IV' block there are per-column entries for responsivity:

<IV>
...
<Responsivity(W/DACfb)_C0> 1.92290e-16 1.92119e-16 0.00000 1.93100e-16 1.92978e-16 1.89769e-16 1.91119e-16 ...
<Responsivity(W/DACfb)_C1> 0.00000 1.82838e-16 0.00000 1.84197e-16 1.84822e-16 1.84447e-16 1.83693e-16 ...
<Responsivity(W/DACfb)_C2> 1.89962e-16 1.88649e-16 0.00000 1.85339e-16 1.84462e-16 1.82965e-16 1.84045e-16 ...
...
</IV>

These can be extracted at once if you pass a printf-style format string to the member function Item2d:

>>> a = rf.Item2d('IV', 'Responsivity(W/DACfb)_C%i', type='float')
>>> print len(a)
32
>>> print len(a[0])
33 
>>> print a[2][1]
1.88649e-16

(i.e. a[2][1] is the responsivity for column 2, row 1.)

Some MCE data are spread across readout cards... e.g. adc_offset is stored in rc# adc_offset#. To accumulate all of these into a single array, use Item2dRC, passing the format in a way such that it can be evaluated first with the RC number and then with the column index... :

>>> a = rf.Item2dRC('HEADER', 'RB rc%i adc_offset%%i', type='int')
>>> len(a)
32
>>> len(a[0])
41

MCEButterworth class

In newer readout card firmware the low pass filter is configurable. The filter parameters can be loaded directly, or parsed from runfiles using the MCEButterworth class.

Parameters may be specified in the constructor; e.g. for a type 1 filter:

>>> import mce_data
>>> filt = mce_data.MCEButterworth([ 32092, 15750, 31238, 14895, 0, 11])

or you can get them from a runfile:

>>> filt = mce_data.MCEButterworth.from_runfile('data_1203.run')

Then you can see the DC response:

>>> filt.gain()
1217.8583042973287

The f3dB frequency:

>>> filt.f3dB(f_samp=15151.)
122.3620800781277

You can also get the whole complex filter transfer function. This example gets the filter gain at frequencies up to 400 Hz, assuming a 15kHz readout:

>>> f = numpy.arange(0., 400)
>>> y = filt.transfer(f, f_samp=15151.)

(The .transfer method used to be called .spectrum.)

To simulate the application of the filter to some time stream data, use apply_filter:

>>> data_freq = 400.
>>> sample_freq = 15151.
>>> data_filtered = filt.apply_filter(data, decimation=data_freq/sample_freq)

To deconvolve the filter from data, use apply_filter with inverse=True:

>>> data_deconvolved = filte.apply_filter(data_filtered, decimation=data_freq/sample_freq, inverse=True)

The application of the FIR to fully sampled data can be simulated with the apply_filter_fir method. See source for documentation.