#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Script to process a data file from a DIY / low cost spectrometer.

Based on software from the DIY particle detector by Oliver Keller.
https://github.com/ozel/DIY_particle_detector

@author: Jeff Moe
@date: January 2022
"""

import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import matplotlib.gridspec as gridspec
import decimal as D
import msgpack
import argparse

mpl.rcParams['font.size']=12

parser = argparse.ArgumentParser()
parser.add_argument("file_name", help="filename to import, e.g. pulses.pkl")
args = parser.parse_args()

def poly1(x, a, b):
    return a*x + b

df = pd.read_pickle(args.file_name)

THL = -300

lp = []
for i,p in df.iterrows():
    if p.pulse.min() < (THL):
        lp.append(p.pulse)

lp = np.asarray(lp[:])

count=0
loopcnt =0
areas = []
peaks = []
min_alpha_peak = 1243

min_g = -20
max_g = -3300
min_length = 44
max_length = 120
min_skip = min_length

xpulses = []
ypulses = []

gmax=0
for i,y in enumerate(lp[:]):
    dydx=np.gradient(y)
    gy= dydx.min()
    gx= dydx.argmin()

    while True:
        if y[0] > THL and y[0] < np.abs(THL) and y.min() < THL and y[gx] <= np.abs(THL) and gy < min_g and gy > max_g:
            trigx=np.where(dydx < min_g)[0]
            trigy=dydx[dydx < min_g]
            peakx1=trigx[0]
            crossing_x1 = (y[peakx1+min_length:] > y[peakx1]).argmax() if (y[peakx1+min_length:] > y[peakx1]).any() else -1
            if crossing_x1 > 0:
                peakx2=peakx1+crossing_x1+min_length
                diff = np.abs(peakx2 - peakx1)
                peak = y[peakx1:peakx2].min()
                if diff >= min_length and diff <= max_length and peak <= THL :
                     area = np.absolute(y[peakx1:peakx2]).sum()
                     areas.append(area)
                     peak = np.absolute(peak) + y[peakx1]
                     peaks.append(peak)
                     ypulse = np.roll(y[peakx1-100:peakx2+100],50-y[peakx1:peakx2].argmin())[130:]
                     if ypulse.size > min_length:
                         xpulses.append(list(range(len(ypulse))))
                         ypulses.append(ypulse)
                     if gy < gmax:
                         gmax=gy
                     count+=1
                else:
                    peakx2=peakx1+min_skip
            else:
                peakx2 = peakx1 + min_skip
        else:
            peakx1=0
            remaining = len(y)
            if remaining <  min_length:
                break
            else:
                peakx2=min_skip
        y = np.delete(y, slice(0,peakx2))
        if len(y) < min_length:
            break
        else:
            dydx=np.delete(dydx, slice(0,peakx2))
            gy= dydx.min()
            gx= dydx.argmin()
        loopcnt+=1

time_diff = df.iloc[-1,0] - df.iloc[0,0]
time_diff = time_diff.total_seconds() / 60.0
cpm = count/time_diff
print("detected pulses:",count, "in", round(time_diff), "minutes ->", round(cpm,3), "CPM")
peaks=np.asarray(peaks)