diff --git a/scripts/acoustics.py b/scripts/acoustics.py
new file mode 100755
index 0000000000000000000000000000000000000000..661cf324ea2b992bee7a7f1910d6c76d0861ca97
--- /dev/null
+++ b/scripts/acoustics.py
@@ -0,0 +1,393 @@
+#!/usr/bin/env python
+# coding=utf-8
+"""
+Online Acoustic Monitoring
+
+Usage:
+ acoustics.py [options]
+ acoustics.py (-h | --help)
+
+Options:
+ -d DET_ID Detector ID.
+ -o PLOT_DIR The directory to save the plot [default: plots].
+ -h --help Show this screen.
+
+"""
+from datetime import datetime
+import os
+import time
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from matplotlib import colors
+import numpy as np
+import km3pipe as kp
+from docopt import docopt
+
+
+def diff(first, second):
+ second = set(second)
+ return [item for item in first if item not in second]
+
+
+def duplicates(lst, item):
+ return [i for i, x in enumerate(lst) if x == item]
+
+
+args = docopt(__doc__)
+
+detid = args['-d']
+directory = args['-o']
+
+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 = kp.hardware.Detector(det_id=detid).dus
+DUS_cycle = list(np.arange(max(DUS)) + 1)
+
+TIT = 600 # Time Interval between Trains of acoustic pulses)
+SSW = 160 # Signal Security Window (Window size with signal)
+
+clbmap = kp.db.CLBMap(detid)
+
+check = True
+while check:
+
+ minrun = None
+ while minrun is None:
+ try:
+ 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()
+ now = now - 600
+ if (now - mintime / 1000) < TIT:
+ minrun = table["RUN"][len(table["RUN"]) - 1] - 1
+ print(now)
+ except:
+ pass
+
+ N_Pulses_Indicator = [
+ ] # Matrix indicating how many pulses each piezo reveals
+ for du in DUS_cycle:
+ N_Pulses_Indicator_DU = [
+ ] # Array indicating for each DU how many pulses each piezo reveals.
+ for dom in DOMS:
+ UTB_MIN = []
+ QF_MAX = []
+ for ab in ACOUSTIC_BEACONS:
+ try:
+ domID = clbmap.omkeys[(du, dom)].dom_id
+ except (KeyError, AttributeError):
+ N_Pulses_Indicator_DU.append(-1.5)
+ continue
+
+ try:
+ 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[0][0]] - SSW / 2
+ UTB_signal_max = UTB_abdom[QF_max_index[0][0]] + SSW / 2
+ temp1 = np.where(UTB_abdom > (UTB_signal_min))
+ temp2 = np.where(UTB_abdom < (UTB_signal_max))
+ inter = np.intersect1d(temp1, temp2)
+ inter = inter.tolist()
+ signal_index = inter
+
+ # Define the signal index if the the pings are splitted in two parts inside the window
+
+ if UTB_signal_min < (now - TIT):
+ temp1 = np.where(UTB_abdom > (now - TIT))
+ temp2 = np.where(UTB_abdom < (UTB_signal_max))
+ inter1 = np.intersect1d(temp1, temp2)
+ inter1 = inter1.tolist()
+ temp11 = np.where(UTB_abdom < (now))
+ temp22 = np.where(UTB_abdom > (now - SSW / 2))
+ inter2 = np.intersect1d(temp11, temp22)
+ inter2 = inter2.tolist()
+ inter = np.union1d(inter1, inter2)
+ inter = inter.tolist()
+ signal_index = inter
+ signal_index = np.array(signal_index)
+ signal_index = signal_index.astype(int)
+ signal_index = signal_index.tolist()
+
+ if UTB_signal_max > now:
+ temp1 = np.where(UTB_abdom < (now))
+ temp2 = np.where(UTB_abdom > (UTB_signal_min))
+ inter1 = np.intersect1d(temp1, temp2)
+ inter1 = inter1.tolist()
+ temp11 = np.where(UTB_abdom > ((now - TIT)))
+ temp22 = np.where(UTB_abdom < ((now - TIT) + SSW / 2))
+ inter2 = np.intersect1d(temp11, temp22)
+ inter2 = inter2.tolist()
+ inter = np.union1d(inter1, inter2)
+ inter = inter.tolist()
+ signal_index = inter
+ signal_index = np.array(signal_index)
+ signal_index = signal_index.astype(int)
+ signal_index = signal_index.tolist()
+
+ 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]
+ TOA_SIGNAL = TOAS_abdom[signal_index]
+ NOISE = QF_abdom[noise_index]
+ NOISElist = NOISE.tolist()
+ NOISElist.sort(reverse=True)
+ NOISE = NOISE.tolist()
+ NOISE.sort(reverse=True)
+ noise_threshold = max(
+ NOISE) # To be sure not to take signal
+
+ # 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 (Delete if Unixtimebase + ToA is the same)
+
+ QF_second = QF_first
+
+ R = []
+ for r in np.arange(len(QF_first)):
+ R.append(np.where(SIGNAL_OLD == QF_first[r])[0][0])
+ UTB_first = np.array(UTB_SIGNAL.tolist())[R]
+ TOA_first = np.array(TOA_SIGNAL.tolist())[R]
+
+ UNIX_TOA = UTB_first + TOA_first
+ UNIX_TOA = UNIX_TOA.tolist()
+
+ UNIX_TOA_index = []
+ for x in set(UNIX_TOA):
+ if UNIX_TOA.count(x) > 1:
+ UNIX_TOA_index.append(duplicates(UNIX_TOA, x))
+
+ ind_del = []
+ for i in range(len(UNIX_TOA_index)):
+ ind_del.append(UNIX_TOA_index[i][0])
+
+ for ide in sorted(ind_del, reverse=True):
+ del QF_second[ide]
+
+ QF_second.sort(reverse=True)
+
+ # Old second filter
+
+ # 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 + (10 * np.std(NOISE)))
+ ]
+
+ # Fifth filter: Check if the clicks are interspersed in the right way
+
+ QF_fifth = QF_fourth
+ Q = []
+ for q in np.arange(len(QF_fifth)):
+ Q.append(np.where(SIGNAL_OLD == QF_fifth[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) > 2
+ and np.mod(
+ (UTB_fourth_l[g] - UTB_fourth_l[0]), 5) < 4)
+ 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_fifth[d]
+
+ # Sixth filter:
+
+ QF_sixth = [
+ k for k in QF_fifth
+ if (abs(k - max(QF_fifth)) < abs(k - noise_threshold))
+ ]
+
+ QF_OK = QF_sixth
+
+ P = []
+ for p in np.arange(len(QF_OK)):
+ P.append(np.where(SIGNAL_OLD == QF_OK[p])[0][0])
+ UTB_OK = np.array(UTB_SIGNAL.tolist())[P]
+ UTB_OK_l = UTB_OK.tolist()
+
+ UTB_MIN.append(min(UTB_OK_l))
+
+ max_QF = max(QF_OK)
+
+ QF_MAX.append(max_QF)
+
+ 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)
+
+ # To avoid to take wrong beacon signals
+
+ dim = np.size(QF_MAX)
+ pulse_inter = 5.04872989654541
+
+ for i in range(dim - 1):
+ if (np.mod((UTB_MIN[i] - UTB_MIN[i + 1]), pulse_inter) == 0
+ or np.mod(
+ (UTB_MIN[i] - UTB_MIN[i + 1]), pulse_inter) > 5):
+ if QF_MAX[i] < QF_MAX[i + 1]:
+ N_Pulses_Indicator_DU[3 * dom + i] = -1.5
+ else:
+ N_Pulses_Indicator_DU[3 * dom + i + 1] = -1.5
+ if i == 0 and dim == 3:
+ if (np.mod(
+ (UTB_MIN[i] -
+ UTB_MIN[i + 2]), pulse_inter) == 0 or np.mod(
+ (UTB_MIN[i] - UTB_MIN[i + 2]), pulse_inter) > 5):
+ if QF_MAX[i] < QF_MAX[i + 2]:
+ N_Pulses_Indicator_DU[3 * dom + i] = -1.5
+ else:
+ N_Pulses_Indicator_DU[3 * dom + i + 2] = -1.5
+
+ N_Pulses_Indicator.append(N_Pulses_Indicator_DU)
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111)
+
+ duab = []
+ DUs = []
+ for du in DUS_cycle:
+ duabdu = []
+ duab1 = (du - 0.2) * np.ones(N_DOMS + 1)
+ duab2 = (du) * np.ones(N_DOMS + 1)
+ duab3 = (du + 0.2) * np.ones(N_DOMS + 1)
+ duabdu.append(duab1)
+ duabdu.append(duab2)
+ duabdu.append(duab3)
+ duab.append(duabdu)
+ DUs.append(np.array(N_Pulses_Indicator[du - 1]))
+
+ iAB1 = []
+ iAB2 = []
+ iAB3 = []
+ for i in DOMS:
+ iAB1.append(3 * i)
+ iAB2.append(3 * i + 1)
+ iAB3.append(3 * i + 2)
+
+ 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 DUS:
+ color = ax.scatter(duab[du - 1][0],
+ DOMS,
+ s=20,
+ c=DUs[du - 1][iAB1],
+ norm=norma,
+ marker='s',
+ cmap=CustomCmap)
+ color = ax.scatter(duab[du - 1][1],
+ DOMS,
+ s=20,
+ c=DUs[du - 1][iAB2],
+ norm=norma,
+ marker='s',
+ cmap=CustomCmap)
+ color = ax.scatter(duab[du - 1][2],
+ DOMS,
+ s=20,
+ c=DUs[du - 1][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, max(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(2 * TIT - check_time))
+ check = True
diff --git a/supervisord_template.conf b/supervisord_template.conf
index ea8a7f69229a753676e23ea51cef764df2d097e7..e0bf5bc280626fb4353b5efb63e32d46c38c011b 100644
--- a/supervisord_template.conf
+++ b/supervisord_template.conf
@@ -79,6 +79,11 @@ priority=10
stdout_logfile=logs/%(program_name)s.out.log ; stdout log path, NONE for none; default AUTO
stderr_logfile=logs/%(program_name)s.err.log ; stderr log path, NONE for none; default AUTO
+[program:acoustics]
+command=python -u scripts/acoustics.py -d %(ENV_DETECTOR_ID)s
+stdout_logfile=logs/%(program_name)s.out.log ; stdout log path, NONE for none; default AUTO
+stderr_logfile=logs/%(program_name)s.err.log ; stderr log path, NONE for none; default AUTO
+
[program:ahrs_calibration]
command=python -u scripts/ahrs_calibration.py -d %(ENV_DETECTOR_ID)s -p %(ENV_MONITORING_LIGIER_PORT)s
;process_name=%(program_name)s ; process_name expr (default %(program_name)s)
@@ -172,7 +177,7 @@ programs=weblog,msg_dumper,chatbot
priority=200
[group:monitoring_process]
-programs=ahrs_calibration,dom_activity,dom_rates,pmt_rates,trigger_rates,triggermap,ztplot,rttc
+programs=acoustics,ahrs_calibration,dom_activity,dom_rates,pmt_rates,trigger_rates,triggermap,ztplot,rttc
priority=500
[group:alerts]