Skip to content
Snippets Groups Projects

Add acoustics monitoring

Merged Add acoustics monitoring
cguidi/km3mon:undefined into master
Merged Carlo Guidi requested to merge cguidi/km3mon:undefined into master
Compare
and
Show latest version
1 file
+ 125
− 126
Compare changes
  • Side-by-side
  • Inline
scripts/acoustics.py
+ 125
− 126
@@ -15,13 +15,14 @@ Options:
"""
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
import km3pipe as kp
parser = argparse.ArgumentParser(description='Online_monitoring')
parser.add_argument('-d', dest='det_id', type=str, nargs=1, required=True,
@@ -49,120 +50,118 @@ DOMS = range(N_DOMS + 1)
N_DUS = 5
DUS = range(1, N_DUS + 1)
TIT=600 # Time Interval between Trains of acoustic pulses)
SSW=160 # Signal Security Window (Window size with signal)
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 == True:
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
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
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.
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
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)
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
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]
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
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()
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):
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):
if (NUM > 7):
N_Pulses_Indicator_DU.append(1.5)
elif (NUM<8 and NUM>3):
elif (NUM < 8 and NUM > 3):
N_Pulses_Indicator_DU.append(0.5)
elif (NUM<4 and NUM>0):
elif (NUM < 4 and NUM > 0):
N_Pulses_Indicator_DU.append(-0.5)
elif (NUM==0):
elif (NUM == 0):
N_Pulses_Indicator_DU.append(-1.5)
except:
stop=[]
stop = []
N_Pulses_Indicator_DU.append(-1.5)
N_Pulses_Indicator.append(N_Pulses_Indicator_DU)
@@ -173,64 +172,64 @@ 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)
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)
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])
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, 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);
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))
@@ -239,7 +238,7 @@ while check==True:
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')
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()
@@ -254,11 +253,11 @@ while check==True:
check=False
check_time=time.time()-now
check_time=time.time() - now
print(check_time)
time.sleep(abs(TIT-check_time))
time.sleep(abs(TIT - check_time))
check=True