path = r'D:\CRASYData\Oct14_09.54 - short' # Directory holding CRASY data files (*.int, *.hdr)
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
# --- 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
# --- 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
# --- 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
# --- 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()
# --- 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()
# --- 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()