removed double entries in e2h, created light version of e2h
This commit is contained in:
1556
events2histogram.py
1556
events2histogram.py
File diff suppressed because it is too large
Load Diff
835
life_histogrammer.py
Executable file
835
life_histogrammer.py
Executable file
@@ -0,0 +1,835 @@
|
||||
__version__ = '2024-03-15'
|
||||
|
||||
import os
|
||||
import sys
|
||||
import subprocess
|
||||
import h5py
|
||||
import glob
|
||||
import numpy as np
|
||||
import argparse
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib as mpl
|
||||
import time
|
||||
import signal
|
||||
import logging
|
||||
from datetime import datetime
|
||||
|
||||
#==============================================================================
|
||||
#==============================================================================
|
||||
class Detector:
|
||||
def __init__(self):
|
||||
self.nBlades = 14 # number of active blades in the detector
|
||||
angle = np.deg2rad( 5.1 ) # deg angle of incidence of the beam on the blades (def: 5.1)
|
||||
self.dZ = 4.0 * np.sin(angle) # mm height-distance of neighboring pixels on one blade
|
||||
self.dX = 4.0 * np.cos(angle) # mm depth-distance of neighboring pixels on one blace
|
||||
self.bladeZ = 10.7 # mm distance between detector blades (consistent with nu!)
|
||||
self.zero = 0.5 * self.nBlades * self.bladeZ # mm vertical center of the detector
|
||||
|
||||
#==============================================================================
|
||||
def pixel2quantity():
|
||||
det = Detector()
|
||||
nPixel = 64 * 32 * det.nBlades
|
||||
pixelID = np.arange(nPixel)
|
||||
(bladeNr, bPixel) = np.divmod(pixelID, 64*32)
|
||||
(bZ, bY) = np.divmod(bPixel, 64)
|
||||
z = det.zero - bladeNr * det.bladeZ - bZ * det.dZ
|
||||
x = (31 - bZ) * det.dX
|
||||
bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*det.bladeZ / detectorDistance) )
|
||||
delta = (det.nBlades/2. - bladeNr) * bladeAngle - np.rad2deg( np.arctan(bZ*det.dZ / ( detectorDistance + bZ * det.dX) ) )
|
||||
quantity = np.vstack((bladeNr.T, bZ.T, bY.T, delta.T, x.T)).T
|
||||
|
||||
return quantity
|
||||
|
||||
#==============================================================================
|
||||
def analyse_ev(event_e, tof_e, yMin, yMax, thetaMin, thetaMax):
|
||||
|
||||
data_e = np.zeros((len(event_e), 10), dtype=float)
|
||||
|
||||
# data_e column description:
|
||||
# 0: wall time / s
|
||||
# 1: pixelID
|
||||
# 2: blade number
|
||||
# 3: z on blade
|
||||
# 4: y on blade
|
||||
# 5: delta / deg = angle on detector
|
||||
# 6: path within detector / mm
|
||||
# 7: lambda / angstrom
|
||||
# 8: theta / deg
|
||||
# 9: q_z / angstrom^-1
|
||||
|
||||
data_e[:,0] = tof_e[:]
|
||||
data_e[:,1] = event_e[:]
|
||||
|
||||
# filter 'strange' tof times > 2 tau
|
||||
if True:
|
||||
filter_e = (data_e[:,0] <= 2*tau)
|
||||
#print(event_e[~filter_e])
|
||||
#print(data_e[~filter_e,0])
|
||||
data_e = data_e[filter_e,:]
|
||||
if np.shape(filter_e)[0]-np.shape(data_e)[0] > 0.5 and verbous:
|
||||
logging.warning(f'## strange times: {np.shape(filter_e)[0]-np.shape(data_e)[0]}')
|
||||
|
||||
pixelLookUp = pixel2quantity()
|
||||
data_e[:,2:7] = pixelLookUp[np.int_(data_e[:,1])-1,:]
|
||||
|
||||
#================================
|
||||
|
||||
# filter y range
|
||||
filter_e = (yMin <= data_e[:,4]) & (data_e[:,4] <= yMax)
|
||||
data_e = data_e[filter_e,:]
|
||||
|
||||
# correct tof for beam size effect at chopper
|
||||
data_e[:,0] -= ( data_e[:,5] / 180. ) * tau
|
||||
|
||||
# effective flight path length
|
||||
#data_e[:,6] = chopperDetectorDistance + data_e[:,6]
|
||||
|
||||
# calculate lambda
|
||||
hdm = 6.626176e-34/1.674928e-27 # h / m
|
||||
data_e[:,7] = 1.e13 * data_e[:,0] * hdm / ( chopperDetectorDistance + data_e[:,6] )
|
||||
|
||||
# theta
|
||||
data_e[:,8] = nu - mu + data_e[:,5]
|
||||
|
||||
# gravity compensation
|
||||
data_e[:,8] += np.rad2deg( np.arctan( 3.07e-10 * ( detectorDistance + data_e[:,6]) * data_e[:,7] * data_e[:,7] ) )
|
||||
|
||||
# filter theta range
|
||||
filter_l = (thetaMin <= data_e[:,8]) & (data_e[:,8] <= thetaMax)
|
||||
data_e = data_e[filter_l,:]
|
||||
|
||||
# q_z
|
||||
data_e[:,9] = 4*np.pi * np.sin( np.deg2rad( data_e[:,8] ) ) / data_e[:,7]
|
||||
|
||||
# filter q_z range
|
||||
#filter_e = (qMin < data_e[:,I7]) & (data_e[:,7] < qMax)
|
||||
#data_e = data_e[filter_e,:]
|
||||
|
||||
return data_e
|
||||
|
||||
#==============================================================================
|
||||
class Meta:
|
||||
# AMOR hdf dataset with associated properties from metadata
|
||||
def __init__(self, filename):
|
||||
self.filename = filename
|
||||
|
||||
fh = h5py.File(filename, 'r', swmr=True)
|
||||
|
||||
# for processing
|
||||
|
||||
self.chopperDistance = float(np.take(fh['/entry1/Amor/chopper/pair_separation'], 0)) # mm
|
||||
# the following is the distance from the sample to the detector entry window, not to the center of rotation
|
||||
self.detectorDistance = float(np.take(fh['/entry1/Amor/detector/transformation/distance'], 0)) # mm
|
||||
self.chopperDetectorDistance = self.detectorDistance - float(np.take(fh['entry1/Amor/chopper/distance'], 0)) # mm
|
||||
|
||||
self.lamdaCut = 2.5 # Aa
|
||||
|
||||
startDate = str(fh['/entry1/start_time'][0].decode('utf-8'))
|
||||
self.startDate = datetime.strptime(startDate, '%Y-%m-%d %H:%M:%S')
|
||||
startDate = datetime.timestamp(self.startDate)
|
||||
self.countingTime = float(np.take(fh['/entry1/Amor/detector/data/event_time_zero'], -1))/1e9 - startDate
|
||||
# not exact for low rates
|
||||
|
||||
ka0 = 0.245 # given inclination of the beam after the Selene guide
|
||||
|
||||
year_date = str(datetime.today()).split(' ')[0].replace("-", "/", 1)
|
||||
|
||||
# deside from where to take the control paralemters
|
||||
try:
|
||||
self.mu = float(np.take(fh['/entry1/Amor/master_parameters/mu/value'], 0))
|
||||
self.nu = float(np.take(fh['/entry1/Amor/master_parameters/nu/value'], 0))
|
||||
self.kap = float(np.take(fh['/entry1/Amor/master_parameters/kap/value'], 0))
|
||||
self.kad = float(np.take(fh['/entry1/Amor/master_parameters/kad/value'], 0))
|
||||
self.div = float(np.take(fh['/entry1/Amor/master_parameters/div/value'], 0))
|
||||
chSp = float(np.take(fh['/entry1/Amor/chopper/rotation_speed/value'], 0))
|
||||
self.chPh = float(np.take(fh['/entry1/Amor/chopper/phase/value'], 0))
|
||||
except (KeyError, IndexError):
|
||||
logging.warning(f" using parameters from nicos cache")
|
||||
#cachePath = '/home/amor/nicosdata/amor/cache/'
|
||||
cachePath = '/home/nicos/amorcache/'
|
||||
value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-mu/'+year_date)).split('\t')[-1]
|
||||
self.mu = float(value)
|
||||
value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-nu/'+year_date)).split('\t')[-1]
|
||||
self.nu = float(value)
|
||||
value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-kap/'+year_date)).split('\t')[-1]
|
||||
self.kap = float(value)
|
||||
value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-kad/'+year_date)).split('\t')[-1]
|
||||
self.kad = float(value)
|
||||
value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-div/'+year_date)).split('\t')[-1]
|
||||
self.div = float(value)
|
||||
value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-ch1_speed/'+year_date)).split('\t')[-1]
|
||||
chSp = float(value)
|
||||
self.chPh = np.nan
|
||||
|
||||
if chSp:
|
||||
self.tau = 30. / chSp
|
||||
else:
|
||||
self.tau = 0
|
||||
|
||||
try: # not yet correctly implemented in nexus template
|
||||
spin = str(fh['/entry1/polarizer/spin_flipper/spin'][0].decode('utf-8'))
|
||||
if spin == "b'p'":
|
||||
self.spin = 'p'
|
||||
elif spin == "b'm'":
|
||||
self.spin = 'm'
|
||||
elif spin == "b'up'":
|
||||
self.spin = 'p'
|
||||
elif spin == "b'down'":
|
||||
self.spin = 'm'
|
||||
elif spin == '?':
|
||||
self.spin = '?'
|
||||
else:
|
||||
self.spin = 'n'
|
||||
except (KeyError, IndexError):
|
||||
self.spin = '?'
|
||||
|
||||
fh.close()
|
||||
|
||||
#==============================================================================
|
||||
def resolveNumbers(dataPath, ident):
|
||||
if ident == '0' or '-' in ident[0]:
|
||||
try:
|
||||
nnr = int(ident)
|
||||
except:
|
||||
logging.error("ERROR: '{}' is no valid file identifier!".format(ident))
|
||||
fileNames = glob.glob(dataPath+'/*.hdf')
|
||||
fileNames.sort()
|
||||
fileName = fileNames[nnr-1]
|
||||
fileName = fileName.split('/')[-1]
|
||||
fileNumber = fileName.split('n')[1].split('.')[0].lstrip('0')
|
||||
numberString = fileNumber
|
||||
numberList = [int(fileNumber)]
|
||||
elif 'a' in ident[0]:
|
||||
fileNumber = fileName.split('n')[1].split('.')[0].lstrip('0')
|
||||
numberString = fileNumber
|
||||
numberList = [int(fileNumber)]
|
||||
else:
|
||||
if 'r' in ident: # substitute 'r' (recent) for the actual number
|
||||
fileName = glob.glob(dataPath+'/*.hdf')[-1]
|
||||
fileName = fileName.split('/')[-1]
|
||||
fileNumber = fileName.split('n')[1].split('.')[0].lstrip('0')
|
||||
ident = ident.replace('r', fileNumber, 1)
|
||||
numberString = ident
|
||||
numberList = get_flist(ident)
|
||||
|
||||
return numberString, numberList
|
||||
|
||||
#==============================================================================
|
||||
def fileNameCreator(dataPath, ident):
|
||||
clas = commandLineArgs()
|
||||
# create path/filename
|
||||
ident=str(ident)
|
||||
if 'a' in ident:
|
||||
fileName = ident
|
||||
|
||||
else:
|
||||
try:
|
||||
nnr = int(ident)
|
||||
except:
|
||||
logging.error("ERROR: '{}' is no valid file identifier!".format(ident))
|
||||
|
||||
if nnr <= 0 :
|
||||
fileName = glob.glob(dataPath+'/*.hdf')[nnr-1]
|
||||
fileName = fileName.split('/')[-1]
|
||||
else:
|
||||
fileName = f'amor{clas.year}n{ident:>06s}'
|
||||
#fileName = 'amor2021n{:>06s}'.format(ident)
|
||||
|
||||
fileName = fileName.split('.')[0]
|
||||
fileName = fileName+'.hdf'
|
||||
fileName = dataPath+fileName
|
||||
#filename = '/home/software/kafka-to-nexus/'+filename
|
||||
#print(fileName)
|
||||
|
||||
fileNumber = fileName.split('n')[1].split('.')[0].lstrip('0')
|
||||
|
||||
return fileName, fileNumber
|
||||
|
||||
#==============================================================================
|
||||
def selectTime(timeMin, timeMax, dataPacketTime_p, dataPacket_p, detPixelID_e, tof_e):
|
||||
dataPacketTime_p = np.array(dataPacketTime_p)/1e9
|
||||
dataPacket_p = np.array(dataPacket_p)
|
||||
detPixelID_e = np.array(detPixelID_e)
|
||||
tof_e = np.array(tof_e)
|
||||
startTime = dataPacketTime_p[0]
|
||||
stopTime = dataPacketTime_p[-1]
|
||||
if timeMin > (stopTime-startTime):
|
||||
logging.error('ERROR: time interval [{} : {}] s outside measurement time range [0 : {}] s'.format(timeMin, timeMax, stopTime-startTime))
|
||||
sys.exit()
|
||||
dataPacket_p = dataPacket_p[(dataPacketTime_p-startTime)<timeMax]
|
||||
dataPacketTime_p = dataPacketTime_p[(dataPacketTime_p-startTime)<timeMax]
|
||||
dataPacket_p = dataPacket_p[(dataPacketTime_p-startTime)>=timeMin]
|
||||
dataPacketTime_p = dataPacketTime_p[(dataPacketTime_p-startTime)>=timeMin]
|
||||
detPixelID_e = detPixelID_e[dataPacket_p[0]:dataPacket_p[-1]]
|
||||
tof_e = tof_e[dataPacket_p[0]:dataPacket_p[-1]]
|
||||
return dataPacket_p, dataPacketTime_p, detPixelID_e, tof_e,
|
||||
|
||||
#==============================================================================
|
||||
class PlotSelection:
|
||||
|
||||
# header / meta data
|
||||
|
||||
def header(self, filename, mu, nu, totalCounts, countingTime, spin):
|
||||
number = filename.split('n')[1].split('.')[0].lstrip('0')
|
||||
if spin == 'p':
|
||||
spin = ' <+|'
|
||||
elif spin == 'm':
|
||||
spin = ' <-|'
|
||||
else:
|
||||
spin = ''
|
||||
header = "#{} \u03bc={:>1.2f} \u03bd={:>1.2f} {} {:>10} cts {:>8.1f} s".format(number, mu+5e-3, nu+5e-3, spin, totalCounts, countingTime)
|
||||
return header
|
||||
|
||||
def headline(self, numberString, totalCounts):
|
||||
headLine = "#{} \u03bc={:>1.2f} \u03bd={:>1.2f} {:>12,} cts {:>8.1f} s".format(numberString, mu+5e-3, nu+5e-3, totalCounts, countingTime)
|
||||
return headLine
|
||||
|
||||
# grids
|
||||
|
||||
def lamda_grid(self):
|
||||
dldl = 0.005 # Delta lambda / lambda
|
||||
lamda_grid = lamdaMin*(1+dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+dldl)+1))
|
||||
return lamda_grid
|
||||
|
||||
def theta_grid(self):
|
||||
det = Detector()
|
||||
|
||||
bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*det.bladeZ / detectorDistance) )
|
||||
blade_grid = np.arctan( np.arange(33) * det.dZ / ( detectorDistance + np.arange(33) * det.dX) )
|
||||
blade_grid = np.rad2deg(blade_grid)
|
||||
stepWidth = blade_grid[1] - blade_grid[0]
|
||||
blade_grid = blade_grid - 0.2 * stepWidth
|
||||
|
||||
delta_grid = []
|
||||
for b in np.arange(det.nBlades-1):
|
||||
delta_grid = np.concatenate((delta_grid, blade_grid), axis=None)
|
||||
blade_grid = blade_grid + bladeAngle
|
||||
delta_grid = delta_grid[delta_grid<blade_grid[0]-0.5*stepWidth]
|
||||
delta_grid = np.concatenate((delta_grid, blade_grid), axis=None)
|
||||
|
||||
theta_grid = nu - mu - np.flip(delta_grid) + 0.5*det.nBlades * bladeAngle
|
||||
|
||||
theta_grid = theta_grid[theta_grid>=thetaMin]
|
||||
theta_grid = theta_grid[theta_grid<=thetaMax]
|
||||
|
||||
return theta_grid
|
||||
|
||||
def q_grid(self):
|
||||
dqdq = 0.010 # Delta q_z / q_z
|
||||
q_grid = qMin*(1.+dqdq)**np.arange(int(np.log(qMax/qMin)/np.log(1+dqdq)))
|
||||
return q_grid
|
||||
|
||||
# create PNG with several plots
|
||||
|
||||
def all(self, numberString, arg, data_e):
|
||||
#cmap='gist_earth'
|
||||
cmap = mpl.cm.gnuplot(np.arange(256))
|
||||
cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1])
|
||||
cmap = mpl.colors.ListedColormap(cmap, name='myColorMap', N=cmap.shape[0])
|
||||
y_grid = np.arange(64)
|
||||
I_yt, bins_y, bins_t = np.histogram2d(data_e[:,4], data_e[:,8], bins = (y_grid, self.theta_grid()))
|
||||
I_lt, bins_l, bins_t = np.histogram2d(data_e[:,7], data_e[:,8], bins = (self.lamda_grid(), self.theta_grid()))
|
||||
I_q, bins_q = np.histogram(data_e[:,9], bins = self.q_grid())
|
||||
q_lim = 4*np.pi*np.array([ max( np.sin(self.theta_grid()[0]*np.pi/180.)/self.lamda_grid()[-1] , 1e-4 ),
|
||||
min( np.sin(self.theta_grid()[-1]*np.pi/180.)/self.lamda_grid()[0] , 0.03 )])
|
||||
if arg == 'lin':
|
||||
#vmin = min(np.min(I_lt), np.min(I_yt))
|
||||
vmin = 0
|
||||
vmax = max(5, np.max(I_lt), np.max(I_yt))
|
||||
else:
|
||||
vmin = 0
|
||||
vmax = max(1, np.log(np.max(I_lt)+.1)/np.log(10)*1.05, np.log(np.max(I_yt)+.1)/np.log(10)*1.05)
|
||||
# I(y, theta)
|
||||
fig = plt.figure()
|
||||
axs = fig.add_gridspec(2,3)
|
||||
myt = fig.add_subplot(axs[0,0])
|
||||
myt.set_title('detector area')
|
||||
myt.set_xlabel('$y ~/~ \\mathrm{bins}$')
|
||||
myt.set_ylabel('$\\theta ~/~ \\mathrm{deg}$')
|
||||
if arg == 'lin':
|
||||
myt.pcolormesh(bins_y, bins_t, I_yt.T, cmap=cmap, vmin=vmin, vmax=vmax)
|
||||
else:
|
||||
myt.pcolormesh(bins_y, bins_t, (np.log(I_yt + 5.e-1) / np.log(10.)).T, cmap=cmap, vmin=vmin, vmax=vmax)
|
||||
# I(lambda, theta)
|
||||
mlt = fig.add_subplot(axs[0,1:])
|
||||
mlt.set_title('angle- and energy disperse')
|
||||
mlt.set_xlabel('$\\lambda ~/~ \\mathrm{\\AA}$')
|
||||
mlt.axes.get_yaxis().set_visible(False)
|
||||
if arg == 'lin':
|
||||
cb = mlt.pcolormesh(bins_l, bins_t, I_lt.T, cmap=cmap, vmin=vmin, vmax=vmax)
|
||||
else:
|
||||
cb = mlt.pcolormesh(bins_l, bins_t, (np.log(I_lt + 5.e-1) / np.log(10.)).T, cmap=cmap, vmin=vmin, vmax=vmax)
|
||||
# I(q_z)
|
||||
lqz = fig.add_subplot(axs[1,:])
|
||||
lqz.set_title('$I(q_z)$')
|
||||
lqz.set_ylabel('$\\log_{10}(\\mathrm{cnts})$')
|
||||
lqz.set_xlabel('$q_z~/~\\mathrm{\\AA}^{-1}$')
|
||||
lqz.set_xlim(q_lim)
|
||||
if arg == 'lin':
|
||||
plt.plot(bins_q[:-1], I_q, color='blue', linewidth=0.5)
|
||||
else:
|
||||
err_q = np.sqrt(I_q+1)
|
||||
low_q = np.where(I_q-err_q>0, I_q-err_q, 0.1)
|
||||
plt.fill_between(bins_q[:-1], np.log(low_q)/np.log(10), np.log(I_q+err_q/2)/np.log(10), color='lightgrey')
|
||||
plt.plot(bins_q[:-1], np.log(I_q+5e-1)/np.log(10), color='blue', linewidth=0.5)
|
||||
lw = I_q[ ((q_lim[0] < bins_q[:-1]) & (bins_q[:-1] < q_lim[1])) ].min()
|
||||
plt.ylim(max(-0.1, np.log(lw+.1)/np.log(10)-0.1), )
|
||||
#
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=2.8, c='r')
|
||||
fig.colorbar(cb, ax=mlt)
|
||||
plt.subplots_adjust(hspace=0.6, wspace=0.1)
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
# create PNG with one plot
|
||||
|
||||
def Iyz(self, numberString, arg, data_e):
|
||||
det = Detector()
|
||||
cmap = mpl.cm.gnuplot(np.arange(256))
|
||||
cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1])
|
||||
cmap = mpl.colors.ListedColormap(cmap, name='myColorMap', N=cmap.shape[0])
|
||||
y_grid = np.arange(64)
|
||||
z_grid = np.arange(det.nBlades*32)
|
||||
I_yz, bins_y, bins_z = np.histogram2d(data_e[:,4], (det.nBlades-data_e[:,2])*32-data_e[:,3], bins = (y_grid, z_grid))
|
||||
if arg == 'log':
|
||||
vmin = 0
|
||||
vmax = max(1, np.log(np.max(I_yt)+.1)/np.log(10)*1.05)
|
||||
plt.pcolormesh(bins_y[:],bins_z[:],(np.log(I_yz+6e-1)/np.log(10)).T, cmap=cmap, vmin=vmin, vmax=vmax)
|
||||
else:
|
||||
plt.pcolormesh(bins_y[:],bins_z[:],I_yz.T, cmap=cmap)
|
||||
plt.xlabel('$y ~/~ \\mathrm{bins}$')
|
||||
plt.ylabel('$z ~/~ \\mathrm{bins}$')
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.colorbar()
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
def Ilt(self, numberString, arg, data_e) :
|
||||
cmap = mpl.cm.gnuplot(np.arange(256))
|
||||
cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1])
|
||||
cmap = mpl.colors.ListedColormap(cmap, name='myColorMap', N=cmap.shape[0])
|
||||
I_lt, bins_l, bins_t = np.histogram2d(data_e[:,7], data_e[:,8], bins = (self.lamda_grid(), self.theta_grid()))
|
||||
if arg == 'log':
|
||||
vmax = max(1, np.log(np.max(I_lt)+.1)/np.log(10)*1.05 )
|
||||
plt.pcolormesh(bins_l, bins_t, (np.log(I_lt+I_lt[I_lt>0].min()/2)/np.log(10.)).T, cmap=cmap, vmin=0, vmax=vmax)
|
||||
else :
|
||||
vmax = max(np.max(I_lt), 5)
|
||||
plt.pcolormesh(bins_l, bins_t, I_lt.T, cmap=cmap, vmin=0, vmax=vmax)
|
||||
plt.xlim(0,)
|
||||
if np.min(bins_t) > 0.01 :
|
||||
plt.ylim(bottom=0)
|
||||
else:
|
||||
plt.ylim(bottom=np.min(bins_t))
|
||||
if np.max(bins_t) < -0.01:
|
||||
plt.ylim(top=0)
|
||||
else:
|
||||
plt.ylim(top=np.max(bins_t))
|
||||
plt.xlabel('$\\lambda ~/~ \\mathrm{\\AA}$')
|
||||
plt.ylabel('$\\theta ~/~ \\mathrm{deg}$')
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.colorbar()
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
def Itz(self, numberString, arg, data_e):
|
||||
det = Detector()
|
||||
cmap = mpl.cm.gnuplot(np.arange(256))
|
||||
cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1])
|
||||
cmap = mpl.colors.ListedColormap(cmap, name='myColorMap', N=cmap.shape[0])
|
||||
time_grid = np.arange(0, tau, 0.0005)
|
||||
z_grid = np.arange(det.nBlades*32+1)
|
||||
|
||||
I_tz, bins_t, bins_z = np.histogram2d(data_e[:,0], 32*det.nBlades-data_e[:,2]*32-data_e[:,3], bins = (time_grid, z_grid))
|
||||
if arg == 'log':
|
||||
vmax = max(2., np.log(np.max(I_tz)+.1)/np.log(10)*1.05 )
|
||||
plt.pcolormesh(bins_t, bins_z, (np.log(I_tz+5.e-1)/np.log(10.)).T, cmap=cmap, vmin=0, vmax=vmax)
|
||||
else :
|
||||
vmax = max(np.max(I_tz), 5)
|
||||
plt.pcolormesh(bins_t, bins_z, I_tz.T, cmap=cmap, vmin=0, vmax=vmax)
|
||||
if True:
|
||||
plt.xlim(0,)
|
||||
plt.ylim(0,)
|
||||
plt.xlabel('$t ~/~ \\mathrm{s}$')
|
||||
plt.ylabel('$z$ pixel row')
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.colorbar()
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
def Iq(self, numberString, arg, data_e):
|
||||
I_q, bins_q = np.histogram(data_e[:,9], bins = self.q_grid())
|
||||
err_q = np.sqrt(I_q+1)
|
||||
q_lim = 4*np.pi*np.array([ max( np.sin(self.theta_grid()[0]*np.pi/180.)/self.lamda_grid()[-1] , 1e-4 ),
|
||||
min( np.sin(self.theta_grid()[-1]*np.pi/180.)/self.lamda_grid()[0] , 0.03 )])
|
||||
if arg == 'log':
|
||||
low_q = np.where(I_q-err_q>0, I_q-err_q, 0.1)
|
||||
plt.fill_between(bins_q[:-1], np.log(low_q)/np.log(10), np.log(I_q+err_q/2)/np.log(10), color='lightgrey')
|
||||
plt.plot(bins_q[:-1], np.log(I_q+5e-1)/np.log(10), color='blue', linewidth=0.5)
|
||||
lw = I_q[ ((q_lim[0] < bins_q[:-1]) & (bins_q[:-1] < q_lim[1])) ].min()
|
||||
plt.ylim(max(-0.1, np.log(lw+.1)/np.log(10)-0.1), )
|
||||
else:
|
||||
plt.plot(bins_q[:-1], I_q, color='blue', linewidth=0.5)
|
||||
plt.ylabel('$\\log_{10}(\\mathrm{cnts})$')
|
||||
plt.xlabel('$q_z ~/~ \\mathrm{\\AA}^{-1}$')
|
||||
plt.xlim(q_lim)
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
def Il(self, numberString, arg, data_e):
|
||||
I_l, bins_l = np.histogram(data_e[:,7], bins = self.lamda_grid())
|
||||
if arg == 'lin':
|
||||
plt.plot(bins_l[:-1], I_l)
|
||||
plt.ylabel('$I ~/~ \\mathrm{cnts}$')
|
||||
else:
|
||||
plt.plot(bins_l[:-1], np.log(I_l+5.e-1)/np.log(10.))
|
||||
plt.ylabel('$\\log_{10} I ~/~ \\mathrm{cnts}$')
|
||||
plt.xlabel('$\\lambda ~/~ \\mathrm{\\AA}$')
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
def It(self, numberString, arg, data_e):
|
||||
I_t, bins_t = np.histogram(data_e[:,8], bins = self.theta_grid())
|
||||
plt.plot( I_t, bins_t[:-1])
|
||||
plt.xlabel('$\\mathrm{cnts}$')
|
||||
plt.ylabel('$\\theta ~/~ \\mathrm{deg}$')
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.savefig(output, format='png', dpi=150)
|
||||
|
||||
def tof(self, numberString, arg, data_e):
|
||||
time_grid = np.arange(0, 1.3*tau, 0.0005)
|
||||
I_t, bins_t = np.histogram(data_e[:,0], bins = time_grid)
|
||||
if arg == 'lin':
|
||||
plt.plot(bins_t[:-1]+tau, I_t)
|
||||
plt.plot(bins_t[:-1], I_t)
|
||||
plt.plot(bins_t[:-1]+2*tau, I_t)
|
||||
else:
|
||||
lI_t = np.log(I_t+5.e-1)/np.log(10.)
|
||||
plt.plot(bins_t[:-1]+tau, lI_t)
|
||||
plt.plot(bins_t[:-1], lI_t)
|
||||
plt.plot(bins_t[:-1]+2*tau, lI_t)
|
||||
plt.ylabel('log(counts)')
|
||||
plt.xlabel('time / s')
|
||||
headline = self.headline(numberString, np.shape(data_e)[0])
|
||||
plt.title(headline, loc='left', y=1.0, c='r')
|
||||
plt.savefig(output, format='png')
|
||||
|
||||
#==============================================================================
|
||||
#==============================================================================
|
||||
def endIt(signal, frame):
|
||||
print('\n# e2h life mode stopped\n')
|
||||
sys.exit(0)
|
||||
|
||||
#==============================================================================
|
||||
def get_flist(flist):
|
||||
# resolve short notation of filenumbers into a list of filenumbers
|
||||
# e.g. '3,4-9:2,12-14' -> '3, 4, 6, 8, 12, 13, 14'
|
||||
out_list=np.array([], dtype=int)
|
||||
if ',' in flist:
|
||||
fsublists = flist.split(',')
|
||||
else:
|
||||
fsublists = [flist]
|
||||
for fsublist in fsublists:
|
||||
if '-' in fsublist:
|
||||
if ':' in fsublist:
|
||||
limits, step = fsublist.split(':', 1)
|
||||
else:
|
||||
step = 1
|
||||
limits = fsublist
|
||||
for number in range(int(limits.split('-', 1)[0]),
|
||||
int(limits.split('-', 1)[1])+1, int(step)):
|
||||
out_list = np.append(out_list, number)
|
||||
else:
|
||||
out_list = np.append(out_list, int(fsublist))
|
||||
|
||||
return out_list
|
||||
|
||||
#==============================================================================
|
||||
def process(dataPath, ident, clas):
|
||||
#================================
|
||||
# constants
|
||||
hdm = 6.626176e-34/1.674928e-27 # h / m
|
||||
#================================
|
||||
# instrument specific parameters
|
||||
#================================
|
||||
global lamdaMin, lamdaMax, qMin, qMax, thetaMin, thetaMax
|
||||
# defaults
|
||||
lamdaCut = 2.5 # Aa used to reshuffle tof
|
||||
# data filtering and folding
|
||||
|
||||
#================================
|
||||
if clas.lambdaRange:
|
||||
lamdaMin = clas.lambdaRange[0]
|
||||
lamdaMax = clas.lambdaRange[1]
|
||||
else:
|
||||
lamdaMin = lamdaCut
|
||||
|
||||
if clas.timeIntervalAbs:
|
||||
timeMin = clas.timeIntervalAbs[0]
|
||||
timeMax = clas.timeIntervalAbs[1]
|
||||
elif clas.timeIntervalInc:
|
||||
timeMin = clas.timeIntervalInc[0] * clas.timeIntervalInc[1]
|
||||
timeMax = clas.timeIntervalInc[0] * (clas.timeIntervalInc[1] + 1.)
|
||||
else:
|
||||
timeMin = 0
|
||||
timeMax = 0
|
||||
|
||||
chopperPhase = clas.chopperPhase
|
||||
tofOffset = clas.TOFOffset
|
||||
thetaMin = clas.thetaRange[0]
|
||||
thetaMax = clas.thetaRange[1]
|
||||
yMin = clas.yRange[0]
|
||||
yMax = clas.yRange[1]
|
||||
qMin = clas.qRange[0]
|
||||
qMax = clas.qRange[1]
|
||||
|
||||
#================================
|
||||
# find and open input file
|
||||
global ev
|
||||
|
||||
data_eSum = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
|
||||
sumTime = 0
|
||||
|
||||
numberString, numberList = resolveNumbers(dataPath, ident)
|
||||
|
||||
for number in numberList:
|
||||
number= str(number)
|
||||
|
||||
filename, fileNumber = fileNameCreator(dataPath, number)
|
||||
|
||||
if verbous:
|
||||
logging.info('events2histogram processing file ->\033[1m {} \033[0m<-'.format(fileNumber))
|
||||
|
||||
for i in range(6):
|
||||
ev = h5py.File(filename, 'r', swmr=True)
|
||||
try:
|
||||
ev['/entry1/Amor/detector/data/event_time_zero'][-1]
|
||||
break
|
||||
except (KeyError, IndexError):
|
||||
ev.close()
|
||||
if i < 5:
|
||||
if verbous:
|
||||
print("no data yet, retrying ({}) ".format(10-2*i), end='\r')
|
||||
time.sleep(2)
|
||||
continue
|
||||
else:
|
||||
if verbous:
|
||||
print("# time-out: no longer waiting for data!\a")
|
||||
return
|
||||
|
||||
# get and process data
|
||||
meta = Meta(filename)
|
||||
|
||||
global mu, nu, tau
|
||||
|
||||
if clas.mu < 98.:
|
||||
mu = clas.mu
|
||||
else:
|
||||
mu = meta.mu + clas.muOffset
|
||||
|
||||
if clas.nu < 98.:
|
||||
nu = clas.nu
|
||||
else:
|
||||
nu = meta.nu
|
||||
|
||||
if clas.chopperSpeed:
|
||||
tau = 30./ clas.chopperSpeed
|
||||
else:
|
||||
tau = meta.tau
|
||||
|
||||
try:
|
||||
chPh
|
||||
except NameError:
|
||||
chPh = meta.chPh
|
||||
spin = meta.spin
|
||||
|
||||
global countingTime, detectorDistance, chopperDetectorDistance
|
||||
detectorDistance = meta.detectorDistance
|
||||
chopperDetectorDistance = meta.chopperDetectorDistance
|
||||
countingTime = meta.countingTime
|
||||
|
||||
if verbous:
|
||||
logging.info(" mu = {:>4.2f} deg, nu = {:>4.2f} deg".format(mu, nu))
|
||||
if spin == 'u':
|
||||
logging.info(' spin <+|')
|
||||
elif spin == 'd':
|
||||
logging.info(' spin <-|')
|
||||
|
||||
try: lamdaMax
|
||||
except NameError: lamdaMax = lamdaMin + tau * hdm/chopperDetectorDistance * 1e13
|
||||
|
||||
tofOffset = tau * chopperPhase / 180. # mismatch of chopper pulse and time-zero
|
||||
tofCut = lamdaCut * chopperDetectorDistance / hdm * 1.e-13 # tof of frame start
|
||||
|
||||
tof_e = np.array(ev['/entry1/Amor/detector/data/event_time_offset'][:], dtype=np.uint64)/1.e9 + tofOffset # tof
|
||||
|
||||
detPixelID_e = np.array(ev['/entry1/Amor/detector/data/event_id'][:], dtype=np.uint64) # pixel index
|
||||
#totalCounts = len(detPixelID_e)
|
||||
|
||||
dataPacket_p = np.array(ev['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) # data packet index
|
||||
|
||||
if timeMax>0: # pick only the time interval defined by `-i` or `-I`
|
||||
dataPacketTime_p = np.array(ev['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.uint64) # data packet index
|
||||
dataPacket_p, dataPacketTime_p, detPixelID_e, tof_e = selectTime(timeMin, timeMax, dataPacketTime_p, dataPacket_p, detPixelID_e, tof_e)
|
||||
countingTime = dataPacketTime_p[-1]-dataPacketTime_p[0]+1
|
||||
|
||||
# command line argument not yet defined
|
||||
#filterThreshold = 0
|
||||
#if filterThreshold:
|
||||
# detPixelID_e, tof_e = filterTof(detPixelID_e, tof_e, dataPacket_p, filterThreshold)
|
||||
|
||||
tof_e = np.where(tof_e<tofCut, tof_e+2.*tau, tof_e)
|
||||
tof_e = np.where(tof_e>tau+tofCut, tof_e-tau, tof_e)
|
||||
|
||||
data_e = analyse_ev(detPixelID_e, tof_e, yMin, yMax, thetaMin, thetaMax)
|
||||
|
||||
ev.close()
|
||||
|
||||
data_eSum = np.append(data_eSum, data_e, axis=0)
|
||||
sumTime += countingTime
|
||||
|
||||
if verbous:
|
||||
logging.info(" total counts = {} in {:6.1f} s".format(np.shape(data_e)[0], sumTime))
|
||||
|
||||
if clas.spy:
|
||||
number = filename.split('n')[1].split('.')[0].lstrip('0')
|
||||
logging.info('chopper speed={:>4.0f} rpm, phase={:>5.3f} deg, tau={} s'.format(30/tau, chPh, tau))
|
||||
logging.info('nr={}, spn={}, cnt={}, tme={}'.format(number, spin, np.shape(data_eSum)[0], sumTime))
|
||||
logging.info('mu={:>1.2f}, nu={:>1.2f}, kap={:>1.2f}, kad={:>1.2f}, div={:>1.2f}'.format(mu, nu, kap, kad, div))
|
||||
sys.exit()
|
||||
|
||||
#================================
|
||||
# plotting data
|
||||
plotType = clas.plot[0]
|
||||
try:
|
||||
arg = clas.plot[1]
|
||||
except IndexError:
|
||||
arg = 'def'
|
||||
plott = PlotSelection()
|
||||
#string = f"plott.{plotType} (numberString, '{arg}', data_e)"
|
||||
try:
|
||||
plotFunction = getattr(plott, plotType)
|
||||
plotFunction(numberString, arg, data_e)
|
||||
plt.close()
|
||||
except Exception as e:
|
||||
logging.error(f"ERROR: '{plotType}' is no known output format!")
|
||||
logging.error(f" original error: {e}")
|
||||
|
||||
#==============================================================================
|
||||
def commandLineArgs():
|
||||
msg = "events2histogram reads the eventstream from an hdf raw file and \
|
||||
creates various histogrammed outputs or plots."
|
||||
clas = argparse.ArgumentParser(description = msg)
|
||||
|
||||
clas.add_argument("-c", "--chopperSpeed",
|
||||
type=float,
|
||||
help ="chopper speed in rpm")
|
||||
clas.add_argument("-d", "--dataPath",
|
||||
help ="relative path to directory with .hdf files")
|
||||
clas.add_argument("-f", "--fileIdent",
|
||||
default='0',
|
||||
help ="file number or offset (if negative)")
|
||||
clas.add_argument("-I", "--timeIntervalInc",
|
||||
nargs=2,
|
||||
type=float,
|
||||
help ="time interval length and increment")
|
||||
clas.add_argument("-i", "--timeIntervalAbs",
|
||||
nargs=2,
|
||||
type=float,
|
||||
help ="absolute time interval to be processed")
|
||||
clas.add_argument("-l", "--lambdaRange",
|
||||
nargs=2,
|
||||
type=float,
|
||||
help ="wavelength range to be used")
|
||||
clas.add_argument("-M", "--muOffset",
|
||||
default=0.,
|
||||
type=float,
|
||||
help ="mu offset")
|
||||
clas.add_argument("-m", "--mu",
|
||||
default=99.,
|
||||
type=float,
|
||||
help ="value of mu")
|
||||
clas.add_argument("-n", "--nu",
|
||||
default=99.,
|
||||
type=float,
|
||||
help ="value of nu")
|
||||
clas.add_argument("-P", "--chopperPhase",
|
||||
default=-5.,
|
||||
type=float,
|
||||
help ="chopper phase offset")
|
||||
clas.add_argument("-p", "--plot",
|
||||
default=['all', 'def'],
|
||||
nargs='+',
|
||||
help ="select what to plot or write")
|
||||
clas.add_argument("-q", "--qRange",
|
||||
default=[0.005, 0.30],
|
||||
nargs=2,
|
||||
type=float,
|
||||
help ="q_z range")
|
||||
clas.add_argument("-T", "--TOFOffset",
|
||||
default=0.0,
|
||||
type=float,
|
||||
help ="TOF zero offset")
|
||||
clas.add_argument("-t", "--thetaRange",
|
||||
default=[-12., 12.],
|
||||
nargs=2,
|
||||
type=float,
|
||||
help ="theta range to be used")
|
||||
clas.add_argument("-Y", "--year",
|
||||
default = str(datetime.today()).split('-')[0],
|
||||
help = "year, the measurement was performed")
|
||||
clas.add_argument("-y", "--yRange",
|
||||
default=[0, 63],
|
||||
nargs=2,
|
||||
type=int,
|
||||
help ="detector y range to be used")
|
||||
|
||||
return clas.parse_args()
|
||||
|
||||
#==============================================================================
|
||||
def get_dataPath(clas):
|
||||
if clas.dataPath:
|
||||
dataPath = clas.dataPath + '/'
|
||||
if not os.path.exists(dataPath):
|
||||
sys.exit('# *** the directory "'+dataPath+'" does not exist ***')
|
||||
elif os.path.exists('./raw'):
|
||||
dataPath = "./raw/"
|
||||
elif os.path.exists('../raw'):
|
||||
dataPath = "../raw/"
|
||||
else:
|
||||
sys.exit('# *** please provide the path to the .hdf data files (-d <path>, default is "./raw")')
|
||||
|
||||
return dataPath
|
||||
|
||||
#==============================================================================
|
||||
def get_directDataPath(clas):
|
||||
#dataPath = clas.dataPath + '/'
|
||||
year = str(datetime.today()).split('-')[0]
|
||||
year_date = str(datetime.today()).split(' ')[0].replace("-", "/", 1)
|
||||
pNr = str(subprocess.getoutput('/usr/bin/grep "proposal\t" /home/amor/nicosdata/amor/cache/nicos-exp/'+year_date)[-1]).split('\'')[1]
|
||||
dataPath = '/home/amor/nicosdata/amor/data/'+year+'/'+pNr+'/'
|
||||
if not os.path.exists(dataPath):
|
||||
sys.exit('# *** the directory "'+dataPath+'" does not exist ***')
|
||||
|
||||
return dataPath
|
||||
|
||||
#==============================================================================
|
||||
def main():
|
||||
global verbous, output, dataPath
|
||||
|
||||
clas = commandLineArgs()
|
||||
|
||||
dataPath = get_dataPath(clas)
|
||||
logging.basicConfig(level=logging.INFO, format='# %(message)s')
|
||||
verbous = True
|
||||
output = 'e2h.png'
|
||||
process(dataPath, clas.fileIdent, clas)
|
||||
|
||||
#==============================================================================
|
||||
#==============================================================================
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
|
||||
|
||||
Reference in New Issue
Block a user