diff --git a/scripts/acoustics.py b/scripts/acoustics.py
new file mode 100644
index 0000000000000000000000000000000000000000..58d56f6647fe77dffaca920694be34314d7d4c17
--- /dev/null
+++ b/scripts/acoustics.py
@@ -0,0 +1,229 @@
+#!/usr/bin/env python
+# coding=utf-8
+"""
+Online Acoustic Monitoring
+
+Usage:
+ acoustics.py [options]
+ acoustics.py (-h | --help)
+
+Options:
+ -d DET_ID Detector ID
+ -n N_DUS Number of DUs
+ -o PLOT_DIR The directory to save the plot
+ -h --help Show this screen.
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+import os
+import matplotlib
+from matplotlib import colors
+from datetime import datetime
+from docopt import docopt
+import km3pipe as kp
+
+
+def diff(first, second):
+ second = set(second)
+ return [item for item in first if item not in second]
+
+
+args = docopt(__doc__)
+
+detid = args['-d']
+directory = args['-o']
+N_DUS = args['-n']
+N_DUS = int(N_DUS)
+
+db = kp.db.DBManager()
+sds = kp.db.StreamDS()
+
+ACOUSTIC_BEACONS = [12, 14, 16]
+N_DOMS = 18
+N_ABS=3
+DOMS = range(N_DOMS + 1)
+DUS = range(1, N_DUS + 1)
+
+TIT = 600 # Time Interval between Trains of acoustic pulses)
+SSW = 160 # Signal Security Window (Window size with signal)
+
+check = True
+while check:
+
+ table = db.run_table(detid)
+ minrun = table["RUN"][len(table["RUN"]) - 1]
+ ind, = np.where((table["RUN"] == minrun))
+ mintime1 = table['UNIXSTARTTIME'][ind]
+ mintime = mintime1.values
+ maxrun = table["RUN"][len(table["RUN"]) - 1]
+ now = time.time()
+ if (now - mintime/1000) < TIT:
+ minrun = table["RUN"][len(table["RUN"]) - 1] - 1
+ print(now)
+
+ N_Pulses_Indicator = [] # Matrix indicating how many pulses each piezo reveals
+ for du in DUS:
+ N_Pulses_Indicator_DU = [] # Array indicating for each DU how many pulses each piezo reveals.
+ for ab in ACOUSTIC_BEACONS:
+ for dom in DOMS:
+ try:
+ domID = db.doms.via_omkey((du, dom), detid).dom_id
+ toas_all = sds.toashort(detid = detid, minrun = minrun, maxrun = maxrun, domid = domID, emitterid = ab)
+
+ QF_abdom = toas_all["QUALITYFACTOR"]
+ UTB_abdom = toas_all["UNIXTIMEBASE"]
+ TOAS_abdom = toas_all["TOA_S"]
+ UTB_abdom = UTB_abdom.values
+ up = np.where(UTB_abdom>(now-TIT))
+ down = np.where(UTB_abdom<(now))
+ intr = np.intersect1d(up,down)
+ UTB_abdom = UTB_abdom[intr]
+ QF_abdom = QF_abdom[intr]
+ QF_abdom = QF_abdom.values
+ QFlist = QF_abdom.tolist()
+ QFlist.sort(reverse = True)
+ QF_max = max(QF_abdom)
+ QF_max_index = np.where(QF_abdom == QF_max)
+ UTB_signal_min = UTB_abdom[QF_max_index] - SSW/2
+ UTB_signal_max = UTB_abdom[QF_max_index] + SSW/2
+ temp1 = np.where(UTB_abdom > (UTB_signal_min[0]))
+ temp2 = np.where(UTB_abdom < (UTB_signal_max[0]))
+ inter = np.intersect1d(temp1, temp2)
+ inter = inter.tolist()
+ signal_index = inter
+ QF_abdom_index = np.where(QF_abdom)
+ all_data_index = QF_abdom_index[0].tolist()
+ noise_index = diff(all_data_index, signal_index)
+ SIGNAL = QF_abdom[signal_index]
+ UTB_SIGNAL = UTB_abdom[signal_index]
+ NOISE = QF_abdom[noise_index]
+ NOISElist = NOISE.tolist()
+ NOISElist.sort(reverse = True)
+ noise_threshold = max(NOISE)
+
+ # First filter: 22 greatest
+ Security_Number = 22 # To be sure to take all the pulses
+
+ SIGNAL = SIGNAL.tolist()
+ SIGNAL_OLD = np.array(SIGNAL)
+ SIGNAL.sort(reverse = True)
+ QF_first = SIGNAL[0:Security_Number]
+
+ # Second filter: delete duplicates
+
+ QF_second = np.unique(QF_first)
+ QF_second = QF_second.tolist()
+ QF_second.sort(reverse = True)
+
+ # Third filter: If there are more than 11 elements I will eliminate the worst
+
+ if len(QF_second) > 11:
+ QF_second = np.array(QF_second)
+ QF_third = [k for k in QF_second if (np.where(QF_second == k)[0][0] < 11)]
+ else:
+ QF_third = QF_second
+
+ # Fourth filter: I remove the data if it is below the maximum noise
+
+ QF_fourth = [k for k in QF_third if k > (noise_threshold + (5*np.std(NOISE)))]
+
+ # Fifth filter: Check if the clicks are interspersed in the right way
+
+ Q = []
+ for q in np.arange(len(QF_fourth)):
+ Q.append(np.where(SIGNAL_OLD == QF_fourth[q])[0][0])
+ UTB_fourth = np.array(UTB_SIGNAL.tolist())[Q]
+ UTB_fourth_l = UTB_fourth.tolist()
+ D = []
+ for g in np.arange(len(UTB_fourth_l)):
+ if ((np.mod((UTB_fourth_l[g] - UTB_fourth_l[0]), 5) > 0.5 and np.mod((UTB_fourth_l[g] - UTB_fourth_l[0]), 5) < 4.5) or (np.mod((UTB_fourth_l[g] - UTB_fourth_l[0]), 5)>5)):
+ D.append(g)
+ for d in sorted(D, reverse = True):
+ del QF_fourth[d]
+ QF_fifth = QF_fourth
+ QF_OK = QF_fifth
+ NUM = len(QF_OK) # Number of pulses
+ print(NUM)
+
+ if (NUM > 7):
+ N_Pulses_Indicator_DU.append(1.5)
+ elif (NUM < 8 and NUM > 3):
+ N_Pulses_Indicator_DU.append(0.5)
+ elif (NUM < 4 and NUM > 0):
+ N_Pulses_Indicator_DU.append(-0.5)
+ elif (NUM == 0):
+ N_Pulses_Indicator_DU.append(-1.5)
+
+
+ except (TypeError, ValueError): # TypeError if no data found for a certain piezo, ValueError if there are zero data for a certain piezo
+ N_Pulses_Indicator_DU.append(-1.5)
+
+ N_Pulses_Indicator.append(N_Pulses_Indicator_DU)
+
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111)
+ N_doms = 18
+ doms = range(N_doms + 1)
+ L = len(doms)
+
+ duab = []
+ DUs = []
+ for du in range(N_DUS):
+ duabdu = []
+ duab1 = (du+0.9)*np.ones(L)
+ duab2 = (du+1)*np.ones(L)
+ duab3 = (du+1.1)*np.ones(L)
+ duabdu.append(duab1)
+ duabdu.append(duab2)
+ duabdu.append(duab3)
+ duab.append(duabdu)
+ DUs.append(np.array(N_Pulses_Indicator[du]))
+
+ ind = np.where(DUs[1] < 1000)
+ iAB1 = np.where(ind[0] < L)
+ iAB2_up = np.where(ind[0] > (L - 1))
+ iAB2_down = np.where(ind[0] < 2*L)
+ iAB2 = np.intersect1d(iAB2_up, iAB2_down)
+ iAB3 = np.where(ind[0] > (2*L - 1))
+
+ colorsList = [(0, 0, 0), (1, 0.3, 0), (1, 1, 0), (0.2, 0.9, 0)]
+ CustomCmap = matplotlib.colors.ListedColormap(colorsList)
+ bounds = [-2, -1, 0, 1, 2]
+ norma = colors.BoundaryNorm(bounds, CustomCmap.N)
+ for du in range(N_DUS):
+ for ab in range(N_ABS):
+ color = ax.scatter(duab[du][ab], doms, s = 20, c = DUs[du][iAB1], norm=norma, marker = 's', cmap = CustomCmap);
+ color = ax.scatter(duab[du][ab], doms, s = 20, c = DUs[du][iAB2], norm=norma, marker = 's', cmap = CustomCmap);
+ color = ax.scatter(duab[du][ab], doms, s = 20, c = DUs[du][iAB3], norm=norma, marker = 's', cmap = CustomCmap);
+
+ cbar = plt.colorbar(color)
+ cbar.ax.get_yaxis().set_ticks([])
+ for j, lab in enumerate(['$0. pings$', '$1-3 pings$', '$4-7 pings$', '$>7. pings$']):
+ cbar.ax.text(3.5, (2*j + 1)/8.0, lab, ha = 'center', va = 'center')
+ cbar.ax.get_yaxis().labelpad = 18
+
+ matplotlib.pyplot.xticks(np.arange(1, N_DUS + 1, step = 1))
+ matplotlib.pyplot.yticks(np.arange(0, 19, step = 1))
+ matplotlib.pyplot.grid(color = 'k', linestyle = '-', linewidth = 0.2)
+ ax.set_xlabel('DUs', fontsize = 18)
+ ax.set_ylabel('Floors', fontsize = 18)
+ ts = now + 3600
+ DATE = datetime.utcfromtimestamp(ts).strftime('%Y-%m-%d %H:%M:%S')
+ ax.set_title(r' %.16s Detection of the pings emitted by autonomous beacons' %DATE, fontsize = 10)
+
+ my_path = os.path.abspath(directory)
+ my_file = 'Online_Acoustic_Monitoring.png'
+
+ fig.savefig(os.path.join(my_path, my_file))
+
+ print(time.time())
+
+ check = False
+ check_time = time.time() - now
+ print(check_time)
+
+ time.sleep(abs(TIT - check_time))
+ check = True