Compare commits

..

No commits in common. "main" and "minor_fixes" have entirely different histories.

5 changed files with 59 additions and 332 deletions

View file

@ -1,28 +1,13 @@
#! /usr/bin/make -f #! /usr/bin/make -f
MUSS_OPT=-h
MUSC_OPT=-h -c2
%.MUSS.bz2: %.dat.xz %.MUSS.bz2: %.dat.xz
xzcat $< \ xzcat $< \
| ./nm64file -n 32 \ | nm64file -n 32 \
| timeout -v 300 python3 mus_sync.py $(MUSS_OPT) \ | timeout -v 300 python3 mus_sync.py -h \
| bzip2 -9c > $@ | bzip2 -9c > $@
%.MUSS.bz2: %.dat %.MUSS.bz2: %.dat
cat $< \ cat $< \
| ./nm64file -n 32 \ | nm64file -n 32 \
| timeout -v 300 python3 mus_sync.py $(MUSS_OPT) \ | timeout -v 300 python3 mus_sync.py -h \
| bzip2 -9c > $@ | bzip2 -9c > $@
%.MUSC: %.dat.xz
xzcat $< \
| ./nm64file -n 32 \
| timeout -v 300 python3 mus_sync.py $(MUSC_OPT) \
> $@
%.MUSC: %.dat
cat $< \
| ./nm64file -n 32 \
| timeout -v 300 python3 mus_sync.py $(MUSC_OPT) \
> $@

View file

@ -1,82 +0,0 @@
#! /usr/bin/python3
import sys
import numpy as np
from matplotlib import pyplot as plt
verbose=True
f = sys.stdin
l=[]
for line in f:
if 'EOF' in line:
continue
if line[0] == 'E':
l.append(list(map(int,line.split()[5:])))
data=np.array(l)
if verbose:
print("data.shape: ", data.shape, file=sys.stderr)
data=data.reshape(data.shape[0],18,3)
print(data[0])
lmus=[]
lnichtmus=[]
for line in data:
if line[5,1]<2 or line[6,1]<2 or line[9,1]<2 or line[10,1]<2:
continue
else:
if line[0,1]>=2:
lmus.append(line)
else:
lnichtmus.append(line)
lnichtmus = np.array(lnichtmus)
i=0
for line in lnichtmus[:,:,0]:
print (line)
i+=1
if i>200:
break
l=[]
for line in lnichtmus[:,5,0]:
l.append(line)
for line in lnichtmus[:,6,0]:
l.append(line)
for line in lnichtmus[:,9,0]:
l.append(line)
for line in lnichtmus[:,10,0]:
l.append(line)
l=np.array(l)
plt.hist(l/14000,bins=100, range=(0.000000001,73))
plt.xlabel('Pulseheight in mV')
plt.ylabel('Counts')
plt.show()
'''
tab = f.read()
f.close
#print(tab)
tab = tab.split('\n')
h = [x for x in tab if 'E' in x]
#print(h)
l = []
for i in h :
l.append([i])
print(l)
'''

View file

@ -1,21 +0,0 @@
import numpy as np
import getopt
import fileinput
import sys
opt,files = getopt.getopt(sys.argv[1:], 'v')
verbose=False
for o,v in opt:
if o=='-v':
verbose=True
data=[]
for line in fileinput.input(files):
if line[0:2]=='EI':
l=line.split()
if int(l[6])>0:
print(l[5],l[6],l[7])

View file

@ -8,9 +8,6 @@ All lines after the last C64 line in the .MUS file are discarded, if the script
If executed with otion -h this script starts writing at the first C64 line in the hour that startet in the .MUS file. The last If executed with otion -h this script starts writing at the first C64 line in the hour that startet in the .MUS file. The last
line printed is the last EMUS line, before the first C64 line in the hour that started approx one hour after the start of the file. line printed is the last EMUS line, before the first C64 line in the hour that started approx one hour after the start of the file.
Thus there are no redundat lines, and no missing lines, when concatenating multiple outputfiles of this script. Thus there are no redundat lines, and no missing lines, when concatenating multiple outputfiles of this script.
Option -v writes verbosity in sys.stderr
Option -c1 adds Coincidences to the outputfile after the housekeeping blocks
Option -c2 has the ame output as -c1 but all EMUS lines are not in the outputfile
''' '''
import numpy as np import numpy as np
import getopt import getopt
@ -19,22 +16,14 @@ import sys
import datetime import datetime
import time import time
opt,files = getopt.getopt(sys.argv[1:], 'h')
opt,files = getopt.getopt(sys.argv[1:], 'hvc:')
hourly=False hourly=False
verbose=False
coincopt=False
for o,v in opt: for o,v in opt:
if o=='-h': if o=='-h':
hourly=True hourly=True
if o=='-v':
verbose=True
if o=='-c':
coincopt=int(v)
print(f'FMUS {" ".join(sys.argv)}')
def fixclockflip(data): def fixclockflip(data):
clockflipch1=0 clockflipch1=0
@ -46,17 +35,11 @@ def fixclockflip(data):
line[2]+=clockflipch1*4294967296 line[2]+=clockflipch1*4294967296
if line[2]<tb1-1000000: if line[2]<tb1-1000000:
clockflipch1+=1 clockflipch1+=1
line[2]+=4294967296
tb1=line[2]
continue
tb1=line[2] tb1=line[2]
else: else:
line[2]+=clockflipch2*4294967296 line[2]+=clockflipch2*4294967296
if line[2]<tb2-1000000: if line[2]<tb2-1000000:
clockflipch2+=1 clockflipch2+=1
line[2]+=4294967296
tb2=line[2]
continue
tb2=line[2] tb2=line[2]
return data return data
@ -74,40 +57,28 @@ def get_first_t(data,ch):
i+=1 i+=1
return data[i][2],i return data[i][2],i
def get_last_t(data,ch):
i=1
if ch==1:
while True:
if data[-i][0]<16:
break
i+=1
if ch==2:
while True:
if data[-i][0]>15:
break
i+=1
return data[-i][2],i
def count_coinc(data): def count_coinc(data):
lastline=np.zeros(3) n=data.shape[0]-1
delta_times=np.zeros(n)
count=0 count=0
for line in data: for i in range (n):
if -0.5<line[2]-lastline[2]<0.5 and ((lastline[0] <=15 and line[0]>=16) or (line[0] <=15 and lastline[0]>=16)): delta_times[i]=data[i+1][2]-data[i][2]
for t in delta_times:
if t>=-1.5 and t<=1.5:
count+=1 count+=1
lastline=line
return count return count
def approx_drift(data): def approx_drift(data):
offset1=0 offset1=0
offset2=0 offset2=0
lastt=0 lastt=0
for t in data[:5000,2]: for t in data[:2000,2]:
if t-lastt<offset1: if t-lastt<offset1:
offset1=t-lastt offset1=t-lastt
lastt=t lastt=t
offset1=-offset1 offset1=-offset1
lastt=0 lastt=0
for t in data[-5000:,2]: for t in data[-2000:,2]:
if t-lastt<offset2: if t-lastt<offset2:
offset2=t-lastt offset2=t-lastt
lastt=t lastt=t
@ -115,24 +86,18 @@ def approx_drift(data):
slope=(offset2-offset1)/(get_first_t(data,1)[0]-get_first_t(data[-2000:,:],1)[0]) slope=(offset2-offset1)/(get_first_t(data,1)[0]-get_first_t(data[-2000:,:],1)[0])
if abs(slope)>4e-6 or abs(slope)<2e-7: if abs(slope)>4e-6 or abs(slope)<2e-7:
if found_weird_lines: if found_weird_lines:
raise Exception('something went wrong, getting the approximate drift of the two clocks. There were Lines in the Data, that seem strange.'+str(slope)) raise Exception('something went wrong, getting the approximate drift of the two clocks. There were Lines in the Data, that seem strange.')
else: else:
raise Exception('something went wrong, getting the approximate drift of the two clocks.'+str(slope)) raise Exception('something went wrong, getting the approximate drift of the two clocks.')
return slope,offset1 return slope,offset1
def fixtimes(data,slope,offset,firsttime): def fixtimes(data,slope,offset):
fixing_times=True fixing_times=True
highest_count=0 highest_count=0
additional=0 additional=0
a=0 a=0
t1=get_first_t(data,unit)[0] t1=get_first_t(data,unit)[0]
t2=get_first_t(data,runit)[0] t2=get_first_t(data,runit)[0]
if abs(abs(t1-t2)-abs(offset))>1000000:
if found_weird_lines:
raise Exception('Could not find a Correction of the Times. There may be a problem in the data. No hits on one Unit for more than 1s less than 2000 lines after EMUS '+str(data[0][0])+' '+str(data[0][1])+' '+str(data[0][2])+'. There were Lines in the Data, that seem strange.')
else:
raise Exception('Could not find a Correction of the Times. There may be a problem in the data, No hits on one Unit for more than 1s less than 2000 lines after EMUS '+str(data[0][0])+' '+str(data[0][1])+' '+str(data[0][2])+'.')
te=get_last_t(data,runit)[0]
for line in data: for line in data:
if unit==1: if unit==1:
if line[0]>15: if line[0]>15:
@ -148,29 +113,28 @@ def fixtimes(data,slope,offset,firsttime):
modi=a%2 modi=a%2
highest_count=count highest_count=count
additional=0 additional=0
offset=te-get_last_t(data,runit)[0] offset=t2-get_first_t(data,runit)[0]
if unit==1: if unit==1:
offset=-offset offset=-offset
if firsttime and a>2000 and a%2==modi: if count*1.1<highest_count and a%2==modi:
if highest_count<15:
if found_weird_lines:
raise Exception('A suspiciously low amount of coincidences between the two units were detected. The script did not find a better fit. There were Lines in the Data, that seem strange')
else:
raise Exception('A suspiciously low amount of coincidences between the two units were detected. The script did not find a better fit.')
firsttime=False
break break
if (not firsttime) and a>=3 and a%2==modi:
break
a+=1
for line in data: for line in data:
if unit==1: if unit==1:
if line[0]>15: if line[0]>15:
line[2]-=a*(-1)**(a%2) line[2]+=a*(-1)**(a%2)
if unit==2: if unit==2:
if line[0]<16: if line[0]<16:
line[2]+=a*(-1)**(a%2) line[2]-=a*(-1)**(a%2)
additional+=1 additional+=1
a+=1
if a>20000:
if found_weird_lines:
raise Exception('Could not find a Correction of the Times. There may be a problem in the data, the script could not work around. There were Lines in the Data, that seem strange')
else:
raise Exception('Could not find a Correction of the Times. There may be a problem in the data, the script could not work around.')
toomuch=additional/2 toomuch=additional/2
for line in data: for line in data:
if unit==1: if unit==1:
if line[0]>15: if line[0]>15:
@ -179,7 +143,7 @@ def fixtimes(data,slope,offset,firsttime):
if line[0]<16: if line[0]<16:
line[2]-=(-1)**(a%2)*toomuch line[2]-=(-1)**(a%2)*toomuch
return offset,data,firsttime return offset,data
unit=None unit=None
@ -190,9 +154,7 @@ ctime=0
lasttu1=0 lasttu1=0
lasttu2=0 lasttu2=0
found_weird_lines=False found_weird_lines=False
if verbose: #startt=time.time()
st=time.time()
startt=time.time()
for line in fileinput.input(files): for line in fileinput.input(files):
if 'EOF' in line: if 'EOF' in line:
continue continue
@ -200,29 +162,17 @@ for line in fileinput.input(files):
numbers=line.split()[1:4] numbers=line.split()[1:4]
if int(numbers[0])<16: if int(numbers[0])<16:
dt=int(numbers[2])-lasttu1 dt=int(numbers[2])-lasttu1
if 8<int(numbers[1])<11 and(lasttu1==0 or abs(dt)<1000000 or 2**32-1000000<int(abs(dt))<2**32+1000000): if 8<int(numbers[1])<11 and(lasttu1==0 or abs(dt)<1000000 or 2**32-1000000<int(abs(dt))<2**32+1000000):#wahrscheinlichkeit 1s ohne event in der selben unit ca. 1e-300
datalist.append([int(numbers[0]), int(numbers[1]), float(numbers[2])]) datalist.append([int(numbers[0]), int(numbers[1]), float(numbers[2])])
lasttu1=int(numbers[2]) lasttu1=int(numbers[2])
false1=0 else:
elif lasttu1!=0 and not (abs(dt)<1000000 or 2**32-1000000<int(abs(dt))<2**32+1000000):
if verbose :
sys.stderr.write(line)
false1+=1
if false1>10000:
raise Exception('Unit 1 did not write data for more than one second.')
found_weird_lines=True found_weird_lines=True
else: else:
dt=int(numbers[2])-lasttu2 dt=int(numbers[2])-lasttu2
if 8<int(numbers[1])<11 and(lasttu2==0 or abs(dt)<1000000 or 2**32-1000000<abs(dt)<2**32+1000000): if 8<int(numbers[1])<11 and(lasttu2==0 or abs(dt)<1000000 or 2**32-1000000<abs(dt)<2**32+1000000):
datalist.append([int(numbers[0]), int(numbers[1]), float(numbers[2])]) datalist.append([int(numbers[0]), int(numbers[1]), float(numbers[2])])
lasttu2=int(numbers[2]) lasttu2=int(numbers[2])
false2=0 else:
elif lasttu2!=0 and not (abs(dt)<1000000 or 2**32-1000000<abs(dt)<2**32+1000000):
if verbose:
sys.stderr.write(line)
false2+=1
if false2>10000:
raise Exception('Unit 2 did not write data for more than one second.')
found_weird_lines=True found_weird_lines=True
elif any(char.isdigit() for char in line): elif any(char.isdigit() for char in line):
if line[0:3]=='C64': if line[0:3]=='C64':
@ -245,92 +195,47 @@ for line in fileinput.input(files):
otherlines[ctime].append(line) otherlines[ctime].append(line)
data=np.array(datalist) data=np.array(datalist)
#sys.stderr.write('reading data took '+str(time.time()-startt)+'s\n')
keys=list(otherlines.keys()) #startt=time.time()
keys.sort()
if len(keys)<60:
sys.exit()
if verbose:
sys.stderr.write('reading data took '+str(time.time()-startt)+'s\n')
startt=time.time()
data=fixclockflip(data) data=fixclockflip(data)
if verbose: #sys.stderr.write('fixing clockflip took '+str(time.time()-startt)+'s\n')
sys.stderr.write('fixing clockflip took '+str(time.time()-startt)+'s\n') #startt=time.time()
startt=time.time()
slope,offset=approx_drift(data) slope,offset=approx_drift(data)
if verbose: #sys.stderr.write('getting approximate drift took '+str(time.time()-startt)+'s\n')
sys.stderr.write(str((slope,offset))) #startt=time.time()
sys.stderr.write('getting approximate drift took '+str(time.time()-startt)+'s\n')
startt=time.time()
def fat(data,slope,offset,biggertu1): def fat(data,slope,offset,biggertu1):
fixing_times=True fixing_times=True
i=0 i=0
firsttime=True
if biggertu1: if biggertu1:
while fixing_times: while fixing_times:
if i<data.shape[0]-4001: if i<data.shape[0]-2001:
if firsttime: offset,data[i:i+2000]=fixtimes(data[i:i+2000],slope,offset-2)
offset,data[i:i+10000],firsttime=fixtimes(data[i:i+10000],slope,offset,firsttime)
i+=10000
else:
offset,data[i:i+2000],firsttime=fixtimes(data[i:i+2000],slope,offset,firsttime)
i+=2000 i+=2000
else: else:
offset,data[i:],firsttime=fixtimes(data[i:],slope,offset,firsttime) offset,data[i:]=fixtimes(data[i:],slope,offset)
fixing_times=False fixing_times=False
else: else:
offset=-offset offset=-offset
while fixing_times: while fixing_times:
if i<data.shape[0]-4001: if i<data.shape[0]-2001:
if firsttime: offset,data[i:i+2000]=fixtimes(data[i:i+2000],-slope,offset-2)
offset,data[i:i+10000],firsttime=fixtimes(data[i:i+10000],-slope,offset,firsttime)
i+=10000
else:
offset,data[i:i+2000],firsttime=fixtimes(data[i:i+2000],-slope,offset,firsttime)
i+=2000 i+=2000
else: else:
offset,data[i:],firsttime=fixtimes(data[i:],-slope,offset,firsttime) offset,data[i:]=fixtimes(data[i:],-slope,offset-2)
fixing_times=False fixing_times=False
biggertu1=get_first_t(data,1)[0]>get_first_t(data,2)[0] biggertu1=get_first_t(data,1)[0]>get_first_t(data,2)[0]
fat(data,slope,offset,biggertu1) fat(data,slope,offset,biggertu1)
if verbose:
sys.stderr.write('fixing times took '+str(time.time()-startt)+'s\n')
startt=time.time()
data=data[np.argsort(data[:,2])]
if verbose:
sys.stderr.write('sorting data took '+str(time.time()-startt)+'s\n')
startt=time.time()
if coincopt: #sys.stderr.write('fixing times took '+str(time.time()-startt)+'s\n')
k=0 #startt=time.time()
coinces=np.zeros((32,32)) data=data[np.argsort(data[:,2])]
lastline=data[0,:] #sys.stderr.write('sorting data took '+str(time.time()-startt)+'s\n')
lastline[2]-=10 #startt=time.time()
llastline=data[0,:]
llastline[2]-=10 keys=list(otherlines.keys())
lllastline=data[0,:] keys.sort()
lllastline[2]-=10
for line in data:
if k>=len(keys):
break
if line[2]>keys[k]:
for i in [8,9,10,11,12,13,14,15,24,25,26,27,28,29,30,31]:
for j in [0,1,2,3,4,5,6,7,16,17,18,19,20,21,22,23]:
c=coinces[i,j]+coinces[j,i]
otherlines[keys[k-1]].append('CMUS '+str(i)+' '+str(j)+' '+str(int(c))+'\n')
coinces=np.zeros((32,32))
k+=1
if lastline[2]-llastline[2]<2 and llastline[2]-lllastline[2]>2 and line[2]-lastline[2]>2:
coinces[int(llastline[0]),int(lastline[0])]+=1
lllastline=llastline
llastline=lastline
lastline=line
if verbose:
sys.stderr.write('finding coinces took '+str(time.time()-startt)+'s\n')
startt=time.time()
def printhour(data,keys): def printhour(data,keys):
i=0 i=0
@ -352,7 +257,7 @@ def printhour(data,keys):
if writing: if writing:
print(''.join(otherlines[keys[i]])) print(''.join(otherlines[keys[i]]))
i+=1 i+=1
if writing and coincopt!=2: if writing:
print('EMUS',int(line[0]),int(line[1]),f'{line[2]:.2f}') print('EMUS',int(line[0]),int(line[1]),f'{line[2]:.2f}')
return return
@ -366,9 +271,6 @@ else:
if line[2]>keys[i]: if line[2]>keys[i]:
print(''.join(otherlines[keys[i]])) print(''.join(otherlines[keys[i]]))
i+=1 i+=1
if coincopt!=2:
print('EMUS',int(line[0]),int(line[1]),f'{line[2]:.2f}') print('EMUS',int(line[0]),int(line[1]),f'{line[2]:.2f}')
if verbose: #sys.stderr.write('printing fixed data took '+str(time.time()-startt)+'s\n')
sys.stderr.write('printing fixed data took '+str(time.time()-startt)+'s\n') sys.stderr.write(str(files)+'done\n')
dt=time.time()-st
sys.stderr.write(str(files)+' done in '+f'{dt:.2f}'+'s\n')

View file

@ -1,57 +0,0 @@
import numpy as np
import getopt
import fileinput
import sys
import matplotlib.pyplot as plt
opt,files = getopt.getopt(sys.argv[1:], 'b:c:t:a')
bins=None
col=0
allmus=False
for o,v in opt:
if o=='-b':
bins=int(v)
if o=='-c':
col=int(v)
if o=='-a':
allmus=True
data=[]
with open(files[0], "r", encoding="us-ascii", errors="ignore") as file:
for line in file:
if line[0:2]=='EI':
try:
data.append(list(map(int,line.split()[5:])))
except:
continue
data=np.array(data)
muonmus=[]
for line in data:
if (line[16]>2 and line[19]>2) and (line[28]>2 and line[31]>2) and line[15]>200000 and line[18]>245000 and line[27]>185000 and line[30]>208000:
muonmus.append(line[col])
muonmus=np.array(muonmus)
if allmus:
hist,bins,_=plt.hist(data[:,0],bins=bins,histtype='step',density=True)
else:
if col==15 or col==18 or col==27 or col==30:
hist,bins,_=plt.hist(muonmus,bins=bins,range=[0,0.75e6],histtype='step',density=True)
else:
hist,bins,_=plt.hist(muonmus,bins=bins,histtype='step',density=True)
dateiname=files[0]+'col'+str(col)+'.hist'
np.savez(dateiname,hist=hist,bins=bins)
#plt.xlabel('A')
#plt.ylabel('Density')
#plt.show()