import re
import mmap
import struct

with open('Zn_run_600ms_61Zn_011.dat', "r+b") as f:
    # pattern is a string that helps locating the binary part in a Eva File
    pattern = b'*----------------here the binary part starts --------------*'
    # memory-map the file, size 0 means whole file
    mm = mmap.mmap(f.fileno(), 0, access=0)
    # Get MCA Parameters from the Header
    mca_pattern = {'section': b'[MCA]', 'micro_sec': b'\xb5s,'}

    selection_from = mm.find(mca_pattern['micro_sec'])
    selection_to = mm.find(mca_pattern['micro_sec']) + 20
    mca_ch = int(re.findall(b'[0-9]+', mm[selection_from: selection_to])[0])

    selection_from = mm.find(mca_pattern['section']) + 7
    selection_to =mm.find(mca_pattern['section']) + 70
    mca_fact = float(re.findall(b'[0-9]+\.[0-9]+',mm[selection_from: selection_to])[0])
    # Get Main Parameters from the the X- and Y-Scan Array
    binary_start = mm[7:].find(pattern) + len(pattern) + 8
    # Find Nr of X Scans
    pointsPerScanX = struct.unpack_from('h', mm[binary_start: binary_start + 4])[0]
    bXvals = binary_start + 4
    my_x_range = bXvals + 8 * pointsPerScanX
    startXscan = struct.unpack_from('d', mm[bXvals: bXvals + 8])[0]
    stopXscan = struct.unpack_from('d', mm[my_x_range - 8:my_x_range])[0]
    stepX = (stopXscan - startXscan) / (pointsPerScanX - 1)
    # Find Nr of Y Scans
    pointsPerScanY = struct.unpack_from('h', mm[my_x_range:my_x_range + 4])[0]
    bYvals = my_x_range + 4
    dataStart = bYvals + 8 * pointsPerScanY
    # leeft is an indicator whether the file was interrupted
    left = mm.size() - dataStart
    bRead = dataStart
    # Initialize total_counts, cycles and scan_Y_Channels counters
    total_counts = 0
    cycles = 1
    scanYCh = 0
    # Create a raw dictionary to populate the acquired data
    raw = []
    # Start reading the Binary part of the file
    while left > 0:
        if not struct.unpack_from('h', mm[bRead:bRead + 4])[0] * pointsPerScanX >= left:
            for scanXCh in range(pointsPerScanX):
                scanYval = struct.unpack_from('d', mm[bYvals + 8 * scanYCh: bYvals + 8 + 8 * scanYCh])[0]
                if left <= 0:
                    break
                else:
                    scanXval = struct.unpack_from('d', mm[bXvals + 8 * scanXCh: bXvals + 8 + 8 * scanXCh])[0]
                    recSize = struct.unpack_from('h', mm[bRead:bRead + 4])[0]
                    bRead += 2
                    time_stamp = struct.unpack_from('I', mm[bRead:bRead + 4])[0]
                    bRead += 4
                    left -= 6
                    nrShort = int ( (recSize - 4) / 2 )
                    if nrShort != 0:
                        if nrShort < mca_ch:
                            Ch_array = []
                            Nr_array = []
                            for i in range(0, nrShort, 2):
                                Ch = struct.unpack_from('h', mm[bRead + 2 * i: bRead + 2 + 2 * i])[0]
                                Nr = struct.unpack_from('h', mm[bRead + 2 * i + 2 : bRead + 2 * i + 4])[0]
                                total_counts += Nr
                                for i in range(Nr):
                                    Ch_array.append(Ch * mca_fact)
                                    Nr_array.append(1)
                            raw.append({'cycles': cycles,
                                        'points_per_scanX': scanXval,
                                        'points_per_scanY': scanYval,
                                        'time_stamp': time_stamp,
                                        'mca_channel': Ch_array,
                                        'counts': Nr_array,
                                        'NrEvents': len(Ch_array)})
                    bRead += 2 * nrShort
                    left -= 2 * nrShort
            if scanYCh == pointsPerScanY - 1:
                scanYCh = 0
            else:
                scanYCh += 1
            cycles += 1
            if left <= 0:
                break
        else:
            break
    # close the map
    mm.close()
    
import numpy as np
import matplotlib.pyplot as plt
plt.close('all')
plt.rcParams.update({'font.size': 8})
plt.rcParams.update({'pdf.fonttype': 42})

cm = 1 / 2.54
figsize = (9.3*cm, 5*cm)

x = np.zeros(len(raw))
time = np.zeros(len(raw))
for i, dictio in enumerate(raw):
    x[i] = dictio['points_per_scanX']
    time[i] = dictio['time_stamp']


f_t = np.zeros((total_counts, 3))
cts = 0
for i, dictio in enumerate(raw):
    # z-cut
    if len(dictio['mca_channel']) <= 3:
        for ch in dictio['mca_channel']:
            f_t[cts, 0] = dictio['points_per_scanX']
            f_t[cts, 1] = ch
            f_t[cts, 2] = dictio['cycles']
            cts += 1
            
f_t = f_t[:cts]
f_t = f_t[f_t[:,1] >= 60]
f_t = f_t[f_t[:,1] <= 130]

print(cts, 'counts')

plt.subplots(figsize=figsize)

frf = np.linspace(f_t[:,0].min(), f_t[:,0].max(), pointsPerScanX)
mean_ToF = np.zeros(pointsPerScanX)
ToF_std = np.zeros(pointsPerScanX)
f_cts = np.zeros(pointsPerScanX)
for i, f in enumerate(frf):
    mean_ToF[i] = f_t[:,1][f_t[:,0].round(3)==f.round(3)].mean()
    ToF_std[i] = f_t[:,1][f_t[:,0].round(3)==f.round(3)].std()
    f_cts[i] = len(f_t[:,1][f_t[:,0].round(3)==f.round(3)])

plt.plot(f_t[:,0]-1494216, f_t[:,1], '.', c='0.1', alpha=0.3, mew=0)
plt.errorbar(frf-1494216, mean_ToF, 10/np.sqrt(f_cts), fmt='_', capsize=2, lw=0.8, c='k', ms=2)

plt.ylim(60,130)
plt.xlim(f_t[:,0].min()-1494216, f_t[:,0].max()-1494216)
plt.xlabel('($\\nu_\mathrm{rf}-1\,494\,216\,$Hz) / Hz')
plt.ylabel('ToF / µs')
plt.minorticks_on()

b = np.loadtxt('fit-fit.txt')
plt.plot(b[:,0]-1494216, b[:,1], c='firebrick', lw=1)

plt.savefig('ToF_ICR_file_11_Zn.pdf', bbox_inches='tight', pad_inches=0.01)
