In [1]:
path = r'D:\CRASYData\Oct14_09.54 - short'              # Directory holding CRASY data files (*.int, *.hdr)
In [2]:
import glob                                             # Utility for file handling
import zlib                                             # Utility for compressing and decompressing data
import numpy as np                                      # Numerical python for array operations
from scipy import sparse, signal                        # Scipy sparse matrix support, signal functions
import matplotlib.pyplot as plt                         # Matplotlib plotting tools
In [3]:
# --- Load mass spectra into sparse matrix
flist = glob.glob(path + '/TS_*[r0|w1].int')            # List alignment data files
datalist = []                                           # Datalist will hold the loaded data
for file in flist:
    with open(file, mode = 'rb') as f:
        data = zlib.decompress(f.read()[10:])           # Load and decompress data
        spec = np.fromstring(data, dtype=np.int16)      # Cast data to 16-bit array
        sspec = sparse.csr_matrix(spec)                 # Convert to sparse data array
        datalist.append(sspec)                          # Append to datalist
csrData = sparse.vstack(datalist, format="csr")         # Sparse matrix of mass versus delay data
delayPts,tofPts = csrData.shape                         # Size of matrix along delay and tof-ms axis
In [4]:
# --- Load delay values from header files.
flist = glob.glob(path + '/TS_*[r0|w1].hdr')      # List header files
delays = np.empty(0)                                    # Array to hold delay values
for file in flist:
    with open(file) as f:
        lines=f.readlines()                             # Read each header file
        delays = np.append(delays, float(lines[1]))     # Keep the delay value (line 2)
delayStep = min(delays[1:] - delays[:-1])               # For sparse data, the smallest delay step should be the nominal stepsize 
In [5]:
# --- For sparse measurements, we must track which delay points were measured.
# --- For continuous measurements, this is unneccessary.
Sparse = False
if Sparse: delayPos = ((delays - delays.min()) / delayStep).astype(int) # Index position of measured delay values
In [6]:
# --- Plot integrated time-of-flight mass spectrum 
tofMS = csrData[0:delayPts:].astype(int).sum(0).A1      # Sum all traces between 0 and delayPts
tofAxis = np.arange(tofPts)*0.5                         # 500 ps binwidth of the digitizer
plt.plot(tofAxis, tofMS)                                # Plot tof-MS
plt.xlabel('time-of-flight (ns)')
plt.show()
In [7]:
# --- Plot integrated delay trace
# --- Change integration boundary for mass-selected traces
trace = csrData[:,0:tofPts].astype(int).sum(1).A1       # Sum all TOF values between 0 and tofPts
plt.plot(delays, trace)                                 # Plot signal as function of time delay
plt.xlabel('delay time (ps)')
plt.show()
In [8]:
# --- Plot Fourier-transformed spectrum
baselinedTrace = signal.detrend(trace)                  # Detrend removes offset and linear trend
if Sparse:                                              # Unsparse trace before FFT
    unsparsed = np.zeros(delayPos[-1]+1)                # create zero-array with length of delay axis
    for value, index in zip(baselinedTrace, delayPos):  # Fill in nonzero values
        unsparsed[index] = value                        # This becomes the unsparsed trace
    baselinedTrace = unsparsed
FFT = np.fft.rfft(baselinedTrace)                       # Fourier transform of real-valued data
freqAxis = np.fft.rfftfreq(baselinedTrace.size, delayStep) # Frequency axis for Fourier transform
plt.plot(freqAxis, abs(FFT)**2)                         # Plot the squared spectrum
plt.xlabel('frequency (THz)')
plt.show()