import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from random import randint
import scipy.io.wavfile as wav
from scipy.signal import butter, filtfilt, find_peaks_cwt

#Add new code to save the time series in csv. Recognize the lines because they have #for save

"""
wavelet function: ramp(points, slope)

- points: the total points of the 1-d numpy array
- slope: the slope function's linear parts
"""

def ramp(points, slope):
    x = np.arange(0, points) - (points - 1.0) / 2
    T = (max(x) - min(x)) / 10
    total = []
    for value in x:
        if value < - T / 2:
            total.append(0)
        elif -T / 2 <= value < - T / 4:
            total.append(slope * (value + T / 2))
        elif - T / 4 <= value < T / 4:
            total.append(- slope * value)
        elif T / 4 <= value < T / 2:
            total.append(slope * (value - T / 2))
        else:
            total.append(0)
    return np.array(total)

# scales of the wavelet family

scales = np.linspace(1.6, 3, 8)

"""

continuous wavelet transform:
    convolves data series with wavelet series resulting from calling wavelet 
    function with different parameters and gives as output a 2-d numpy array of 
    rows corresponding to wavelet variations containing the convoluted time 
    series

"""

def cwt(data, wavelet, family):
    output = np.zeros((len(family), len(data)))
    for index, values in enumerate(family):
        wavelet_data = wavelet(values * 10, values)[::-1]
        output[index] = np.convolve(data, wavelet_data, mode='same')
    return output

"""
read wav, fish, mc simulation files
"""

# read filtered pulses: so far mc simulation creates a tiny amount of corrupted
# pulses, so a prefiltering is applied to the mc simulated pulses file.
#pulses = pd.read_csv('./Simulations/mc_pulses.csv')
pulses = pd.read_csv('/sps/km3net/users/vtsourap/Acoustics/features_for_NN_paper/pulses_Em2_comb_pr.csv')#, delimiter=',', error_bad_lines=False)
# cutting the extra info and keeping only the pressure time series
pulsesonly = pulses.iloc[:, :146]
En = pulses.iloc[:,-3]
x0 = pulses.iloc[:,-6]
y0 = pulses.iloc[:,-5]
z0 = pulses.iloc[:,-4]

# read fishlist: The main, if not only, discrete background sound source appears
# to be dolphin clicks. A preprocessing is applied to the sound file separating
# it to 5ms chunks and determining which ones contain dolphin clicks.
fish = np.loadtxt("/sps/km3net/users/vtsourap/Acoustics/features_for_NN_paper/fish_10.txt").astype(np.int64)

# read noise file and define sampling frequency
sf, y = wav.read("/sps/km3net/users/vtsourap/Acoustics/data1hydro//201650198.180208055818.wav")
# dc offset computation. Normalization is applied according to number of bits.
dc_offset = np.mean(y) / (2 ** 15)

"""

sound file calibration

"""

hydroCal = - 201.0  # hydrophone calibration dBV re 1uPascal
boardGain = 33.0  # board gain dB full-scale referenced to 1.0 V
sensitivityDb = hydroCal + boardGain
sensitivityLinear = np.power(10, sensitivityDb / 20.0)  # calculate V / uPa

"""

highpass filter design:
    sixth order transfer function of butterworth type (butterworth filter)

"""

# filter sound and cut low frequence, flattening
freq = [5000, 10000, 50000, 60000]
a, b = butter(6, freq[0], btype='highpass', fs=sf)
c, d = butter(4, freq[1], btype='highpass', fs=sf)
e, f = butter(6, [freq[1], freq[2]], btype='bandpass', fs=sf)
g, h = butter(6, [freq[0], freq[3]], btype='bandpass', fs=sf)

"""

separation into chunks

"""

# chunks time duration in seconds. Must be the same as the dolphin click script.
chunks_dur = 0.005
space = int(sf * chunks_dur)

"""

wav chunks iteration

"""

# as many iterations as the double of our simulation file
i_end = 2 * len(pulsesonly)

"""
initialize feature's dataframe for classification through decision tree:

features: -1: max over mean of the summed columns of the corelation coefficients matrix resulting from the cwt
          -2: max over mean of the sum of the squared values of each column
          -3: max over mean of the product of the columns
          -4: max of the biggest spike of the columns product function
          -5: number of spikes for the columns product
          -5: frobenius norm of the matrix
          -6: l-1 norm of the matrix
          -7: l-2 norm of the matrix
          -8: l-oo norm of the matrix
          -9: max over mean of the raw function
          -10: max over mean of the filtered function
          (default filtering: 10kHz)
          -11: max over mean of the filtered by 5kHz
          -12: max over mean of the filtered by 10 - 50 kHz
          -13: max over mean of the filtered by 5 - 60 kHz
          -14: max over mean of the summed squared columns of the correlation coefficients for 10 - 50 kHz
          -15: - // - for the columns product for 10-50 filtering
          -16: max of the spike of the columns product for 10-50 filtering
          -17: number of spikes for the columns product
"""

col_sum, sq_col_sum, col_prod, max_spike, frob, l1, l2, loo, raw, filt, filt_5, filt_10_50, filt_5_60, \
    sq_col_sum_10_50, col_prod_10_50, max_spike_10_50, x_coord, y_coord, z_coord, energy, label = \
    np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), \
    np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), \
    np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end), np.zeros(i_end),  \
    np.zeros(i_end), np.zeros(i_end), []

time_series_signal = [] #for save
time_series_bkg = [] #for save

for i in range(i_end):
     
    if i % round(i_end * 0.01) == 0:
        print(str(round((i / i_end) * 100)) + '%')

    # avoid fishes: if a chunk contains a fish it is replaced by another until
    # it is not in the fish list.
    # in order to avoid chunk overlap an iterator above the number of the total
    # iterations is used.
    # a sound file much bigger than the double amount of simulated pulses would
    # be optimal.

    i_true = i
    while np.isin(i_true, fish):
        i_true = i_true + i_end

    # create chunk
    chunk = y[space * i_true: space * (i_true + 1)]
    chunk = chunk / 32768.0  # normalization of raw file
    chunk = chunk - dc_offset
    yMicroPascals = chunk / sensitivityLinear
    yPascals = yMicroPascals / 1000000

    # background
    if i % 2 == 0:
        yPascalsPulse = yPascals
        # filter background, cut low frequencies
        Filtered = filtfilt(c, d, yPascals)
        Filtered_5 = filtfilt(a, b, yPascals)
        Filtered_10_50 = filtfilt(e, f, yPascals)
        Filtered_5_60 = filtfilt(g, f, yPascals)

        # wavelet transform
        cc = cwt(Filtered, ramp, scales)
        cc = cc[:, 4: -4]  # exclude false peaks at margins

        cc_10_50 = cwt(Filtered_10_50, ramp, scales)
        cc_10_50 = cc_10_50[:, 4: -4]  # exclude false peaks at margins

        x_coord[i] = 0.0
        y_coord[i] = 0.0
        z_coord[i] = 0.0
        energy[i] = 0.0
        # define label
        label.append("noise")

        time_series_bkg.append( np.concatenate((Filtered,[x_coord[i],y_coord[i],y_coord[i], energy[i]])) ) #for save
        print("saved bkg chunk:", i) #for save

        

    else:
        k = i // 2  # index on csv of filtered pulses

        # define pulse time series
        chunkpulse = pulsesonly.iloc[k, :]

        # add pulse in a random position inside the chunk
        yPascalsPulse = np.zeros(len(yPascals))
        r = randint(50, 550)
        for j in range(len(yPascals)):
            if r <= j < r + len(chunkpulse):
                yPascalsPulse[j] = yPascals[j] + chunkpulse[j - r]
            else:
                yPascalsPulse[j] = yPascals[j]

        # filter background + signal
        Filtered = filtfilt(c, d, yPascalsPulse)
        Filtered_5 = filtfilt(a, b, yPascalsPulse)
        Filtered_10_50 = filtfilt(e, f, yPascalsPulse)
        Filtered_5_60 = filtfilt(g, f, yPascalsPulse)

        # wavelet transform
        cc = cwt(Filtered, ramp, scales)
        cc = cc[:, 4: -4]  # exclude false peaks at margins

        cc_10_50 = cwt(Filtered_10_50, ramp, scales)
        cc_10_50 = cc_10_50[:, 4: -4]  # exclude false peaks at margins

        x_coord[i] = x0.iloc[k]
        y_coord[i] = y0.iloc[k]
        z_coord[i] = z0.iloc[k]
        energy[i] = En.iloc[k]
        # define label
        label.append("signal")

        time_series_signal.append( np.concatenate((Filtered,[x_coord[i],y_coord[i],y_coord[i], energy[i]])) ) #for save
        print("saved signal chunk:", i) #for save


    yPascalsPulse_abs = abs(yPascalsPulse)
    Filtered_abs = abs(Filtered)
    Filtered_5_abs = abs(Filtered_5)
    Filtered_10_50_abs = abs(Filtered_10_50)
    Filtered_5_60_abs = abs(Filtered_5_60)

    # sum on each column
    scc = []
    for col in range(len(cc[1, :])):
        scc.append(sum(abs(cc[:, col])))
    scc_max = max(scc)
    for times in range(5):
        scc.remove(max(scc))
    col_sum[i] = scc_max / np.mean(scc)

    # squared sum
    s2cc = []
    for col in range(len(cc[1, :])):
        summed = 0
        for row in range(len(cc[:, 1])):
            summed = summed + cc[row, col]**2
        s2cc.append(summed)
    s2cc_max = max(s2cc)
    for times in range(5):
        s2cc.remove(max(s2cc))
    sq_col_sum[i] = s2cc_max / np.mean(s2cc)

    s2cc_10_50 = []
    for col in range(len(cc_10_50[1, :])):
        summed = 0
        for row in range(len(cc_10_50[:, 1])):
            summed = summed + cc_10_50[row, col]**2
        s2cc_10_50.append(summed)
    s2cc_10_50_max = max(s2cc_10_50)
    for times in range(5):
        s2cc_10_50.remove(max(s2cc_10_50))
    sq_col_sum_10_50[i] = s2cc_10_50_max / np.mean(s2cc_10_50)

    # multiply on each column
    mcc = []
    for col in range(len(cc[1, :])):
        prod = 1
        for row in range(len(cc[:, 1])):
            prod = prod * cc[row, col]
        mcc.append(prod)
    mcc_max = max(mcc)
    for times in range(5):
        mcc.remove(max(mcc))
    col_prod[i] = mcc_max / np.mean(mcc)

    # multiply on each column
    mcc_10_50 = []
    for col in range(len(cc_10_50[1, :])):
        prod = 1
        for row in range(len(cc_10_50[:, 1])):
            prod = prod * cc_10_50[row, col]
        mcc_10_50.append(prod)
    mcc_10_50_max = max(mcc_10_50)
    for times in range(5):
        mcc_10_50.remove(max(mcc_10_50))
    col_prod_10_50[i] = mcc_10_50_max / np.mean(mcc_10_50)

    # mcc peaks
    '''
    peaks = list(find_peaks_cwt(np.array(mcc), 60))
    max_peak = 0
    for peak in peaks:
        max_peak = max_peak + mcc[peak]
    max_spike[i] = max_peak / len(peaks)
    '''
    max_mcc = max(mcc)
    mcc.remove(max_mcc)
    max_spike[i] = max_mcc / max(mcc)

    # mcc peaks
    '''
    peaks = list(find_peaks_cwt(np.array(mcc_10_50), 60))
    max_peak = 0
    for peak in peaks:
        max_peak = max_peak + mcc_10_50[peak]
    max_spike_10_50[i] = max_peak / len(peaks)
    '''
    max_mcc_10_50 = max(mcc_10_50)
    mcc_10_50.remove(max_mcc_10_50)
    max_spike_10_50[i] = max_mcc_10_50 / max(mcc_10_50)

    # frobenius norm
    frob[i] = np.linalg.norm(cc, ord='fro')

    # l1 norm
    l1[i] = np.linalg.norm(cc, ord=1)

    # l2 norm
    l2[i] = np.linalg.norm(cc, ord=2)

    # loo norm
    loo[i] = np.linalg.norm(cc, ord=np.inf)

    # max over mean of raw function
    raw[i] = max(yPascalsPulse_abs) / np.mean(yPascalsPulse_abs)

    # max over mean of filtered function
    filt[i] = max(Filtered_abs) / np.mean(Filtered_abs)
    filt_5[i] = max(Filtered_5_abs) / np.mean(Filtered_5_abs)
    filt_10_50[i] = max(Filtered_10_50_abs) / np.mean(Filtered_10_50_abs)
    filt_5_60[i] = max(Filtered_5_60_abs) / np.mean(Filtered_5_60_abs)

features = pd.DataFrame({"col_sum": col_sum, "sq_col_sum": sq_col_sum, "col_prod": col_prod, "max_spike": max_spike,
                         "frob_norm": frob, "l1_norm": l1, "l2_norm": l2, "loo_norm": loo, "raw_signal": raw,
                         "filt_signal": filt, "filt_5": filt_5, "filt_10_50": filt_10_50, "filt_5_60": filt_5_60,
                         "sq_col_sum_10_50": sq_col_sum_10_50, "col_prod_10_50": col_prod_10_50,
                         "max_spike_10_50": max_spike_10_50, "x":x_coord, "y":y_coord, "z":z_coord, "energy":energy, "label": label}) 

#features.to_csv('features_Em2_comb_pr.csv', index=False)
features.to_csv('features_Em2_pEnDist.csv', index=False)

pulses["p2p_new"] = np.zeros(len(pulses))
for i in range(len(pulses)):
    pulses.loc[i, "p2p_new"] = max(pulses.iloc[i, :147]) - min(pulses.iloc[i, :147])


save_signal = pd.DataFrame(time_series_signal) #for save
save_bkg = pd.DataFrame(time_series_bkg) #for save
save_signal.to_csv('pulseInNoise.csv', index=False, header=False) #for save
save_bkg.to_csv('bkgOnly.csv', index=False, header=False) #for save

'''
"plots"
i = 4
i_true = 10
k = 10
chunk = y[space * i_true: space * (i_true + 1)]
chunk = chunk / 32768.0  # normalization of raw file
chunk = chunk - dc_offset
yMicroPascals = chunk / sensitivityLinear
yPascals = yMicroPascals / 1000000
chunkpulse = pulsesonly.iloc[k, :]

# add pulse in a random position inside the chunk
yPascalsPulse = np.zeros(len(yPascals))
r = randint(50, 550)
for j in range(len(yPascals)):
    if r <= j < r + len(chunkpulse):
        yPascalsPulse[j] = yPascals[j] + chunkpulse[j - r]
    else:
        yPascalsPulse[j] = yPascals[j]

Filtered = filtfilt(c, d, yPascals)
Filtered_pulse = filtfilt(c, d, yPascalsPulse)
cc = cwt(Filtered, ramp, scales)
cc = cc[:, 4: -4]  # exclude false peaks at margins
cc_pulses = cwt(Filtered_pulse, ramp, scales)
cc_pulses = cc_pulses[:, 4: -4]
mcc = []
for col in range(len(cc[1, :])):
    prod = 1
    for row in range(len(cc[:, 1])):
        prod = prod * cc[row, col]
    mcc.append(prod)

scc = []
for col in range(len(cc[1, :])):
    scc.append(sum(abs(cc[:, col])))

scc_pulses = []
for col in range(len(cc[1, :])):
    scc_pulses.append(sum(abs(cc_pulses[:, col])))

s2cc = []
for col in range(len(cc[1, :])):
    summed = 0
    for row in range(len(cc[:, 1])):
        summed = summed + cc[row, col]**2
    s2cc.append(summed)

s2cc_pulses = []
for col in range(len(cc_pulses[1, :])):
    summed = 0
    for row in range(len(cc_pulses[:, 1])):
        summed = summed + cc_pulses[row, col]**2
    s2cc_pulses.append(summed)


mcc_pulses = []
for col in range(len(cc_pulses[1, :])):
    prod = 1
    for row in range(len(cc_pulses[:, 1])):
        prod = prod * cc_pulses[row, col]
    mcc_pulses.append(prod)

"""
fish
"""
chunk_fish = y[space * i: space * (i + 1)]
chunk_fish = chunk_fish / 32768.0  # normalization of raw file
chunk_fish = chunk_fish - dc_offset
yMicroPascals = chunk_fish / sensitivityLinear
yPascals_fish = yMicroPascals / 1000000
plt.plot(yPascals_fish)
plt.xlabel("time samples [-]")
plt.ylabel("pressure [Pa]")
plt.grid()
plt.show()

"""
fish vs filtered
"""
b, a = butter(7, [22000, 45000], btype='bandstop', fs=sf)
plt.figure(1)
yPascalsFiltered_fish = filtfilt(b, a, yPascals_fish)
plt.plot(yPascals_fish, color="tab:blue", linewidth=0.8)
plt.plot(yPascalsFiltered_fish, color="tab:orange", linewidth=0.8)
plt.xlabel("time samples [-]", fontsize=17)
plt.ylabel("pressure [Pa]", fontsize=17)
plt.grid()
plt.legend(["raw", "filtered"], fontsize=20)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

"""
fish vs filtered frequency
"""
from numpy.fft import rfft, rfftfreq
f = rfftfreq(n=len(yPascals_fish), d=1/144000)
yPascals_fish_freq = abs(rfft(yPascals_fish))
yPascalsFiltered_fish_freq = abs(rfft(yPascalsFiltered_fish))
dB = 20*np.log10(yPascals_fish_freq/(1e-6))
dB_filtered = 20*np.log10(yPascalsFiltered_fish_freq/(1e-6))
def runningMean(x, N):
    y = np.zeros((len(x),))
    for ctr in range(len(x)):
         y[ctr] = np.sum(x[ctr:(ctr+N)])
    return y/N
plt.figure(2)
plt.plot(f[:-20], runningMean(dB, 3)[:-20], linewidth=1.9, color="tab:blue")
plt.plot(f[:-20], runningMean(dB_filtered, 3)[:-20],  linewidth=1.9, color="tab:orange")
plt.xlabel("frequency [Hz]", fontsize=17)
plt.ylabel("sound pressure level [dB]", fontsize=17)
plt.grid()
plt.legend(["raw", "filtered"], fontsize=20)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

"""
pulse raw vs filtered
"""
plt.plot(yPascalsPulse, color='tab:orange', linewidth=1.9)
plt.plot(yPascals, color='tab:blue', linewidth=1.9)
plt.plot(Filtered_pulse, color='tab:green', linewidth=1.9)
plt.xlabel("time samples [-]",  fontsize=17)
plt.ylabel("pressure [Pa]",  fontsize=17)
plt.grid()
plt.legend(["noise + signal", "noise", "filtered"], fontsize=17)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

"""
3 x 2 plot
"""
fig, ax = plt.subplots(2, 2)
ax[0, 0].plot(yPascals, color='tab:blue', linewidth=1.9)
ax[0, 0].plot(Filtered, color='tab:green', linewidth=1.9)
ax[0, 0].grid()
ax[0, 0].set_ylabel("pressure [Pa]",  fontsize=17)
ax[0, 0].set_xlim([0, 700])
ax[0, 0].set_ylim([-0.20, 1.8])
ax[0, 0].legend(["noise", "filtered"], fontsize=10)
ax[0, 0].tick_params(axis='both', which='major', labelsize=14)

ax[1, 0].contourf(abs(cc), vmin=0, vmax=np.max(cc_pulses))
ax[1, 0].set_ylabel("wavelet scales [-]",  fontsize=17)
ax[1, 0].set_xlim([0, 700])
ax[1, 0].tick_params(axis='both', which='major', labelsize=14)
ax[1, 0].set_xlabel("time samples [-]",  fontsize=17)

"""
ax[2, 0].plot(mcc, linewidth=1.9)
ax[2, 0].grid()
ax[2, 0].set_xlim([0, 700])
ax[2, 0].set_ylabel("product of matrix [-]",  fontsize=17)
ax[2, 0].set_ylim([-0.1e-7, 2e-7])
ax[2, 0].set_xlabel("time samples [-]",  fontsize=17)
ax[2, 0].set_xticks(fontsize=14)
ax[2, 0].set_yticks(fontsize=14)
"""

ax[0, 1].plot(yPascalsPulse, color='tab:orange')
ax[0, 1].plot(yPascals, color='tab:blue')
ax[0, 1].plot(Filtered_pulse, color='tab:green')
ax[0, 1].grid()
ax[0, 1].set_xlim([0, 700])
ax[0, 1].set_ylim([-0.20, 1.8])
ax[0, 1].legend(["noise + signal", "noise", "filtered"], fontsize=10)
ax[0, 1].tick_params(axis='both', which='major', labelsize=14)

ax[1, 1].set_xlim([0, 700])
ax[1, 1].contourf(abs(cc_pulses), vmin=0, vmax=np.max(cc_pulses))
ax[1, 1].tick_params(axis='both', which='major', labelsize=14)
ax[1, 1].set_xlabel("time samples [-]",  fontsize=17)

"""
ax[2, 1].plot(mcc_pulses)
ax[2, 1].grid()
ax[2, 1].set_xlim([0, 700])
ax[2, 1].set_ylabel("product of matrix [-]")
ax[2, 1].set_ylim([-0.1e-7, 2e-7])
ax[2, 1].set_xlabel("time samples [-]")
"""

from scipy import signal

fig, ax = plt.subplots(2, 1)
ax[0].plot([i/sf for i in range(len(yPascals))], yPascals)
ax[0].set_ylabel("pressure [Pa]")
ax[0].grid()
ax[0].set_xlim([0, 0.005])
f, t, Sxx = signal.spectrogram(yPascals, sf, nperseg=120, noverlap=80)
t = np.array([0] + list(t) + [0.005])
Sxx = np.concatenate((np.array([[0 for i in range(len(f))]]).T,
                      np.concatenate((Sxx, np.array([[0 for i in range(len(f))]]).T), axis=1)), axis=1)
cmesh = ax[1].pcolormesh(t, f, Sxx, shading='gouraud')
ax[1].set_ylabel('frequency [Hz]')
ax[1].set_xlabel('time [sec]')
ax[1].set_xlim([0, 0.005])
plt.show()

"""
plots with mother wavelet
"""

test = pulsesonly.iloc[5477, :]
time = np.array([i for i in range(int(len(test)/2))]) / sf
time = np.append(-time[::-1], time)

wav1 = ramp(16, 1.6)
time1 = np.array([i for i in range(int(len(wav1)/2))]) / sf
time1 = np.append(-time1[::-1], time1)

wav2 = ramp(18, 1.8)
time2 = np.array([i for i in range(int(len(wav2)/2))]) / sf
time2 = np.append(-time2[::-1], time2)

wav3 = ramp(20, 2.0)
time3 = np.array([i for i in range(int(len(wav3)/2))]) / sf
time3 = np.append(-time3[::-1], time3)

wav4 = ramp(22, 2.2)
time4 = np.array([i for i in range(int(len(wav4)/2))]) / sf
time4 = np.append(-time4[::-1], time4)

wav5 = ramp(24, 2.4)
time5 = np.array([i for i in range(int(len(wav5)/2))]) / sf
time5 = np.append(-time5[::-1], time5)

wav6 = ramp(26, 2.6)
time6 = np.array([i for i in range(int(len(wav6)/2))]) / sf
time6 = np.append(-time6[::-1], time6)

wav7 = ramp(28, 2.8)
time7 = np.array([i for i in range(int(len(wav7)/2))]) / sf
time7 = np.append(-time7[::-1], time7)

wav8 = ramp(30, 3.0)
time8 = np.array([i for i in range(int(len(wav8)/2))]) / sf
time8 = np.append(-time8[::-1], time8)

plt.plot(time1, wav1)
plt.plot(time2, wav2)
plt.plot(time3, wav3)
plt.plot(time4, wav4)
plt.plot(time5, wav5)
plt.plot(time6, wav6)
plt.plot(time7, wav7)
plt.plot(time8, wav8)
plt.grid()

fig, ax = plt.subplots(3, 2)

ax[0, 0].plot(yPascals, color='tab:orange')
ax[0, 0].plot(Filtered, color='tab:green')
ax[0, 0].grid()
ax[0, 0].set_ylabel("pressure [Pa]")
ax[0, 0].set_xlim([0, 700])
ax[0, 0].legend(["noise", "filtered"])

ax[1, 0].contourf(abs(cc), vmin=0, vmax=np.max(cc_pulses))
ax[1, 0].set_ylabel("wavelet scales [-]")
ax[1, 0].set_xlim([0, 700])

ax[2, 0].plot(mcc)
ax[2, 0].set_ylabel("correlation coefficients procuct [-]")
ax[2, 0].set_xlim([0, 700])
ax[2, 0].set_ylim(0.8*min(mcc_pulses), 1.2*max(mcc_pulses))
ax[2, 0].grid()
ax[2, 0].set_xlabel("time_samples [-]")

"""
ax[2, 0].plot(mcc, linewidth=1.9)
ax[2, 0].grid()
ax[2, 0].set_xlim([0, 700])
ax[2, 0].set_ylabel("product of matrix [-]",  fontsize=17)
ax[2, 0].set_ylim([-0.1e-7, 2e-7])
ax[2, 0].set_xlabel("time samples [-]",  fontsize=17)
ax[2, 0].set_xticks(fontsize=14)
ax[2, 0].set_yticks(fontsize=14)
"""

ax[0, 1].plot(yPascalsPulse, color='tab:blue')
ax[0, 1].plot(yPascals, color='tab:orange')
ax[0, 1].plot(Filtered_pulse, color='tab:green')
ax[0, 1].grid()
ax[0, 1].set_xlim([0, 700])
ax[0, 1].legend(["noise + signal", "noise", "filtered"])

ax[1, 1].set_xlim([0, 700])
ax[1, 1].contourf(abs(cc_pulses), vmin=0, vmax=np.max(cc_pulses))

ax[2, 1].plot(mcc_pulses)
ax[2, 1].set_xlim([0, 700])
ax[2, 1].set_ylim(0.8*min(mcc_pulses), 1.2*max(mcc_pulses))
ax[2, 1].grid()
ax[2, 1].set_xlabel("time_samples [-]")

fig, ax = plt.subplots(3, 2)

ax[0, 0].plot(scc)
ax[0, 0].grid()
ax[0, 0].set_ylabel("sum [-]")
ax[0, 0].set_xlim([0, 700])
ax[0, 0].set_ylim([min(scc_pulses), max(scc_pulses)])

ax[0, 1].plot(scc_pulses)
ax[0, 1].plot([0, 700], [np.mean(scc), np.mean(scc)], linestyle='dashdot')
ax[0, 1].grid()
ax[0, 1].set_xlim([0, 700])
ax[0, 1].set_ylim([min(scc_pulses), max(scc_pulses)])

ax[1, 0].plot(s2cc)
ax[1, 0].grid()
ax[1, 0].set_ylabel("squared sum [-]")
ax[1, 0].set_xlim([0, 700])
ax[1, 0].set_ylim([min(s2cc_pulses), max(s2cc_pulses)])

ax[1, 1].plot(s2cc_pulses)
ax[1, 1].grid()
ax[1, 1].set_xlim([0, 700])
ax[1, 1].set_ylim([min(s2cc_pulses), max(s2cc_pulses)])
ax[1, 1].plot([0, 700], [np.mean(s2cc), np.mean(s2cc)], linestyle='dashdot')

ax[2, 0].plot(mcc)
ax[2, 0].grid()
ax[2, 0].set_ylabel("product [-]")
ax[2, 0].set_xlabel("time samples [-]")
ax[2, 0].set_xlim([0, 700])
ax[2, 0].set_ylim([min(mcc_pulses), max(mcc_pulses)])

ax[2, 1].plot(mcc_pulses)
ax[2, 1].grid()
ax[2, 1].set_xlabel("time samples [-]")
ax[2, 1].set_xlim([0, 700])
ax[2, 1].set_ylim([min(mcc_pulses), max(mcc_pulses)])
ax[2, 1].plot([0, 700], [np.mean(mcc), np.mean(mcc)], linestyle='dashdot')

fig, ax1 = plt.subplots()
color = 'tab:blue'
ax1.set_xlabel('peak to peak width [Pa]')
ax1.set_ylabel('count [-]', color=color)
ax1.hist(pulses["width_resampled"], bins=np.linspace(min(pulses["width_resampled"]),max(pulses["width_resampled"]),500))
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:red'
ax2.set_ylabel('cumulative count [-]', color=color)  # we already handled the x-label with ax1
ax2.plot(bins[:-1], hist1_cum, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.grid(which="both")
ax1.set_xticks(np.linspace(0, 2, 21), minor=True)
ax2.grid(which="minor")
ax1.grid(which='both', axis='x', linestyle='--')

fig, ax1 = plt.subplots()

color = 'tab:blue'
ax1.set_xlabel('shower energy ' + r'$\left[\log{E\left[GeV\right]}\right]$' )
ax1.set_ylabel('count [-]', color=color)
ax1.hist(pulses["energy"], bins=np.linspace(10, 11, 50))
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:red'
ax2.set_ylabel('cumulative count [%]', color=color)  # we already handled the x-label with ax1
ax2.plot(hist1[1][:-1], hist1_cum, color=color)
ax2.tick_params(axis='y', labelcolor=color)
ax1.set_xticks(np.linspace(10, 11, 21), minor=True)
ax2.grid(which="minor")
ax1.grid(which='both', axis='x', linestyle='--')

fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.grid(which="both")

bins = np.linspace(0, 200, 1000)
color = 'tab:blue'
fig, ax = plt.subplots()
signal = features[features["label"] == "signal"]["sq_col_sum"]
noise = features[features["label"] == "noise"]["sq_col_sum"]
has = ax.hist(signal, bins)
no = ax.hist(noise, bins, alpha=0.5)
ax.set_ylabel("count [-]", color=color)
ax.tick_params(axis='y', labelcolor=color)
ax.set_xlabel(r"$k_1=\frac{\max_{j}{\sigma_j}}{\overline{\sigma_j}}~\left[-\right]$", fontsize=13)
ax.grid()
ax.legend(["with neutrino", "without neutrino"])
"""
ax.axvline(signal.mean(), color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = ax.set_ylim()
ax.text(signal.mean()*1.1, max_ylim*0.9, 'Mean: {:.2f}'.format(signal.mean()))
ax.axvline(noise.mean(), color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = ax.set_ylim()
ax.text(noise.mean()*1.1, max_ylim*0.9, 'Mean: {:.2f}'.format(noise.mean()))
"""
eff = []
cont = []
has = np.array(has[0])/sum(np.array(has[0]))
no = np.array(no[0])/sum(np.array(no[0]))
for i in range(len(has)):
    eff.append(sum(has[i:]))
    cont.append(sum(no[i:])/sum(has[i:]))

eff = np.array(eff)/eff[0]
pur = 1-np.array(cont)
eff_noise = []
for i in range(len(no[:])):
    eff_noise.append(sum(no[i:]))
eff_noise = np.array(eff_noise)/eff_noise[0]
ax2 = ax.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:red'
ax2.set_ylabel('efficiency, purity, significance [-]', color=color)  # we already handled the x-label with ax1
ax2.tick_params(axis='y', labelcolor=color)
ax2.plot(bins[:-1], eff)
ax2.plot(bins[:-1], pur)
sig = eff/np.sqrt(eff+eff_noise)
ax2.plot(bins[:-1], sig)
ax2.legend(["efficiency", "purity", "significance"])
ax2.axvline(bins[79], color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = ax2.set_ylim()
ax2.text(bins[79]*1.1, max_ylim*0.75, r'$k_1^*$'+'= {:.2f}'.format(bins[79]))
plt.show()



bins = np.linspace(0, 2000, 1000)
signal = features[features["label"] == "signal"]["sq_col_sum"]
noise = features[features["label"] == "noise"]["sq_col_sum"]
plt.hist(signal, bins)
plt.hist(noise, bins, alpha=0.5)
plt.ylabel("count [-]")
plt.xlabel(r"$\frac{\max_{j}{\sigma_{2,j}}}{\overline{\sigma_{2,j}}}~\left[-\right]$", fontsize=17)
plt.grid()
plt.legend(["with neutrino", "without neutrino"])
plt.axvline(signal.mean(), color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = plt.ylim()
plt.text(signal.mean()*1.1, max_ylim*0.9, 'Mean: {:.2f}'.format(signal.mean()))
plt.axvline(noise.mean(), color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = plt.ylim()
plt.text(noise.mean()*1.1, max_ylim*0.9, 'Mean: {:.2f}'.format(noise.mean()))

bins = np.linspace(0, max(signal), 100000)
signal = features[features["label"] == "signal"]["col_prod"]
noise = features[features["label"] == "noise"]["col_prod"]
plt.hist(signal, bins)
plt.hist(noise, bins, alpha=0.5)
plt.ylabel("count [-]")
plt.xlabel(r"$\frac{\max_{j}{\pi_{j}}}{\overline{\pi_{j}}}~\left[-\right]$", fontsize=17)
plt.grid()
plt.legend(["with neutrino", "without neutrino"])
plt.axvline(10000, color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = plt.ylim()
plt.text(10000*0.82, max_ylim*0.9, 'Mean: {:.2f}'.format(signal.mean()))
plt.axvline(noise.mean(), color='k', linestyle='dashed', linewidth=1)
min_ylim, max_ylim = plt.ylim()
plt.text(noise.mean()*1.1, max_ylim*0.9, 'Mean: {:.2f}'.format(noise.mean()))
'''