Skip to content
Snippets Groups Projects

Add acoustics monitoring

Merged Add acoustics monitoring
cguidi/km3mon:undefined into master
1 unresolved thread
Merged Carlo Guidi requested to merge cguidi/km3mon:undefined into master
1 unresolved thread
Compare
and
Show latest version
1 file
+ 173
− 183
Compare changes
  • Side-by-side
  • Inline
scripts/acoustics.py
+ 173
− 183
#!/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
-h --help Show this screen.
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import colors
import km3pipe as kp
import argparse
import time
import os
import matplotlib
from matplotlib import colors
from datetime import datetime
import argparse
parser = argparse.ArgumentParser(description='Online_monitoring')
parser.add_argument('-d', dest='det_id', type=str, nargs=1, required=True,
help='The detector ID (e.g. D0ARCA001)')
parser.add_argument('-o', dest='plot_dir', type=str, nargs=1, required=True,
help='Directory in which the plot is saved')
import km3pipe as kp
parser = argparse.ArgumentParser(description = 'Online_monitoring')
parser.add_argument('-d', dest = 'det_id', type = str, nargs = 1, required = True,
help = 'The detector ID (e.g. D_ORCA005)')
parser.add_argument('-o', dest = 'plot_dir', type = str, nargs = 1, required = True,
help = 'Directory in which the plot is saved')
args = parser.parse_args()
detid = args.det_id[0]
@@ -22,137 +38,132 @@ def diff(first, second):
second = set(second)
return [item for item in first if item not in second]
#detid='D_ORCA005'
#print(time.time())
db = kp.db.DBManager()
sds = kp.db.StreamDS()
table=db.run_table(detid)
table = db.run_table(detid)
AB=[12,14,16] # Acoustic Beacons
DOM=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18] # DOMs
DU=[1,2,3,4,5] # DUs
ACOUSTIC_BEACONS = [12, 14, 16]
N_DOMS = 18
DOMS = range(N_DOMS + 1)
N_DUS = 5
DUS = range(1, N_DUS + 1)
#AB=[16] # Acoustic Beacons
#DOM=[7] #
#DU=[3] # DUs
TIT = 600 # Time Interval between Trains of acoustic pulses)
SSW = 160 # Signal Security Window (Window size with signal)
check=True
while check==True:
check = True
while check:
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)<600:
minrun=table["RUN"][len(table["RUN"])-1]-1
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)
COL=[]
for i in DU:
DUx=[]
for j in AB:
for k in DOM:
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:
macaddress = db.doms.via_omkey((i,k), detid).dom_id
toas_all = sds.toashort(detid=detid, minrun=minrun, maxrun=maxrun, domid=macaddress, emitterid=j) # Prendere i dati basandosi sul tempo e non sul run????
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-600))
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]-80
UTB_signal_max=UTB_abdom[QF_max_index]+80
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)
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:22]
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)
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)]
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
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)))]
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
# QF_fifth=[k for k in QF_fourth if (abs(k-max(QF_fourth))<abs(k-noise_threshold))]
Q=[]
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()
# UTB_fourth_l.sort()
D=[]
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)):
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):
for d in sorted(D, reverse = True):
del QF_fourth[d]
QF_fifth=QF_fourth
QF_OK=QF_fifth
# print(QF_OK)
NUM=len(QF_OK)
QF_fifth = QF_fourth
QF_OK = QF_fifth
NUM = len(QF_OK) # Number of pulses
print(NUM)
if (NUM>7):
DUx.append(1.5)
elif (NUM<8 and NUM>3):
DUx.append(0.5)
elif (NUM<4 and NUM>0):
DUx.append(-0.5)
elif (NUM==0):
DUx.append(-1.5)
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:
stop=[]
DUx.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)
COL.append(DUx)
N_Pulses_Indicator.append(N_Pulses_Indicator_DU)
@@ -160,74 +171,75 @@ while check==True:
fig = plt.figure()
ax = fig.add_subplot(111)
dom=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]
l=len(dom)
du1AB1=0.9*np.ones(l)
du1AB2=1*np.ones(l)
du1AB3=1.1*np.ones(l)
du2AB1=1.9*np.ones(l)
du2AB2=2*np.ones(l)
du2AB3=2.1*np.ones(l)
du3AB1=2.9*np.ones(l)
du3AB2=3*np.ones(l)
du3AB3=3.1*np.ones(l)
du4AB1=3.9*np.ones(l)
du4AB2=4*np.ones(l)
du4AB3=4.1*np.ones(l)
du5AB1=4.9*np.ones(l)
du5AB2=5*np.ones(l)
du5AB3=5.1*np.ones(l)
N_doms = 18
doms = range(N_doms + 1)
L = len(doms)
du1AB1 = 0.9*np.ones(L)
du1AB2 = 1*np.ones(L)
du1AB3 = 1.1*np.ones(L)
du2AB1 = 1.9*np.ones(L)
du2AB2 = 2*np.ones(L)
du2AB3 = 2.1*np.ones(L)
du3AB1 = 2.9*np.ones(L)
du3AB2 = 3*np.ones(L)
du3AB3 = 3.1*np.ones(L)
du4AB1 = 3.9*np.ones(L)
du4AB2 = 4*np.ones(L)
du4AB3 = 4.1*np.ones(L)
du5AB1 = 4.9*np.ones(L)
du5AB2 = 5*np.ones(L)
du5AB3 = 5.1*np.ones(L)
DU1=np.array(COL[0])
DU2=np.array(COL[1])
DU3=np.array(COL[2])
DU4=np.array(COL[3])
DU5=np.array(COL[4])
DU1 = np.array(N_Pulses_Indicator[0])
DU2 = np.array(N_Pulses_Indicator[1])
DU3 = np.array(N_Pulses_Indicator[2])
DU4 = np.array(N_Pulses_Indicator[3])
DU5 = np.array(N_Pulses_Indicator[4])
ind=np.where(DU2<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))
ind = np.where(DU2 < 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]
bounds=[-2, -1, 0, 1, 2]
norma = colors.BoundaryNorm(bounds, CustomCmap.N)
color=ax.scatter(du1AB1,dom,s=20,c=DU1[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du1AB2,dom,s=20,c=DU1[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du1AB3,dom,s=20,c=DU1[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du2AB1,dom,s=20,c=DU2[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du2AB2,dom,s=20,c=DU2[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du2AB3,dom,s=20,c=DU2[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du3AB1,dom,s=20,c=DU3[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du3AB2,dom,s=20,c=DU3[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du3AB3,dom,s=20,c=DU3[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du4AB1,dom,s=20,c=DU4[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du4AB2,dom,s=20,c=DU4[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du4AB3,dom,s=20,c=DU4[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du5AB1,dom,s=20,c=DU5[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du5AB2,dom,s=20,c=DU5[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du5AB3,dom,s=20,c=DU5[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du1AB1, doms, s = 20, c = DU1[iAB1], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du1AB2, doms, s = 20, c = DU1[iAB2], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du1AB3, doms, s = 20, c = DU1[iAB3], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du2AB1, doms, s = 20, c = DU2[iAB1], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du2AB2, doms, s = 20, c = DU2[iAB2], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du2AB3, doms, s = 20, c = DU2[iAB3], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du3AB1, doms, s = 20, c = DU3[iAB1], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du3AB2, doms, s = 20, c = DU3[iAB2], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du3AB3, doms, s = 20, c = DU3[iAB3], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du4AB1, doms, s = 20, c = DU4[iAB1], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du4AB2, doms, s = 20, c = DU4[iAB2], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du4AB3, doms, s = 20, c = DU4[iAB3], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du5AB1, doms, s = 20, c = DU5[iAB1], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du5AB2, doms, s = 20, c = DU5[iAB2], norm=norma, marker = 's', cmap = CustomCmap);
color = ax.scatter(du5AB3, doms, s = 20, c = DU5[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')
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, 6, 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)
matplotlib.pyplot.xticks(np.arange(1, 6, 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)
# plt.show()
@@ -240,32 +252,10 @@ while check==True:
print(time.time())
check=False
check_time=time.time()-now
check = False
check_time=time.time() - now
print(check_time)
time.sleep(abs(600-check_time))
check=True
time.sleep(abs(TIT - check_time))
check = True