diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..237fb3f --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +*.hdf filter=lfs diff=lfs merge=lfs -text diff --git a/.github/workflows/unit_tests.yml b/.github/workflows/unit_tests.yml index f46b7db..be278bc 100644 --- a/.github/workflows/unit_tests.yml +++ b/.github/workflows/unit_tests.yml @@ -38,5 +38,4 @@ jobs: - name: Test with pytest run: | - cd tests - python -m pytest --pyargs . + python -m pytest tests diff --git a/.gitignore b/.gitignore index 7f4f5b5..d7fd4ed 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ raw test_results build dist +profile_test.prof +amor_eos.log diff --git a/e2h_new.py b/e2h_new.py deleted file mode 100644 index a33ae52..0000000 --- a/e2h_new.py +++ /dev/null @@ -1,1082 +0,0 @@ -__version__ = '2024-03-15' - -# essential changes with regard to 2022 version -# - imprved orso header -# - nexus compatible -# - new theta grid -# - new content in data_e (angleas rather than distances) -# - bug fixes: tof correction within detector - -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 -from orsopy import fileio - -#============================================================================== -#============================================================================== -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 - # t_cor = (delta / 180 deg) * tau - 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 + np.rad2deg( np.arctan2(data_e[:,5], detectorDistance) ) - 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 - if mu > -0.25: - data_e[:,9] = 4*np.pi * np.sin( np.deg2rad( data_e[:,8] ) ) / data_e[:,7] - else: - 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(self.hdf['/entry1/Amor/instrument_control_parameters/mu'], 0)) - self.nu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/nu'], 0)) - self.kap = 0.245 - #self.kap = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kap/value'], 0)) - self.kad = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kad'], 0)) - self.div = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/div'], 0)) - self.chopperSpeed = float(np.take(self.hdf['/entry1/Amor/chopper/rotation_speed'], 0)) - self.chopperPhase = float(np.take(self.hdf['/entry1/Amor/chopper/phase'], 0)) - self.ch1TriggerPhase = float(np.take(self.hdf['/entry1/Amor/chopper/ch1_trigger_phase'], 0)) - self.ch2TriggerPhase = float(np.take(self.hdf['/entry1/Amor/chopper/ch2_trigger_phase'], 0)) - polarizationConfigLabel = float(np.take(self.hdf['/entry1/Amor/polarization/configuration/value'], 0)) - polarizationConfig = polarizationConfigs[int(polarizationConfigLabel)] - - - 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/' - cachePath = '/home/amor/cache/' - 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 = '?' - - - # for .ort header - - self.title = str(fh['entry1/title'][0].decode('utf-8')) - - self.proposal_id = str(fh['entry1/proposal_id'][0].decode('utf-8')) - - self.userName = str(fh['entry1/user/name'][0].decode('utf-8')) - try: - self.userEmail = str(fh['entry1/user/email'][0].decode('utf-8')) - except: - self.userEmail = None - - self.sampleName = str(fh['entry1/sample/name'][0].decode('utf-8')) - self.sampleModel = str(fh['entry1/sample/model'][0].decode('utf-8')) - - #!!! self.instrumentName = str(fh['entry1/Amor/name'][0]) - self.instrumentName = 'Amor' - #!!! self.instrumentType = str(fh['entry1/Amor/type'][0]) - self.source = str(fh['entry1/Amor/source/name'][0]) - #!!! self.sourceType = str(fh['entry1/Amor/source/type'][0]) - #!!! self.sourceProbe = str(fh['entry1/Amor/source/probe'][0]) - self.sourceProbe = 'neutron' - - 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)=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, - -#============================================================================== -def ort_header(fileName): - - #print(fileName) - meta = Meta(fileName) - ''' - Build information object for ORSO file headers. - ''' - - #for fidx in self.get_flist(self.Files): - # fdata = data_vault['tmp%s'%fidx] - # fdate = datetime.strptime(fdata['timestamp'], '%Y-%m-%d %H:%M:%S') - # datafiles.append(fileio.File(file = fdata['srcname'], timestamp = fdate)) - - inst = fileio.InstrumentSettings( - incident_angle = fileio.ValueRange( - mu + meta.kap + meta.kad - meta.div/2, - mu + meta.kap + meta.kad + meta.div/2, - 'deg' - ), - wavelength = fileio.ValueRange( - lamdaMin, - lamdaMax, - 'angstrom' - ), - ) - inst.mu = fileio.Value(mu, 'deg', comment = 'sample angle to horizon') - inst.nu = fileio.Value(nu, 'deg', comment = 'detector angle to horizon') - inst.kap = fileio.Value(meta.kap, 'deg', comment = 'nominal beam inclination') - inst.kad = fileio.Value(meta.kad, 'deg', comment = 'offset of beam inclination') - inst.div = fileio.Value(meta.div, 'deg', comment = 'incoming beam divergence') - mess = fileio.Measurement( - instrument_settings = inst, - data_files = fileName.split('/')[-1], - counting_time = fileio.Value(meta.countingTime, 's'), - scheme = 'angle- and energy-dispersive') - se_dict = None - smpl = fileio.Sample( - name = meta.sampleName, - model = meta.sampleModel, - sample_parameters = se_dict - ) - experiment = fileio.Experiment( - title = meta.title, - instrument = meta.instrumentName, - #instrument_type = meta.instrumentType, - start_date = meta.startDate, - probe = meta.sourceProbe, - #facility = meta.sourceName, - #source_type = meta.sourceType, - proposalID = meta.proposal_id, - ) - ds = fileio.DataSource( - fileio.Person( - meta.userName, - None, - contact = meta.userEmail - ), - experiment, - smpl, - mess - ) - red = fileio.Reduction( - software = fileio.Software( - 'events2histogram', - version = __version__ - ), - timestamp = datetime.now(), - creator = None, - corrections = None, - computer = None, - call = 'python '+' '.join(sys.argv) - ) - cols = [ - fileio.Column('Qz', '1/angstrom'), - fileio.Column('R', comment = 'uncorrected intensity'), - fileio.Column('sR', comment = 'sigma of gaussian probability function'), - fileio.Column('sQz', '1/angstrom', comment = 'resolution based only on sigma_lambda!') - ] - - header = fileio.Orso(ds, red, cols) - - return header - -#============================================================================== -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 - if foldback: - lamda_grid = lamdaMin*(1+dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+dldl)+1)) - else: - lamda_grid = np.arange(0.01, 2.*lamdaMax-2.*lamdaMin, 0.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=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 )]) - q_lim = np.array([qMin, qMax]) - 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) - #plt.close() - - # 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 == 'file': - print('# y z conts') - for y in range(len(bins_y)-1): - for z in range(len(bins_z)-1): - print(" %6.3f %6.4f %10.3e" %(bins_y[y], bins_z[z], I_yz[y,z])) - print("") - return - elif arg == 'log': - vmin = 0 - vmax = max(1, np.log(np.max(I_yz)+.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) - #plt.close() - - 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 == 'file': - print('# lambda theta conts') - for l in range(len(bins_l)-1): - for t in range(len(bins_t)-1): - print(" %6.3f %6.4f %10.3e" %(bins_l[l], bins_t[t], I_lt[l,t]/ bins_t[t])) - print("") - return - elif 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) - #plt.close() - - 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]) - if foldback: - time_grid = np.arange(0, tau, 0.0005) - else: - time_grid = np.arange(0, 2.*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 == 'file': - print('# time z conts') - for t in range(len(bins_t)-1): - for z in range(len(bins_z)-1): - print(" %6.3f %6.4f %10.3e" %(bins_t[t], bins_z[z], I_tz[t,z])) - print("") - return - elif 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) - #plt.close() - - 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 )]) - q_lim = [qMin, qMax] - print(q_lim) - if arg == 'file': - header = '# q counts' - I_q = np.vstack((bins_q[:-1], I_q, err_q)) - np.savetxt(sys.stdout, I_q.T, header=header) - logging.info(' use `-p ort` instead of `-p Iq`!') - return - elif 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) - #plt.close() - - def Il(self, numberString, arg, data_e): - I_l, bins_l = np.histogram(data_e[:,7], bins = self.lamda_grid()) - if arg == 'file': - header = '# lambda counts' - I_l = np.vstack((bins_l[:-1], I_l)) - np.savetxt(sys.stdout, I_l.T, header=header) - return - elif 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) - #plt.close() - - def It(self, numberString, arg, data_e): - I_t, bins_t = np.histogram(data_e[:,8], bins = self.theta_grid()) - if arg == 'file': - header = '# 2theta counts' - I_t = np.vstack((bins_t[:-1], I_t)) - np.savetxt(sys.stdout, I_t.T, header=header) - return - else: - 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) - #plt.close() - - def tof(self, numberString, arg, data_e): - if foldback: - time_grid = np.arange(0, 1.3*tau, 0.0005) - else: - time_grid = np.arange(0, 2.*tau, 0.0005) - I_t, bins_t = np.histogram(data_e[:,0], bins = time_grid) - if arg == 'file': - header = '# time counts' - I_t = np.vstack((bins_t[:-1], I_t)) - np.savetxt(sys.stdout, I_t.T, header=header) - return - elif 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 ort(self, numberString, arg, data_e): - I_q, bins_q = np.histogram(data_e[:,9], bins = self.q_grid()) - sI_q = np.sqrt(I_q) - sq_q = bins_q[:-1]*0.022 - I_q = np.vstack((bins_q[:-1], I_q, sI_q, sq_q)) - - datasets = [] - fileNumber = get_flist(numberString)[-1] - fileName = fileNameCreator(dataPath, fileNumber)[0] - header = ort_header(fileName) - orso_data=fileio.OrsoDataset(header, I_q.T) - datasets.append(orso_data) - fileio.save_orso(datasets, 'e2h.ort', data_separator='\n') - -#============================================================================== -#============================================================================== -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] - global foldback - foldback = not clas.noTOFCorrection - - #================================ - # 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) - - if foldback: - tof_e = np.where(tof_etau+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("-b", "--noTOFCorrection", - action='store_true', - help ="do not correct tof of seond neutron pulse") - 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=-7.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("-r", "--iDonno", - action='store_true', - help ="no idea") - clas.add_argument("-s", "--spy", - action='store_true', - help ="report a few key values, no plotting or writing") - 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("-u", "--update", - action='store_true', - help ="update output every 5 seconds") - 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 , 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() - - if clas.update: - verbous = False - logging.basicConfig(level=logging.ERROR, format='# %(message)s') - logging.error('e2h: recreating "e2h_life.png" every 5 s') - delay = 5 - output = 'e2h_life.png' - oldtime = 0 - while True: - dataPath = get_directDataPath(clas) - fileName = fileNameCreator(dataPath, 0)[0] - newtime = os.path.getmtime(fileName) - if newtime-oldtime: - print('\r# processing (press ^C to stop)', end='', flush=True) - process(dataPath, '0', clas) - oldtime = newtime - for i in range(5): - print('\r# waiting'+'.'*i+' '*(5-i)+' (press ^C to stop)', end='', flush=True) - time.sleep(delay/5) - signal.signal(signal.SIGINT, endIt) - else: - 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() - - diff --git a/eos.py b/eos.py deleted file mode 100644 index 61a2232..0000000 --- a/eos.py +++ /dev/null @@ -1,46 +0,0 @@ -#!/usr/bin/env python -""" -eos reduces measurements performed on Amor@SINQ, PSI - -Author: Jochen Stahn (algorithms, python draft), - Artur Glavic (structuring and optimisation of code) - -conventions (not strictly followed, yet): -- array names end with the suffix '_x[y]' with the meaning - _e = events - _tof - _l = lambda - _t = theta - _z = detector z - _lz = (lambda, detector z) - _q = q_z -""" - -import logging - -from libeos.command_line import command_line_options -from libeos.logconfig import setup_logging -from libeos.reduction import AmorReduction - -#===================================================================================================== -# TODO: -# - calculate resolution using the chopperPhase -# - deal with background correction -# - format of 'call' + add '-Y' if not supplied -#===================================================================================================== - -def main(): - setup_logging() - logging.warning('######## eos - data reduction for Amor ########') - - # read command line arguments and generate classes holding configuration parameters - config = command_line_options() - # Create reducer with these arguments - reducer = AmorReduction(config) - # Perform actual reduction - reducer.reduce() - - logging.info('######## eos - finished ########') - -if __name__ == '__main__': - main() diff --git a/libeos/__init__.py b/eos/__init__.py similarity index 51% rename from libeos/__init__.py rename to eos/__init__.py index 983b1c7..8c93e50 100644 --- a/libeos/__init__.py +++ b/eos/__init__.py @@ -1,7 +1,6 @@ """ -Package to handle data redction at AMOR instrument to be used by eos.py script. +Package to handle data redction at AMOR instrument to be used by __main__.py script. """ -__version__ = '2.2.0' -__date__ = '2025-09-16' - +__version__ = '3.0.0' +__date__ = '2025-10-06' diff --git a/eos/__main__.py b/eos/__main__.py new file mode 100644 index 0000000..ee4bf45 --- /dev/null +++ b/eos/__main__.py @@ -0,0 +1,41 @@ +""" +eos reduces measurements performed on Amor@SINQ, PSI + +Author: Jochen Stahn (algorithms, python draft), + Artur Glavic (structuring and optimisation of code) +""" + +import logging + +# need to do absolute import here as pyinstaller requires it +from eos.options import EOSConfig, ReaderConfig, ExperimentConfig, ReductionConfig, OutputConfig +from eos.command_line import commandLineArgs +from eos.logconfig import setup_logging, update_loglevel + +def main(): + setup_logging() + + # read command line arguments and generate classes holding configuration parameters + clas = commandLineArgs([ReaderConfig, ExperimentConfig, ReductionConfig, OutputConfig], + 'eos') + update_loglevel(clas.verbose) + + reader_config = ReaderConfig.from_args(clas) + experiment_config = ExperimentConfig.from_args(clas) + reduction_config = ReductionConfig.from_args(clas) + output_config = OutputConfig.from_args(clas) + config = EOSConfig(reader_config, experiment_config, reduction_config, output_config) + + logging.warning('######## eos - data reduction for Amor ########') + + # only import heavy module if sufficient command line parameters were provided + from eos.reduction import AmorReduction + # Create reducer with these arguments + reducer = AmorReduction(config) + # Perform actual reduction + reducer.reduce() + + logging.info('######## eos - finished ########') + +if __name__ == '__main__': + main() diff --git a/eos/command_line.py b/eos/command_line.py new file mode 100644 index 0000000..12ad312 --- /dev/null +++ b/eos/command_line.py @@ -0,0 +1,39 @@ +import argparse + +from typing import List +from .options import ArgParsable + + +def commandLineArgs(config_items: List[ArgParsable], program_name=None): + """ + Process command line argument. + The type of the default values is used for conversion and validation. + """ + msg = "eos reads data from (one or several) raw file(s) of the .hdf format, \ + performs various corrections, conversations and projections and exports\ + the resulting reflectivity in an orso-compatible format." + clas = argparse.ArgumentParser(prog=program_name, + description = msg, formatter_class=argparse.ArgumentDefaultsHelpFormatter) + + clas.add_argument('-v', '--verbose', action='count', default=0) + + clas_groups = {} + + all_arguments = [] + for cls in config_items: + all_arguments += cls.get_commandline_parameters() + + all_arguments.sort() # parameters are sorted alphabetically, unless they have higher priority + for cpc in all_arguments: + if not cpc.group in clas_groups: + clas_groups[cpc.group] = clas.add_argument_group(cpc.group) + if cpc.short_form: + clas_groups[cpc.group].add_argument( + f'-{cpc.short_form}', f'--{cpc.argument}', **cpc.add_argument_args + ) + else: + clas_groups[cpc.group].add_argument( + f'--{cpc.argument}', **cpc.add_argument_args + ) + + return clas.parse_args() diff --git a/libeos/const.py b/eos/const.py similarity index 84% rename from libeos/const.py rename to eos/const.py index 27cd24a..b7dc7ce 100644 --- a/libeos/const.py +++ b/eos/const.py @@ -4,3 +4,4 @@ Constants used in data reduction. hdm = 6.626176e-34/1.674928e-27 # h / m lamdaCut = 2.5 # Aa +lamdaMax = 15.0 # Aa \ No newline at end of file diff --git a/eos/event_analysis.py b/eos/event_analysis.py new file mode 100644 index 0000000..250434e --- /dev/null +++ b/eos/event_analysis.py @@ -0,0 +1,121 @@ +""" +Define an event dataformat that performs reduction actions like wavelength calculation on per-event basis. +With large number of events these actions can be time consuming so they use numba based functions. +""" +import numpy as np +import logging + +from typing import Tuple + +from . import const +from .event_data_types import EventDataAction, EventDatasetProtocol, append_fields, EVENT_BITMASKS +from .helpers import filter_project_x, merge_frames, extract_walltime +from .instrument import Detector +from .options import IncidentAngle +from .header import Header + +class ExtractWalltime(EventDataAction): + def perform_action(self, dataset: EventDatasetProtocol) ->None: + wallTime = extract_walltime(dataset.data.events.tof, + dataset.data.packets.start_index, + dataset.data.packets.Time) + logging.debug(f' expending event stream by wallTime') + new_events = append_fields(dataset.data.events, [('wallTime', wallTime.dtype)]) + new_events.wallTime = wallTime + dataset.data.events = new_events + +class MergeFrames(EventDataAction): + def perform_action(self, dataset: EventDatasetProtocol)->None: + tofCut = const.lamdaCut*dataset.geometry.chopperDetectorDistance/const.hdm*1e-13 + total_offset = (tofCut + + dataset.timing.tau * (dataset.timing.ch1TriggerPhase + dataset.timing.chopperPhase/2)/180) + dataset.data.events.tof = merge_frames(dataset.data.events.tof, tofCut, dataset.timing.tau, total_offset) + + +class AnalyzePixelIDs(EventDataAction): + def __init__(self, yRange: Tuple[int, int]): + self.yRange = yRange + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + d = dataset.data + (detZ, detXdist, delta, mask) = filter_project_x( + Detector.pixelLookUp, d.events.pixelID, self.yRange[0], self.yRange[1] + ) + ana_events = append_fields(d.events, [ + ('detZ', detZ.dtype), ('detXdist', detXdist.dtype), ('delta', delta.dtype)]) + # add analysis per event + ana_events.detZ = detZ + ana_events.detXdist = detXdist + ana_events.delta = delta + ana_events.mask += np.logical_not(mask)*EVENT_BITMASKS['yRange'] + d.events = ana_events + +class CalculateWavelength(EventDataAction): + def __init__(self, lambdaRange: Tuple[float, float]): + self.lambdaRange = lambdaRange + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + d = dataset.data + if not 'detXdist' in dataset.data.events.dtype.names: + raise ValueError("CalculateWavelength requires dataset with analyzed pixels, perform AnalyzePixelIDs first") + + #lamdaMax = const.lamdaCut+1.e13*dataset.timing.tau*const.hdm/(dataset.geometry.chopperDetectorDistance+124.) + + # lambda + # TODO: one of the most time consuming actions, could be implemented in numba, instead? + lamda = (1.e13*const.hdm)*d.events.tof/(dataset.geometry.chopperDetectorDistance+d.events.detXdist) + + final_events = append_fields(d.events, [('lamda', np.float64)]) + # add analysis per event + final_events.lamda = lamda + final_events.mask += EVENT_BITMASKS["LamdaRange"]*( + (self.lambdaRange[0]>lamda) | (lamda>self.lambdaRange[1])) + d.events = final_events + +class CalculateQ(EventDataAction): + def __init__(self, incidentAngle: IncidentAngle): + self.incidentAngle = incidentAngle + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + d = dataset.data + if not 'lamda' in dataset.data.events.dtype.names: + raise ValueError("CalculateQ requires dataset with analyzed wavelength, perform CalculateWavelength first") + + lamda = d.events.lamda + + final_events = append_fields(d.events, [('qz', np.float64)]) + + # alpha_f + # q_z + if self.incidentAngle == IncidentAngle.alphaF: + alphaF_e = dataset.geometry.nu - dataset.geometry.mu + d.events.delta + final_events.qz = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/lamda) + elif self.incidentAngle == IncidentAngle.nu: + alphaF_e = (dataset.geometry.nu + d.events.delta + dataset.geometry.kap + dataset.geometry.kad) / 2. + final_events.qz = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/lamda) + else: + alphaF_e = dataset.geometry.nu - dataset.geometry.mu + d.events.delta + alphaI = dataset.geometry.kap + dataset.geometry.kad + dataset.geometry.mu + final_events.qz = 2*np.pi * ((np.sin(np.deg2rad(alphaF_e)) + np.sin(np.deg2rad(alphaI)))/lamda) + final_events = append_fields(final_events, [('qx', np.float64)]) + final_events.qx = 2*np.pi * ((np.cos(np.deg2rad(alphaF_e)) - np.cos(np.deg2rad(alphaI)))/lamda) + + dataset.data.events = final_events + + def update_header(self, header: Header): + if self.incidentAngle == IncidentAngle.alphaF: + header.measurement_scheme = 'angle- and energy-dispersive' + else: + header.measurement_scheme = 'energy-dispersive' + +class FilterQzRange(EventDataAction): + def __init__(self, qzRange: Tuple[float, float]): + self.qzRange = qzRange + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + d = dataset.data + if not 'qz' in dataset.data.events.dtype.names: + raise ValueError("FilterQzRange requires dataset with qz values per events, perform WavelengthAndQ first") + + if self.qzRange[1]<0.5: + d.events.mask += EVENT_BITMASKS["qRange"]*((self.qzRange[0]>d.events.qz) | (d.events.qz>self.qzRange[1])) diff --git a/eos/event_data_types.py b/eos/event_data_types.py new file mode 100644 index 0000000..2faf297 --- /dev/null +++ b/eos/event_data_types.py @@ -0,0 +1,144 @@ +""" +Specify the data type and protocol used for event datasets. +""" +from typing import List, Optional, Protocol, Tuple +from dataclasses import dataclass +from .header import Header +from abc import ABC, abstractmethod +from hashlib import sha256 +import numpy as np +import logging + +@dataclass +class AmorGeometry: + mu:float + nu:float + kap:float + kad:float + div:float + + chopperSeparation: float + detectorDistance: float + chopperDetectorDistance: float + +@dataclass +class AmorTiming: + ch1TriggerPhase: float + ch2TriggerPhase: float + chopperSpeed: float + chopperPhase: float + tau: float + +# Structured datatypes used for event streams +EVENT_TYPE = np.dtype([('tof', np.float64), ('pixelID', np.uint32), ('mask', np.int32)]) +PACKET_TYPE = np.dtype([('start_index', np.uint32), ('Time', np.int64)]) +PULSE_TYPE = np.dtype([('time', np.int64), ('monitor', np.float32)]) +PC_TYPE = np.dtype([('current', np.float32), ('time', np.int64)]) + +# define the bitmask for individual event filters +EVENT_BITMASKS = { + 'MonitorThreshold': 1, + 'StrangeTimes': 2, + 'yRange': 4, + 'LamdaRange': 8, + 'qRange': 16, + } + +def append_fields(input: np.recarray, new_fields: List[Tuple[str, np.dtype]]): + # add one ore more fields to a recarray, numpy functions seems to fail + flds = [(name, dtypei[0]) for name, dtypei in input.dtype.fields.items()] + flds += new_fields + output = np.recarray(len(input), dtype=flds) + + for field in input.dtype.fields.keys(): + output[field] = input[field] + return output + +@dataclass +class AmorEventStream: + events: np.recarray # EVENT_TYPE + packets: np.recarray # PACKET_TYPE + pulses: np.recarray # PULSE_TYPE + proton_current: np.recarray # PC_TYPE + +class EventDatasetProtocol(Protocol): + """ + Minimal attributes a dataset needs to provide to work with EventDataAction + """ + geometry: AmorGeometry + timing: AmorTiming + data: AmorEventStream + + def append(self, other): + # Should define a way to add events from other to own + ... + + def update_header(self, header:Header): + # update a header with the information read from file + ... + +class EventDataAction(ABC): + """ + Abstract base class used for actions applied to an EventDatasetProtocol based objects. + Each action can optionally modify the header information. + Actions can be combined using the pipe operator | (OR). + """ + + def __call__(self, dataset: EventDatasetProtocol)->None: + logging.debug(f" Enter action {self.__class__.__name__} on {dataset!r}") + self.perform_action(dataset) + + @abstractmethod + def perform_action(self, dataset: EventDatasetProtocol)->None: ... + + def update_header(self, header:Header)->None: + if hasattr(self, 'action_name'): + header.reduction.corrections.append(getattr(self, 'action_name')) + + def __or__(self, other:'EventDataAction')->'CombinedAction': + return CombinedAction([self, other]) + + def __repr__(self): + output = self.__class__.__name__+'(' + for key,value in self.__dict__.items(): + output += f'{key}={value}, ' + return output.rstrip(', ')+')' + + def action_hash(self)->bytes: + # generate a unique hash that encodes this action with its configuration parameters + mh = sha256() + mh.update(self.__class__.__name__.encode()) + for key,value in sorted(self.__dict__.items()): + mh.update(repr(value).encode()) + return mh.hexdigest() + +class CombinedAction(EventDataAction): + """ + Used to perform multiple actions in one call. Stores a sequence of actions + that are then performed individually one after the other. + """ + def __init__(self, actions: List[EventDataAction]) -> None: + self._actions = actions + + def perform_action(self, dataset: EventDatasetProtocol)->None: + for action in self._actions: + action(dataset) + + def update_header(self, header:Header)->None: + for action in self._actions: + action.update_header(header) + + def __or__(self, other:'EventDataAction')->'CombinedAction': + return CombinedAction(self._actions+[other]) + + def __repr__(self): + output = repr(self._actions[0]) + for ai in self._actions[1:]: + output += ' | '+repr(ai) + return output + + def action_hash(self)->bytes: + mh = sha256() + for action in self._actions: + mh.update(action.action_hash().encode()) + return mh.hexdigest() diff --git a/eos/event_handling.py b/eos/event_handling.py new file mode 100644 index 0000000..0ad806a --- /dev/null +++ b/eos/event_handling.py @@ -0,0 +1,179 @@ +""" +Calculations performed on AmorEventData. +This module contains actions that do not need the numba base helper functions. Other actions are in event_analysis +""" +import logging +import os +import numpy as np + +from .header import Header +from .options import ExperimentConfig, MonitorType +from .event_data_types import EventDatasetProtocol, EventDataAction, EVENT_BITMASKS + +class ApplyPhaseOffset(EventDataAction): + def __init__(self, chopperPhaseOffset: float): + self.chopperPhaseOffset=chopperPhaseOffset + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + logging.debug( + f' replaced ch1TriggerPhase = {dataset.timing.ch1TriggerPhase} ' + f'with {self.chopperPhaseOffset}') + dataset.timing.ch1TriggerPhase = self.chopperPhaseOffset + +class ApplyParameterOverwrites(EventDataAction): + def __init__(self, config: ExperimentConfig): + self.config=config + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + if self.config.muOffset: + logging.debug(f' set muOffset = {self.config.muOffset}') + dataset.geometry.mu += self.config.muOffset + if self.config.mu: + logging.debug(f' replaced mu = {dataset.geometry.mu} with {self.config.mu}') + dataset.geometry.mu = self.config.mu + if self.config.nu: + logging.debug(f' replaced nu = {dataset.geometry.nu} with {self.config.nu}') + dataset.geometry.nu = self.config.nu + logging.info(f' mu = {dataset.geometry.mu:6.3f}, ' + f'nu = {dataset.geometry.nu:6.3f}, ' + f'kap = {dataset.geometry.kap:6.3f}, ' + f'kad = {dataset.geometry.kad:6.3f}') + + def update_header(self, header:Header) ->None: + if self.config.sampleModel: + import yaml + from orsopy.fileio.model_language import SampleModel + if self.config.sampleModel.endswith('.yml') or self.config.sampleModel.endswith('.yaml'): + if os.path.isfile(self.config.sampleModel): + with open(self.config.sampleModel, 'r') as model_yml: + model = yaml.safe_load(model_yml) + else: + logging.warning(f' ! the file {self.config.sampleModel}.yml does not exist. Ignored!') + return + else: + model = dict(stack=self.config.sampleModel) + logging.debug(f' set sample.model = {self.config.sampleModel}') + header.sample.model = SampleModel.from_dict(model) + + +class CorrectChopperPhase(EventDataAction): + def perform_action(self, dataset: EventDatasetProtocol) ->None: + dataset.data.events.tof += dataset.timing.tau*(dataset.timing.ch1TriggerPhase-dataset.timing.chopperPhase/2)/180 + + +class CorrectSeriesTime(EventDataAction): + def __init__(self, seriesStartTime): + self.seriesStartTime = np.int64(seriesStartTime) + + def perform_action(self, dataset: EventDatasetProtocol)->None: + if not 'wallTime' in dataset.data.events.dtype.names: + raise ValueError("CorrectTimeSeries requires walltTime to be extracted, please run ExtractWalltime first") + dataset.data.pulses.time -= self.seriesStartTime + dataset.data.events.wallTime -= self.seriesStartTime + dataset.data.proton_current.time -= self.seriesStartTime + start, stop = dataset.data.events.wallTime[0], dataset.data.events.wallTime[-1] + logging.debug(f' wall time from {start/1e9:6.1f} s to {stop/1e9:6.1f} s, ' + f'series time = {self.seriesStartTime/1e9:6.1f}') + + +class AssociatePulseWithMonitor(EventDataAction): + def __init__(self, monitorType:MonitorType, lowCurrentThreshold:float): + self.monitorType = monitorType + self.lowCurrentThreshold = lowCurrentThreshold + + def perform_action(self, dataset: EventDatasetProtocol)->None: + logging.debug(f' using monitor type {self.monitorType}') + if self.monitorType in [MonitorType.proton_charge or MonitorType.debug]: + if not 'wallTime' in dataset.data.events.dtype.names: + raise ValueError( + "AssociatePulseWithMonitor requires walltTime to be extracted, please run ExtractWalltime first") + monitorPerPulse = self.get_current_per_pulse(dataset.data.pulses.time, + dataset.data.proton_current.time, + dataset.data.proton_current.current)\ + * 2*dataset.timing.tau * 1e-3 + # filter low-current pulses + dataset.data.pulses.monitor = np.where( + monitorPerPulse > 2*dataset.timing.tau * self.lowCurrentThreshold * 1e-3, + monitorPerPulse, 0) + elif self.monitorType==MonitorType.time: + dataset.data.pulses.monitor = 2*dataset.timing.tau + else: # pulses + dataset.data.pulses.monitor = 1 + + if self.monitorType == MonitorType.debug: + cpp, t_bins = np.histogram(dataset.data.events.wallTime, dataset.data.pulses.time) + np.savetxt('tme.hst', np.vstack((dataset.data.pulses.time[:-1], cpp, dataset.data.pulses.monitor[:-1])).T) + + if self.monitorType in [MonitorType.proton_charge or MonitorType.debug]: + self.monitor_threshold(dataset) + + def monitor_threshold(self, dataset): + goodTimeS = dataset.data.pulses.time[dataset.data.pulses.monitor!=0] + filter_e = np.logical_not(np.isin(dataset.data.events.wallTime, goodTimeS)) + dataset.data.events.mask += EVENT_BITMASKS['MonitorThreshold']*filter_e + logging.info(f' low-beam (<{self.lowCurrentThreshold} mC) rejected pulses: ' + f'{dataset.data.pulses.monitor.shape[0]-goodTimeS.shape[0]} ' + f'out of {dataset.data.pulses.monitor.shape[0]}') + logging.info(f' with {filter_e.sum()} events') + if goodTimeS.shape[0]: + logging.info(f' average counts per pulse = {dataset.data.events.shape[0]/goodTimeS.shape[0]:7.1f}') + else: + logging.info(f' average counts per pulse = undefined') + + @staticmethod + def get_current_per_pulse(pulseTimeS, currentTimeS, currents): + # add currents for early pulses and current time value after last pulse (j+1) + currentTimeS = np.hstack([[0], currentTimeS, [pulseTimeS[-1]+1]]) + currents = np.hstack([[0], currents]) + pulseCurrentS = np.zeros(pulseTimeS.shape[0], dtype=float) + j = 0 + for i, ti in enumerate(pulseTimeS): + # find the last current item that was before this pulse + while ti >= currentTimeS[j+1]: + j += 1 + pulseCurrentS[i] = currents[j] + return pulseCurrentS + + +class FilterStrangeTimes(EventDataAction): + def perform_action(self, dataset: EventDatasetProtocol)->None: + filter_e = np.logical_not(dataset.data.events.tof<=2*dataset.timing.tau) + dataset.data.events.mask += EVENT_BITMASKS['StrangeTimes']*filter_e + if filter_e.any(): + logging.warning(f' strange times: {filter_e.sum()}') + + +class TofTimeCorrection(EventDataAction): + def __init__(self, correct_chopper_opening: bool = True): + self.correct_chopper_opening = correct_chopper_opening + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + d = dataset.data + if self.correct_chopper_opening: + d.events.tof -= ( d.events.delta / 180. ) * dataset.timing.tau + else: + d.events.tof -= ( dataset.geometry.kad / 180. ) * dataset.timing.tau + +class ApplyMask(EventDataAction): + def __init__(self, bitmask_filter=None): + self.bitmask_filter = bitmask_filter + + def perform_action(self, dataset: EventDatasetProtocol) ->None: + # TODO: why is this action time consuming? + d = dataset.data + pre_filter = d.events.shape[0] + if logging.getLogger().level == logging.DEBUG: + # only run this calculation if debug level is actually active + filtered_by_mask = {} + for key, value in EVENT_BITMASKS.items(): + filtered_by_mask[key] = ((d.events.mask & value)!=0).sum() + logging.debug(f" Removed by filters: {filtered_by_mask}") + if self.bitmask_filter is None: + d.events = d.events[d.events.mask==0] + else: + # remove the provided bitmask_filter bits from the events + # this means that all bits that are set in bitmask_filter will NOT be used to filter events + fltr = (d.events.mask & (~self.bitmask_filter)) == 0 + d.events = d.events[fltr] + post_filter = d.events.shape[0] + logging.info(f' number of events: total = {pre_filter:7d}, filtered = {post_filter:7d}') diff --git a/eos/file_reader.py b/eos/file_reader.py new file mode 100644 index 0000000..598647e --- /dev/null +++ b/eos/file_reader.py @@ -0,0 +1,345 @@ +""" +Reading of Amor NeXus data files to extract metadata and event stream. +""" +from typing import BinaryIO, List, Union + +import h5py +import numpy as np +import platform +import logging +import subprocess + +from datetime import datetime + +from orsopy import fileio +from orsopy.fileio.model_language import SampleModel + +from . import const +from .header import Header +from .event_data_types import AmorGeometry, AmorTiming, AmorEventStream, PACKET_TYPE, EVENT_TYPE, PULSE_TYPE, PC_TYPE + +try: + import zoneinfo +except ImportError: + # for python versions < 3.9 try to use the backports version + from backports import zoneinfo + + +# Time zone used to interpret time strings +AMOR_LOCAL_TIMEZONE = zoneinfo.ZoneInfo(key='Europe/Zurich') + +if platform.node().startswith('amor'): + NICOS_CACHE_DIR = '/home/amor/nicosdata/amor/cache/' + GREP = '/usr/bin/grep "%s"' +else: + NICOS_CACHE_DIR = None + +class AmorEventData: + """ + Read one amor NeXus datafile and extract relevant header information. + + Implements EventDatasetProtocol + """ + file_list: List[str] + first_index: int + last_index: int = -1 + EOF: bool = False + max_events: int + owner: fileio.Person + experiment: fileio.Experiment + sample: fileio.Sample + instrument_settings: fileio.InstrumentSettings + geometry: AmorGeometry + timing: AmorTiming + data: AmorEventStream + + eventStartTime: np.int64 + + def __init__(self, fileName:Union[str, h5py.File, BinaryIO], first_index:int=0, max_events:int=100_000_000): + if type(fileName) is str: + self.file_list = [fileName] + self.hdf = h5py.File(fileName, 'r', swmr=True) + elif type(fileName) is h5py.File: + self.file_list = [fileName.filename] + self.hdf = fileName + else: + self.file_list = [repr(fileName)] + self.hdf = h5py.File(fileName, 'r') + self.first_index = first_index + self.max_events = max_events + + self.read_header_info() + self.read_instrument_configuration() + self.read_event_stream() + + # actions applied to any dataset + self.read_chopper_trigger_stream() + self.read_proton_current_stream() + + if type(fileName) is str: + # close the input file to free memory, only if the file was opened in this object + self.hdf.close() + del(self.hdf) + + def _replace_if_missing(self, key, nicos_key, dtype=float): + try: + return dtype(self.hdf[f'/entry1/Amor/{key}'][0]) + except(KeyError, IndexError): + if NICOS_CACHE_DIR: + try: + logging.warning(f" using parameter {nicos_key} from nicos cache") + year_date = self.fileDate.strftime('%Y') + value = str(subprocess.getoutput(f'{GREP} {NICOS_CACHE_DIR}nicos-{nicos_key}/{year_date}')).split('\t')[-1] + return dtype(value) + except Exception: + logging.error("Couldn't get value from nicos cache", exc_info=True) + return dtype(0) + else: + logging.warning(f" parameter {key} not found, relpace by zero") + return dtype(0) + + def read_header_info(self): + # read general information and first data set + title = self.hdf['entry1/title'][0].decode('utf-8') + proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') + user_name = self.hdf['entry1/user/name'][0].decode('utf-8') + user_affiliation = 'unknown' + user_email = self.hdf['entry1/user/email'][0].decode('utf-8') + user_orcid = None + sampleName = self.hdf['entry1/sample/name'][0].decode('utf-8') + model = self.hdf['entry1/sample/model'][0].decode('utf-8') + if 'stack:' in model: + import yaml + model = yaml.safe_load(model) + else: + model = dict(stack=model) + instrumentName = 'Amor' + source = self.hdf['entry1/Amor/source/name'][0].decode('utf-8') + sourceProbe = 'neutron' + start_time = self.hdf['entry1/start_time'][0].decode('utf-8') + # extract start time as unix time, adding UTC offset of 1h to time string + start_date = datetime.fromisoformat(start_time) + self.fileDate = start_date.replace(tzinfo=AMOR_LOCAL_TIMEZONE) + + self.owner = fileio.Person( + name=user_name, + affiliation=user_affiliation, + contact=user_email, + ) + if user_orcid: + self.owner.orcid = user_orcid + + self.experiment = fileio.Experiment( + title=title, + instrument=instrumentName, + start_date=start_date, + probe=sourceProbe, + facility=source, + proposalID=proposal_id + ) + self.sample = fileio.Sample( + name=sampleName, + model=SampleModel.from_dict(model), + sample_parameters=None, + ) + + def read_instrument_configuration(self): + chopperSeparation = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0)) + detectorDistance = float(np.take(self.hdf['entry1/Amor/detector/transformation/distance'], 0)) + chopperDetectorDistance = detectorDistance-float(np.take(self.hdf['entry1/Amor/chopper/distance'], 0)) + + polarizationConfigs = ['unpolarized', 'unpolarized', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm'] + + mu = self._replace_if_missing('instrument_control_parameters/mu', 'mu', float) + nu = self._replace_if_missing('instrument_control_parameters/nu', 'nu', float) + kap = self._replace_if_missing('instrument_control_parameters/kappa', 'kappa', float) + kad = self._replace_if_missing('instrument_control_parameters/kappa_offset', 'kad', float) + div = self._replace_if_missing('instrument_control_parameters/div', 'div', float) + ch1TriggerPhase = self._replace_if_missing('chopper/ch1_trigger_phase', 'ch1_trigger_phase', float) + ch2TriggerPhase = self._replace_if_missing('chopper/ch2_trigger_phase', 'ch2_trigger_phase', float) + try: + chopperTriggerTime = (float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][7]) \ + -float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][0])) \ + /7 + chopperTriggerTimeDiff = float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][2]) + except (KeyError, IndexError): + logging.debug(' chopper speed and phase taken from .hdf file') + chopperSpeed = self._replace_if_missing('chopper/rotation_speed', 'chopper_phase', float) + chopperPhase = self._replace_if_missing('chopper/phase', 'chopper_phase', float) + tau = 30/chopperSpeed + else: + tau = int(1e-6*chopperTriggerTime/2+0.5)*(1e-3) + chopperTriggerPhase = 180e-9*chopperTriggerTimeDiff/tau + chopperSpeed = 30/tau + chopperPhase = chopperTriggerPhase+ch1TriggerPhase-ch2TriggerPhase + + self.geometry = AmorGeometry(mu, nu, kap, kad, div, + chopperSeparation, detectorDistance, chopperDetectorDistance) + self.timing = AmorTiming(ch1TriggerPhase, ch2TriggerPhase, chopperSpeed, chopperPhase, tau) + + polarizationConfigLabel = self._replace_if_missing('polarization/configuration/average_value', 'polarizer_config_label', int) + polarizationConfig = fileio.Polarization(polarizationConfigs[polarizationConfigLabel]) + logging.debug(f' polarization configuration: {polarizationConfig} (index {polarizationConfigLabel})') + + + self.instrument_settings = fileio.InstrumentSettings( + incident_angle = fileio.ValueRange(round(mu+kap+kad-0.5*div, 3), + round(mu+kap+kad+0.5*div, 3), + 'deg'), + wavelength = fileio.ValueRange(const.lamdaCut, const.lamdaMax, 'angstrom'), + #polarization = fileio.Polarization.unpolarized, + polarization = fileio.Polarization(polarizationConfig) + ) + self.instrument_settings.mu = fileio.Value( + round(mu, 3), + 'deg', + comment='sample angle to horizon') + self.instrument_settings.nu = fileio.Value( + round(nu, 3), + 'deg', + comment='detector angle to horizon') + self.instrument_settings.div = fileio.Value( + round(div, 3), + 'deg', + comment='incoming beam divergence') + self.instrument_settings.kap = fileio.Value( + round(kap, 3), + 'deg', + comment='incoming beam inclination') + if abs(kad)>0.02: + self.instrument_settings.kad = fileio.Value( + round(kad, 3), + 'deg', + comment='incoming beam angular offset') + + + def update_header(self, header:Header): + """ + Add dataset information into an existing header. + """ + logging.info(f' meta data from: {self.file_list[0]}') + header.owner = self.owner + header.experiment = self.experiment + header.sample = self.sample + header.measurement_instrument_settings = self.instrument_settings + + def read_event_stream(self): + """ + Read the actual event data from file. If file is too large, find event index from packets + that allow splitting of file smaller than self.max_events. + """ + packets = np.recarray(self.hdf['/entry1/Amor/detector/data/event_index'].shape, dtype=PACKET_TYPE) + packets.start_index = self.hdf['/entry1/Amor/detector/data/event_index'][:] + packets.Time = self.hdf['/entry1/Amor/detector/data/event_time_zero'][:] + try: + # packet index that matches first event index + start_packet = int(np.where(packets.start_index==self.first_index)[0][0]) + except IndexError: + raise IndexError(f'No event packet found starting at event #{self.first_index}') + packets = packets[start_packet:] + + nevts = self.hdf['/entry1/Amor/detector/data/event_time_offset'].shape[0] + if (nevts-self.first_index)>self.max_events: + end_packet = np.where(packets.start_index<=(self.first_index+self.max_events))[0][-1] + self.last_index = packets.start_index[end_packet]-1 + packets = packets[:end_packet] + else: + self.last_index = nevts-1 + self.EOF = True + nevts = self.last_index+1-self.first_index + + # adapte packet to event index relation + packets.start_index -= self.first_index + + events = np.recarray(nevts, dtype=EVENT_TYPE) + events.tof = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][self.first_index:self.last_index+1])/1.e9 + events.pixelID = self.hdf['/entry1/Amor/detector/data/event_id'][self.first_index:self.last_index+1] + events.mask = 0 + + pulses = self.read_chopper_trigger_stream() + current = self.read_proton_current_stream() + self.data = AmorEventStream(events, packets, pulses, current) + + if self.first_index>0 and not self.EOF: + # label the file name if not all events were used + self.file_list[0] += f'[{self.first_index}:{self.last_index}]' + + def read_chopper_trigger_stream(self): + chopper1TriggerTime = np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][:-2], dtype=np.int64) + #self.chopper2TriggerTime = self.chopper1TriggerTime + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time'][:-2], dtype=np.int64) + # + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][:], dtype=np.int64) + if np.shape(chopper1TriggerTime)[0] > 2: + startTime = chopper1TriggerTime[0] + pulseTimeS = chopper1TriggerTime + else: + logging.warn(' no chopper trigger data available, using event steram instead') + startTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][0], dtype=np.int64) + stopTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-2], dtype=np.int64) + pulseTimeS = np.arange(startTime, stopTime, self.timing.tau*1e9, dtype=np.int64) + pulses = np.recarray(pulseTimeS.shape, dtype=PULSE_TYPE) + pulses.time = pulseTimeS + pulses.monitor = 1. # default is monitor pulses as it requires no calculation + # apply filter in case the events were filtered + if self.first_index>0 or not self.EOF: + pulses = pulses[(pulses.time>=self.data.packets.Time[0])&(pulses.time<=self.data.packets.Time[-1])] + self.eventStartTime = startTime + return pulses + + def read_proton_current_stream(self): + proton_current = np.recarray(self.hdf['entry1/Amor/detector/proton_current/time'].shape, dtype=PC_TYPE) + proton_current.time = self.hdf['entry1/Amor/detector/proton_current/time'][:] + proton_current.current = self.hdf['entry1/Amor/detector/proton_current/value'][:,0] + if self.first_index>0 or not self.EOF: + proton_current = proton_current[(proton_current.time>=self.data.packets.Time[0])& + (proton_current.time<=self.data.packets.Time[-1])] + return proton_current + + def info(self): + output = "" + for key in ['owner', 'experiment', 'sample', 'instrument_settings']: + value = repr(getattr(self, key)).replace("\n","\n ") + output += f'\n{key}={value},' + output += '\n' + return output + + def append(self, other): + """ + Append event streams from another file to this one. Adjusts the event indices in the + packets to stay valid. + """ + new_events = np.concatenate([self.data.events, other.data.events]).view(np.recarray) + new_pulses = np.concatenate([self.data.pulses, other.data.pulses]).view(np.recarray) + new_proton_current = np.concatenate([self.data.proton_current, other.data.proton_current]).view(np.recarray) + new_packets = np.concatenate([self.data.packets, other.data.packets]).view(np.recarray) + new_packets.start_index[self.data.packets.shape[0]:] += self.data.events.shape[0] + self.data = AmorEventStream(new_events, new_packets, new_pulses, new_proton_current) + # Indicate that this is amodified dataset, basically counts number of appends as negative indices + self.last_index = min(self.last_index-1, -1) + self.file_list += other.file_list + + def __repr__(self): + output = (f"AmorEventData({self.file_list!r}) # {self.data.events.shape[0]} events, " + f"{self.data.pulses.shape[0]} pulses") + + return output + + def get_timeslice(self, start, end)->'AmorEventData': + # return a new dataset with just events that occured in given time slice + if not 'wallTime' in self.data.events.dtype.names: + raise ValueError("This dataset is missing a wallTime that is required for time slicing") + # convert from seconds to epoch integer values + start , end = start*1e9, end*1e9 + event_filter = self.data.events.wallTime>=start + event_filter &= self.data.events.wallTime=start + pulse_filter &= self.data.pulses.time=self.qzRange[0]] return q_grid + @cache def lamda(self): lamdaMax = 16 - lamdaMin = self.lamdaCut + lamdaMin = const.lamdaCut lamda_grid = lamdaMin*(1+self.dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+self.dldl)+1)) return lamda_grid + @cache def z(self): return np.arange(Detector.nBlades*Detector.nWires+1) + @cache def lz(self): - return np.ones(( np.shape(self.lamda()[:-1])[0], np.shape(self.z()[:-1])[0] )) + return np.ones(( self.lamda().shape[0]-1, self.z().shape[0]-1)) + @cache def delta(self, detectorDistance): # unused for now bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*Detector.bladeZ / detectorDistance) ) diff --git a/libeos/logconfig.py b/eos/logconfig.py similarity index 94% rename from libeos/logconfig.py rename to eos/logconfig.py index 6f80945..b93b104 100644 --- a/libeos/logconfig.py +++ b/eos/logconfig.py @@ -33,10 +33,10 @@ def setup_logging(): logfile.setLevel(logging.DEBUG) logger.addHandler(logfile) -def update_loglevel(verbose=False, debug=False): - if verbose: +def update_loglevel(verbose=0): + if verbose==1: logging.getLogger().handlers[0].setLevel(logging.INFO) - if debug: + if verbose>1: console = logging.getLogger().handlers[0] console.setLevel(logging.DEBUG) formatter = logging.Formatter('%(levelname).1s %(message)s') diff --git a/eos/normalisation.py b/eos/normalisation.py new file mode 100644 index 0000000..b1e7b33 --- /dev/null +++ b/eos/normalisation.py @@ -0,0 +1,79 @@ +""" +Defines how to normalize a focusing reflectometry dataset by a reference measurement. +""" +import logging +import os +import numpy as np +from typing import List, Optional + +from .event_data_types import EventDatasetProtocol +from .header import Header +from .options import NormalisationMethod +from .instrument import Detector, LZGrid + + +class LZNormalisation: + file_list = List[str] + angle: float + monitor: float + norm: np.ndarray + + def __init__(self, reference:EventDatasetProtocol, normalisationMethod: NormalisationMethod, grid: LZGrid): + self.angle = reference.geometry.nu-reference.geometry.mu + lamda_e = reference.data.events.lamda + detZ_e = reference.data.events.detZ + self.monitor = np.sum(reference.data.pulses.monitor) + norm_lz, _, _ = np.histogram2d(lamda_e, detZ_e, bins=(grid.lamda(), grid.z())) + norm_lz = np.where(norm_lz>2, norm_lz, np.nan) + if normalisationMethod==NormalisationMethod.direct_beam: + # TODO: move flipping to projection + self.norm = np.flip(norm_lz, 1) + else: + # correct for reference sm reflectivity + lamda_l = grid.lamda() + theta_z = self.angle+Detector.delta_z + lamda_lz = (grid.lz().T*lamda_l[:-1]).T + theta_lz = grid.lz()*theta_z + qz_lz = 4.0*np.pi*np.sin(np.deg2rad(theta_lz))/lamda_lz + # TODO: introduce variable for `m` and propably for the slope + # Correct reflectivity of m=5 supermirror + Rsm_lz = np.ones(np.shape(qz_lz)) + Rsm_lz = np.where(qz_lz>0.0217, 1-(qz_lz-0.0217)*(0.0625/0.0217), Rsm_lz) + Rsm_lz = np.where(qz_lz>0.0217*5, np.nan, Rsm_lz) + self.norm = norm_lz/Rsm_lz + self.file_list = [os.path.basename(entry) for entry in reference.file_list] + + @classmethod + def from_file(cls, filename, check_hash=None) -> Optional['LZNormalisation']: + self = super().__new__(cls) + with open(filename, 'rb') as fh: + hash = str(np.load(fh, allow_pickle=True)) + self.file_list = np.load(fh, allow_pickle=True) + self.angle = np.load(fh, allow_pickle=True) + self.norm = np.load(fh, allow_pickle=True) + self.monitor = np.load(fh, allow_pickle=True) + if check_hash is not None and hash != check_hash: + logging.info(' file hash does not match this reduction configuration') + raise ValueError('file hash does not match this reduction configuration') + return self + + @classmethod + def unity(cls, grid:LZGrid) -> 'LZNormalisation': + logging.warning(f'normalisation is unity') + self = super().__new__(cls) + self.norm = grid.lz() + self.file_list = [] + self.angle = 1. + self.monitor = 1. + return self + + def safe(self, filename, hash): + with open(filename, 'wb') as fh: + np.save(fh, hash, allow_pickle=False) + np.save(fh, np.array(self.file_list), allow_pickle=False) + np.save(fh, np.array(self.angle), allow_pickle=False) + np.save(fh, self.norm, allow_pickle=False) + np.save(fh, self.monitor, allow_pickle=False) + + def update_header(self, header:Header): + header.measurement_additional_files = self.file_list diff --git a/eos/options.py b/eos/options.py new file mode 100644 index 0000000..3975994 --- /dev/null +++ b/eos/options.py @@ -0,0 +1,576 @@ +""" +Classes for stroing various configurations needed for reduction. +""" +import argparse +from dataclasses import dataclass, field, Field, fields, MISSING +from enum import StrEnum +from typing import get_args, get_origin, List, Optional, Tuple, Union +from datetime import datetime +from os import path +import numpy as np + +import logging + +@dataclass +class CommandlineParameterConfig: + argument: str # default parameter for command line resutign ins "--argument" + add_argument_args: dict # all arguments that will be passed to add_argument method + short_form: Optional[str] = None + group: str = 'misc' + priority: int = 0 + + def __gt__(self, other): + """ + Sort required arguments first, then use priority, then name + """ + return (not self.add_argument_args.get('required', False), -self.priority, self.argument)>( + not other.add_argument_args.get('required', False), -other.priority, other.argument) + +class ArgParsable: + def __init_subclass__(cls): + # create a nice documentation string that takes help into account + cls.__doc__ = cls.__name__ + " Parameters:\n" + for key, typ in cls.__annotations__.items(): + if get_origin(typ) is Union and type(None) in get_args(typ): + optional = True + typ = get_args(typ)[0] + else: + optional = False + + value = getattr(cls, key, None) + cls.__doc__ += f" {key} ({typ.__name__})" + if isinstance(value, Field): + if value.default is not MISSING: + cls.__doc__ += f" = {value.default}" + if 'help' in value.metadata: + cls.__doc__ += f" - {value.metadata['help']}" + elif value is not None: + cls.__doc__ += f" = {value}" + if optional: + cls.__doc__ += " [Optional]" + cls.__doc__ += "\n" + return cls + + @classmethod + def get_commandline_parameters(cls) -> List[CommandlineParameterConfig]: + """ + Return a list of arguments used in building the command line parameters. + + Union types besides Optional are not supported. + """ + output = [] + for field in fields(cls): + args={} + if field.default is not MISSING: + args['default'] = field.default + args['required'] = False + elif field.default_factory is not MISSING: + args['default'] = field.default_factory() + args['required'] = False + else: + args['required'] = True + if get_origin(field.type) is Union and type(None) in get_args(field.type): + # optional argument + typ = get_args(field.type)[0] + del(args['default']) + else: + typ = field.type + if get_origin(typ) is list: + args['nargs'] = '+' + typ = get_args(typ)[0] + elif get_origin(typ) is tuple: + args['nargs'] = len(get_args(typ)) + typ = get_args(typ)[0] + + if issubclass(typ, StrEnum): + args['choices'] = [ci.value for ci in typ] + if field.default is not MISSING: + args['default'] = field.default.value + typ = str + + if typ is bool: + args['action'] = 'store_false' if field.default else 'store_true' + else: + args['type'] = typ + + if 'help' in field.metadata: + args['help'] = field.metadata['help'] + + output.append(CommandlineParameterConfig( + field.name, + add_argument_args=args, + group=field.metadata.get('group', 'misc'), + short_form=field.metadata.get('short', None), + priority=field.metadata.get('priority', 0), + )) + return output + + @classmethod + def get_default(cls, key): + """ + Return the default argument for an attribute, None if it doesn't exist. + """ + for field in fields(cls): + if field.name != key: + continue + if field.default is not MISSING: + return field.default + elif field.default_factory is not MISSING: + return field.default_factory() + return None + + def is_default(self, key): + value = getattr(self, key) + return value == self.get_default(key) + + @classmethod + def from_args(cls, args: argparse.Namespace): + """ + Create the child class from the command line argument Namespace object. + All attributes that are not needed for this class are ignored. + """ + inpargs = {} + for field in fields(cls): + value = getattr(args, field.name) + typ = field.type + if get_origin(field.type) is Union and type(None) in get_args(field.type): + # optional argument + typ = get_args(field.type)[0] + + if isinstance(typ, type) and issubclass(typ, StrEnum): + # convert str to enum + try: + value = typ(value) + except ValueError: + choices = [ci.value for ci in typ] + raise ValueError(f"Parameter --{field.name} has to be one of {choices}") + + inpargs[field.name] = value + return cls(**inpargs) + +# definition of command line arguments + +@dataclass +class ReaderConfig(ArgParsable): + year: int = field( + default=datetime.now().year, + metadata={ + 'short': 'Y', + 'group': 'input data', + 'help': 'year the measurement was performed', + }, + ) + rawPath: List[str] = field( + default_factory=lambda: ['.', path.join('.','raw'), path.join('..','raw'), path.join('..','..','raw')], + metadata={ + 'short': 'rp', + 'group': 'input data', + 'help': 'search paths for hdf files', + }, + ) + startTime: Optional[float] = field( + default = None, + metadata={ + 'short': 'st', + 'group': 'data manicure', + 'help': 'set time zero other than the start of the data aquisition', + }, + ) + +class IncidentAngle(StrEnum): + alphaF = 'alphaF' + mu = 'mu' + nu = 'nu' + +class MonitorType(StrEnum): + auto = 'a' + proton_charge = 'p' + time = 't' + neutron_monitor = 'n' + debug = 'x' + +@dataclass +class ExperimentConfig(ArgParsable): + chopperPhase: float = field( + default=0, + metadata={ + 'short': 'cp', + 'group': 'instrument settings', + 'help': 'phase between opening of chopper 1 and closing of chopper 2 window', + }, + ) + chopperPhaseOffset: float = field( + default=-5, + metadata={ + 'short': 'co', + 'group': 'instrument settings', + 'help': 'phase between chopper 1 index pulse and closing edge', + }, + ) + chopperSpeed: float = field( + default=500, + metadata={ + 'short': 'cs', + 'group': 'instrument settings', + 'help': 'rotation speed of the chopper disks in rpm', + }, + ) + yRange: Tuple[int, int] = field( + default=(18, 48), + metadata={ + 'short': 'y', + 'group': 'region of interest', + 'help': 'horizontal pixel range on the detector to be used', + }, + ) + lambdaRange: Tuple[float, float] = field( + default_factory=lambda: [3, 12.5], + metadata={ + 'short': 'l', + 'group': 'region of interest', + 'help': 'wavelength range to be used (in angstrom)', + }, + ) + lowCurrentThreshold: float = field( + default=50, + metadata={ + 'short': 'pt', + 'group': 'instrument settings', + 'help': 'proton current below which the events are ignored (per chopper pulse)', + }, + ) + + incidentAngle: IncidentAngle = field( + default=IncidentAngle.alphaF, + metadata={ + 'short': 'ai', + 'group': 'instrument settings', + 'help': 'calculate alphaI = [alphaF], [mu]+kappa+delta_kappa or ([nu]+kappa+delta_kappa)/2', + }, + ) + alphaF = 'alphaF' + sampleModel: Optional[str] = field( + default=None, + metadata={ + 'short': 'sm', + 'group': 'sample', + 'help': 'orso type string to describe the sample in one line', + }, + ) + mu: Optional[float] = field( + default=None, + metadata={ + 'short': 'mu', + 'group': 'sample', + 'help': 'inclination of the sample surface w.r.t. the instrument horizon', + }, + ) + nu: Optional[float] = field( + default=None, + metadata={ + 'short': 'nu', + 'group': 'sample', + 'help': 'inclination of the detector w.r.t. the instrument horizon', + }, + ) + muOffset: Optional[float] = field( + default=0, + metadata={ + 'short': 'm', + 'group': 'sample', + 'help': 'correction offset for mu misalignment (mu_real = mu_file + mu_offset)', + }, + ) + monitorType: MonitorType = field( + default=MonitorType.proton_charge, + metadata={ + 'short': 'mt', + 'group': 'instrument settings', + 'help': 'one of [a]uto, [p]rotonCurrent, [t]ime or [n]eutronMonitor', + }, + ) + +class NormalisationMethod(StrEnum): + direct_beam = 'd' + over_illuminated = 'o' + under_illuminated = 'u' + +@dataclass +class ReductionConfig(ArgParsable): + fileIdentifier: List[str] = field( + metadata={ + 'short': 'f', + 'priority': 100, + 'group': 'input data', + 'help': 'file number(s) or offset (if < 1)', + }, + ) + + qResolution: float = field( + default=0.01, + metadata={ + 'short': 'r', + 'group': 'data manicure', + 'help': 'output resolution of q-scale Delta q / q', + }, + ) + qzRange: Tuple[float, float] = field( + default_factory=lambda: [0.005, 0.51], + metadata={ + 'short': 'q', + 'group': 'region of interest', + 'help': '?', + }, + ) + thetaRange: Tuple[float, float] = field( + default_factory=lambda: [-12., 12.], + metadata={ + 'short': 't', + 'group': 'region of interest', + 'help': 'absolute theta region of interest', + }, + ) + thetaRangeR: Tuple[float, float] = field( + default_factory=lambda: [-0.75, 0.75], + metadata={ + 'short': 'T', + 'group': 'region of interest', + 'help': 'theta region of interest w.r.t. beam center', + }, + ) + normalisationMethod: NormalisationMethod = field( + default=NormalisationMethod.over_illuminated, + metadata={ + 'short': 'nm', + 'priority': 90, + 'group': 'input data', + 'help': 'normalisation method: [o]verillumination, [u]nderillumination, [d]irect_beam'}) + scale: List[float] = field( + default_factory=lambda: [1.], + metadata={ + 'short': 's', + 'group': 'data manicure', + 'help': '(list of) scaling factors, if less elements than files use the last one', + }, + ) + autoscale: Tuple[float, float] = field( + default=None, + metadata={ + 'short': 'S', + 'group': 'data manicure', + 'help': '', + }, + ) + subtract: Optional[str] = field( + default=None, + metadata={ + 'short': 'sub', + 'group': 'input data', + 'help': 'File with R(q_z) curve to be subtracted (in .Rqz.ort format)'}) + normalisationFileIdentifier: Optional[List[str]] = field( + default_factory=lambda: [None], + metadata={ + 'short': 'n', + 'priority': 90, + 'group': 'input data', + 'help': 'file number(s) of normalisation measurement'}) + timeSlize: Optional[List[float]] = field( + default= None, + metadata={ + 'short': 'ts', + 'group': 'region of interest', + 'help': 'time slizing ,[ [,stop]]', + }, + ) + + def _expand_file_list(self, short_notation:str): + """Evaluate string entry for file number lists""" + file_list=[] + for i in short_notation.split(','): + if '-' in i: + if ':' in i: + step = i.split(':', 1)[1] + file_list += range(int(i.split('-', 1)[0]), + int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, + int(step)) + else: + step = 1 + file_list += range(int(i.split('-', 1)[0]), + int(i.split('-', 1)[1])+1, + int(step)) + else: + file_list += [int(i)] + file_list.sort() + return file_list + + def data_files(self): + # get input files from expanding fileIdentifier + return list(map(self._expand_file_list, self.fileIdentifier)) + + def normal_files(self): + return list(map(self._expand_file_list, self.normalisationFileIdentifier)) + + +class OutputFomatOption(StrEnum): + Rqz_ort = "Rqz.ort" + Rqz_orb = "Rqz.orb" + Rlt_ort = "Rlt.ort" + Rlt_orb = "Rlt.orb" + ort = "ort" + orb = "orb" + Rqz = "Rqz" + Rlt = "Rlt" + + +@dataclass +class OutputConfig(ArgParsable): + outputFormats: List[OutputFomatOption] = field( + default_factory=lambda: ['Rqz.ort'], + metadata={ + 'short': 'of', + 'group': 'output', + 'help': 'one of "Rqz[.ort]", "Rlt[.ort]" or both with "ort"', + }, + ) + outputName: str = field( + default='fromEOS', + metadata={ + 'short': 'o', + 'group': 'output', + 'help': '?', + }, + ) + outputPath: str = field( + default='.', + metadata={ + 'short': 'op', + 'group': 'output', + 'help': '?', + }, + ) + + def _output_format_list(self, outputFormat): + format_list = [] + if OutputFomatOption.ort in outputFormat\ + or OutputFomatOption.Rqz_ort in outputFormat\ + or OutputFomatOption.Rqz in outputFormat: + format_list.append(OutputFomatOption.Rqz_ort) + if OutputFomatOption.ort in outputFormat\ + or OutputFomatOption.Rlt_ort in outputFormat\ + or OutputFomatOption.Rlt in outputFormat: + format_list.append(OutputFomatOption.Rlt_ort) + if OutputFomatOption.orb in outputFormat\ + or OutputFomatOption.Rqz_orb in outputFormat\ + or OutputFomatOption.Rqz in outputFormat: + format_list.append(OutputFomatOption.Rqz_orb) + if OutputFomatOption.orb in outputFormat\ + or OutputFomatOption.Rlt_orb in outputFormat\ + or OutputFomatOption.Rlt in outputFormat: + format_list.append(OutputFomatOption.Rlt_orb) + return sorted(format_list, reverse=True) + + def __post_init__(self): + self.outputFormats = self._output_format_list(self.outputFormats) + + +# =================================== + +@dataclass +class EOSConfig: + reader: ReaderConfig + experiment: ExperimentConfig + reduction: ReductionConfig + output: OutputConfig + + _call_string_overwrite=None + + #@property + #def call_string(self)->str: + # if self._call_string_overwrite: + # return self._call_string_overwrite + # else: + # return self.calculate_call_string() + + def call_string(self): + base = 'eos' + + inpt = '' + if self.reader.year: + inpt += f' -Y {self.reader.year}' + else: + inpt += f' -Y {datetime.now().year}' + if np.shape(self.reader.rawPath)[0] == 1: + inpt += f' --rawPath {self.reader.rawPath}' + if self.reduction.subtract: + inpt += f' -subtract {self.reduction.subtract}' + if self.reduction.normalisationFileIdentifier: + inpt += f' -n {" ".join(self.reduction.normalisationFileIdentifier)}' + if self.reduction.fileIdentifier: + inpt += f' -f {" ".join(self.reduction.fileIdentifier)}' + + otpt = '' + if self.reduction.qResolution: + otpt += f' -r {self.reduction.qResolution}' + if self.output.outputPath != '.': + inpt += f' --outputdPath {self.output.outputPath}' + if self.output.outputName: + otpt += f' -o {self.output.outputName}' + if self.output.outputFormats != ['Rqz.ort']: + otpt += f' -of {" ".join(self.output.outputFormats)}' + + mask = '' + # TODO: Check if you want these parameters for the case of default call + mask += f' -y {" ".join(str(ii) for ii in self.experiment.yRange)}' + mask += f' -l {" ".join(str(ff) for ff in self.experiment.lambdaRange)}' + mask += f' -t {" ".join(str(ff) for ff in self.reduction.thetaRange)}' + mask += f' -T {" ".join(str(ff) for ff in self.reduction.thetaRangeR)}' + mask += f' -q {" ".join(str(ff) for ff in self.reduction.qzRange)}' + + para = '' + # TODO: Check if we want these parameters for defaults + para += f' --chopperPhase {self.experiment.chopperPhase}' + para += f' --chopperPhaseOffset {self.experiment.chopperPhaseOffset}' + if self.experiment.mu: + para += f' --mu {self.experiment.mu}' + elif self.experiment.muOffset: + para += f' --muOffset {self.experiment.muOffset}' + if self.experiment.nu: + para += f' --nu {self.experiment.nu}' + + modl = '' + if self.experiment.sampleModel: + modl += f" --sampleModel '{self.experiment.sampleModel}'" + + acts = '' + if self.reduction.autoscale: + acts += f' --autoscale {" ".join(str(ff) for ff in self.reduction.autoscale)}' + # TODO: Check if should be shown if not default + acts += f' --scale {self.reduction.scale}' + if self.reduction.timeSlize: + acts += f' --timeSlize {" ".join(str(ff) for ff in self.reduction.timeSlize)}' + + mlst = base + inpt + otpt + if mask: + mlst += mask + if para: + mlst += para + if acts: + mlst += acts + if modl: + mlst += modl + + if len(mlst) > 70: + mlst = base + ' ' + inpt + ' ' + otpt + if mask: + mlst += ' ' + mask + if para: + mlst += ' ' + para + if acts: + mlst += ' ' + acts + if modl: + mlst += ' ' + modl + + logging.debug(f'Argument list build in EOSConfig.call_string: {mlst}') + return mlst + + diff --git a/eos/path_handling.py b/eos/path_handling.py new file mode 100644 index 0000000..e401d22 --- /dev/null +++ b/eos/path_handling.py @@ -0,0 +1,49 @@ +""" +Defines how file paths are resolved from short_notation, year and number to filename. +""" +import os +from typing import List + + +class PathResolver: + def __init__(self, year, rawPath): + self.year = year + self.rawPath = rawPath + + def resolve(self, short_notation): + return list(map(self.get_path, self.expand_file_list(short_notation))) + + @staticmethod + def expand_file_list(short_notation)->List[int]: + """Evaluate string entry for file number lists""" + file_list = [] + for i in short_notation.split(','): + if '-' in i: + if ':' in i: + step = i.split(':', 1)[1] + file_list += range(int(i.split('-', 1)[0]), + int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, + int(step)) + else: + step = 1 + file_list += range(int(i.split('-', 1)[0]), + int(i.split('-', 1)[1])+1, + int(step)) + else: + file_list += [int(i)] + return list(sorted(file_list)) + + def get_path(self, number): + fileName = f'amor{self.year}n{number:06d}.hdf' + path = '' + for rawd in self.rawPath: + if os.path.exists(os.path.join(rawd, fileName)): + path = rawd + break + if not path: + if os.path.exists( + f'/afs/psi.ch/project/sinqdata/{self.year}/amor/{int(number/1000)}/{fileName}'): + path = f'/afs/psi.ch/project/sinqdata/{self.year}/amor/{int(number/1000)}' + else: + raise FileNotFoundError(f'# ERROR: the file {fileName} can not be found in {self.rawPath}') + return os.path.join(path, fileName) diff --git a/eos/projection.py b/eos/projection.py new file mode 100644 index 0000000..e6a7c5b --- /dev/null +++ b/eos/projection.py @@ -0,0 +1,285 @@ +""" +Classes used to calculate projections/binnings from event data onto given grids. +""" + +import logging +import numpy as np +from dataclasses import dataclass + +from .event_data_types import EventDatasetProtocol +from .instrument import Detector, LZGrid +from .normalisation import LZNormalisation + +@dataclass +class ProjectedReflectivity: + R: np.ndarray + dR: np.ndarray + Q: np.ndarray + dQ: np.ndarray + + @property + def data(self): + """ + Return combined data compatible with storing as columns in orso file. + Q, R, dR, dQ + """ + return np.array([self.Q, self.R, self.dR, self.dQ]).T + + def data_for_time(self, time): + tme = np.ones(np.shape(self.Q))*time + return np.array([self.Q, self.R, self.dR, self.dQ, tme]).T + + def scale(self, factor): + self.R *= factor + self.dR *= factor + + def autoscale(self, range): + filter_q = (range[0]<=self.Q) & (self.Q<=range[1]) + filter_q &= self.dR>0 + if filter_q.sum()>0: + scale = (self.R[filter_q]/self.dR[filter_q]).sum()/(self.R[filter_q]**2/self.dR[filter_q]).sum() + self.scale(scale) + logging.info(f' scaling factor = {scale}') + return scale + else: + logging.warning(' automatic scaling not possible') + return 1.0 + + def stitch(self, other: 'ProjectedReflectivity'): + # find scaling factor between two reflectivities at points both are not zero + filter_q = np.logical_not(np.isnan(other.R*self.R)) + filter_q &= self.R>0 + filter_q &= other.R>0 + R1 = self.R[filter_q] + dR1 = self.dR[filter_q] + R2 = other.R[filter_q] + dR2 = other.dR[filter_q] + if len(R1)>0: + scale = (R1**2*R2**2/(dR1**2*dR2**2)).sum() / (R1**3*R2/(dR1**2*dR2**2)).sum() + self.scale(scale) + logging.info(f' scaling factor = {scale}') + return scale + else: + logging.warning(' automatic scaling not possible') + return 1.0 + + def subtract(self, R, dR): + # subtract another dataset with same q-points + self.R -= R + self.dR = np.sqrt(self.dR**2+dR**2) + +class LZProjection: + grid: LZGrid + lamda: np.ndarray + alphaF: np.ndarray + is_normalized: bool + + data: np.recarray + + def __init__(self, tthh: float, grid: LZGrid): + self.grid = grid + self.is_normalized = False + + alphaF_z = tthh + Detector.delta_z + lamda_l = self.grid.lamda() + lamda_c = (lamda_l[:-1]+lamda_l[1:])/2 + + lz_shape = self.grid.lz() + + self.lamda = lz_shape*lamda_c[:, np.newaxis] + self.alphaF = lz_shape*alphaF_z[np.newaxis, :] + self.data = np.zeros(self.alphaF.shape, dtype=[ + ('I', np.float64), + ('mask', bool), + ('ref', np.float64), + ('err', np.float64), + ('res', np.float64), + ('qz', np.float64), + ('qx', np.float64), + ('norm', np.float64), + ]).view(np.recarray) + self.data.mask = True + self.monitor = 0. + + @classmethod + def from_dataset(cls, dataset: EventDatasetProtocol, grid: LZGrid, has_offspecular=False): + tthh = dataset.geometry.nu - dataset.geometry.mu + output = cls(tthh, grid) + output.correct_gravity(dataset.geometry.detectorDistance) + if has_offspecular: + alphaI_lz = grid.lz()*(dataset.geometry.mu+dataset.geometry.kap+dataset.geometry.kad) + output.calculate_q(alphaI_lz) + else: + output.calculate_q() + return output + + def correct_gravity(self, detector_distance): + self.alphaF += np.rad2deg( np.arctan( 3.07e-10 * detector_distance * self.lamda**2 ) ) + + def calculate_q(self, alphaI=None): + if alphaI is None: + self.data.qz = 4.0*np.pi*np.sin(np.deg2rad(self.alphaF))/self.lamda + self.data.qx = 0.*self.data.qz + else: + self.data.qz = 2.0*np.pi*(np.sin(np.deg2rad(self.alphaF))+np.sin(np.deg2rad(alphaI)))/self.lamda + self.data.qx = 2.0*np.pi*(np.cos(np.deg2rad(self.alphaF))-np.cos(np.deg2rad(alphaI)))/self.lamda + + if self.data.qz[0,self.data.qz.shape[1]//2] < 0: + # assuming a 'measurement from below' when center of detector at negative qz + self.data.qz *= -1 + + self.calculate_q_resolution() + + def calculate_q_resolution(self): + res_lz = self.grid.lz() * 0.022**2 + res_lz = res_lz + (0.008/self.alphaF)**2 + self.data.res = self.data.qz * np.sqrt(res_lz) + + def apply_theta_filter(self, theta_range): + # Filters points within theta range + self.data.mask &= (self.alphaFtheta_range[1]) + + def apply_theta_mask(self, theta_range): + # Mask points outside theta range + self.data.mask &= self.alphaF>=theta_range[0] + self.data.mask &= self.alphaF<=theta_range[1] + + def apply_lamda_mask(self, lamda_range): + # Mask points outside lambda range + self.data.mask &= self.lamda>=lamda_range[0] + self.data.mask &= self.lamda<=lamda_range[1] + + def apply_norm_mask(self, norm: LZNormalisation): + # Mask points where normliazation is nan + self.data.mask &= np.logical_not(np.isnan(norm.norm)) + + def project(self, dataset: EventDatasetProtocol, monitor: float): + """ + Project dataset on grid and add to intensity. + Can be called multiple times to sequentially add events. + """ + e = dataset.data.events + int_lz, *_ = np.histogram2d(e.lamda, e.detZ, bins = (self.grid.lamda(), self.grid.z())) + self.data.I += int_lz + self.monitor += monitor + # in case the intensity changed one needs to normalize again + self.is_normalized = False + + @property + def I(self): + output = self.data.I[:] + output[np.logical_not(self.data.mask)] = np.nan + return output / self.monitor + + def calc_error(self): + # calculate error bars for resulting intensity after normalization + self.data.err = self.data.ref * np.sqrt( 1./(self.data.I+.1) + 1./self.data.norm ) + + def normalize_over_illuminated(self, norm: LZNormalisation): + """ + Normalize the dataaset and take into account a difference in + detector angle for measurement and reference. + """ + norm_lz = norm.norm + thetaN_z = Detector.delta_z+norm.angle + thetaN_lz = np.ones_like(norm_lz)*thetaN_z + thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) + self.data.mask &= (np.absolute(thetaN_lz)>5e-3) + ref_lz = (self.data.I*np.absolute(thetaN_lz))/(norm_lz*np.absolute(self.alphaF)) + ref_lz *= norm.monitor/self.monitor + ref_lz[np.logical_not(self.data.mask)] = np.nan + self.data.norm = norm_lz + self.data.ref = ref_lz + self.calc_error() + self.is_normalized = True + + def normalize_no_footprint(self, norm: LZNormalisation): + norm_lz = norm.norm + ref_lz = (self.data.I/norm_lz) + ref_lz *= norm.monitor/self.monitor + ref_lz[np.logical_not(self.data.mask)] = np.nan + self.data.norm = norm_lz + self.data.ref = ref_lz + self.calc_error() + self.is_normalized = True + + def scale(self, factor: float): + if not self.is_normalized: + raise ValueError("Dataset needs to be normalized, first") + self.data.ref *= factor + self.data.err *= factor + + def project_on_qz(self): + if not self.is_normalized: + raise ValueError("Dataset needs to be normalized, first") + q_q = self.grid.q() + weights_lzf = self.data.norm[self.data.mask] + q_lzf = self.data.qz[self.data.mask] + R_lzf = self.data.ref[self.data.mask] + dR_lzf = self.data.err[self.data.mask] + dq_lzf = self.data.res[self.data.mask] + + N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + + R_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf * R_lzf )[0] + R_q = R_q / N_q + + dR_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dR_lzf)**2 )[0] + dR_q = np.sqrt( dR_q ) / N_q + + # TODO: different error propagations for dR and dq! + # this is what should work: + #dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0] + #dq_q = np.sqrt( dq_q ) / N_q + # and this actually works: + N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf**2 )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0] + dq_q = np.sqrt( dq_q / N_q ) + + return ProjectedReflectivity(R_q, dR_q, (q_q[1:]+q_q[:-1])/2., dq_q) + + ############## potential speedup not used right now, needs to be tested #################### + @classmethod + def histogram2d_lz(cls, lamda_e, detZ_e, bins): + """ + Perform binning operation equivalent to numpy bin2d for the sepcial case + of the second dimension using integer positions (pre-defined pixels). + Based on the devide_bin algorithm below. + """ + dimension = bins[1].shape[0]-1 + if not (np.array(bins[1])==np.arange(dimension+1)).all(): + raise ValueError("histogram2d_lz requires second bin dimension to be contigous integer range") + binning = cls.devide_bin(lamda_e, detZ_e.astype(np.int64), bins[0], dimension) + return np.array(binning), bins[0], bins[1] + + @classmethod + def devide_bin(cls, lambda_e, position_e, lamda_edges, dimension): + ''' + Use a divide and conquer strategy to bin the data. For the actual binning the + numpy bincount function is used, as it is much faster than histogram for + counting of integer values. + + :param lambda_e: Array of wavelength for each event + :param position_e: Array of positional indices for each event + :param lamda_edges: The edges of bins to be used for the histogram + :param dimension: position number of buckets in output arrray + + :return: 2D list of dimensions (lambda, x) of counts + ''' + if len(lambda_e)==0: + # no more events in range, return empty bins + return [np.zeros(dimension, dtype=np.int64).tolist()]*(len(lamda_edges)-1) + if len(lamda_edges)==2: + # deepest recursion reached, all items should be within the two ToF edges + return [np.bincount(position_e, minlength=dimension).tolist()] + # split all events into two time of flight regions + split_idx = len(lamda_edges)//2 + left_region = lambda_e1: + self.experiment_config.monitorType = MonitorType.proton_charge + logging.debug(' monitor type set to "proton current"') + else: + self.experiment_config.monitorType = MonitorType.time + logging.debug(' monitor type set to "time"') + # update actions to sue selected monitor + self.prepare_actions() + # reload normalization to make sure the monitor matches + if self.reduction_config.normalisationFileIdentifier: + self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0]) + + self.dataevent_time_correction.seriesStartTime = self.dataset.eventStartTime + self.dataevent_actions(self.dataset) + self.dataset.update_header(self.header) + self.dataevent_actions.update_header(self.header) + for fi in file_list[1:]: + di = AmorEventData(fi) + self.dataevent_actions(di) + self.dataset.append(di) + + for fileName in file_list: + self.header.measurement_data_files.append(fileio.File( file=fileName.split('/')[-1], + timestamp=self.dataset.fileDate)) + + + if self.reduction_config.timeSlize: + if i>0: + logging.warning(" time slizing should only be used for on set of datafiles, check parameters") + self.analyze_timeslices(i) + else: + self.analyze_unsliced(i) + + def analyze_unsliced(self, i): + self.monitor = np.sum(self.dataset.data.pulses.monitor) + logging.warning(f' monitor = {self.monitor:8.2f} {MONITOR_UNITS[self.experiment_config.monitorType]}') + + proj:LZProjection = self.project_on_lz() + try: + scale = self.reduction_config.scale[i] + except IndexError: + scale = self.reduction_config.scale[-1] + proj.scale(scale) + + if 'Rqz.ort' in self.output_config.outputFormats: + headerRqz = self.header.orso_header() + headerRqz.data_set = f'Nr {i} : mu = {self.dataset.geometry.mu:6.3f} deg' + + # projection on qz-grid + result = proj.project_on_qz() + + if self.reduction_config.autoscale: + if i==0: + result.autoscale(self.reduction_config.autoscale) + else: + result.stitch(self.last_result) + + if self.subtract: + if len(result.Q)==len(self.sq_q): + result.subtract(self.sR_q, self.sdR_q) + else: + logging.warning( + f'backgroung file {self.sFileName} not compatible with q_z scale ({len(self.sq_q)} vs. {len(result.Q)})') + + orso_data = fileio.OrsoDataset(headerRqz, result.data) + self.last_result = result + self.datasetsRqz.append(orso_data) + if 'Rlt.ort' in self.output_config.outputFormats: + columns = [ + fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), + fileio.Column('R', '', 'specular reflectivity'), + fileio.ErrorColumn(error_of='R', error_type='uncertainty', value_is='sigma'), + fileio.ErrorColumn(error_of='Qz', error_type='resolution', value_is='sigma'), + fileio.Column('lambda', 'angstrom', 'wavelength'), + fileio.Column('alpha_f', 'deg', 'final angle'), + fileio.Column('l', '', 'index of lambda-bin'), + fileio.Column('t', '', 'index of theta bin'), + fileio.Column('intensity', '', 'filtered neutron events per pixel'), + fileio.Column('norm', '', 'normalisation matrix'), + fileio.Column('mask', '', 'pixels used for calculating R(q_z)'), + fileio.Column('Qx', '1/angstrom', 'parallel momentum transfer'), + ] + + ts, zs = proj.data.shape + lindex_lz = np.tile(np.arange(1, ts+1), (zs, 1)).T + tindex_lz = np.tile(np.arange(1, zs+1), (ts, 1)) + + j = 0 + for item in zip( + proj.data.qz.T, + proj.data.ref.T, + proj.data.err.T, + proj.data.res.T, + proj.lamda.T, + proj.alphaF.T, + lindex_lz.T, + tindex_lz.T, + proj.data.I.T, + proj.data.norm.T, + proj.data.mask.T, + proj.data.qx.T, + ): + data = np.array(list(item)).T + headerRlt = self.header.orso_header(columns=columns) + headerRlt.data_set = f'dataset_{i}_{j+1} : alpha_f = {proj.alphaF[0, j]:6.3f} deg' + orso_data = fileio.OrsoDataset(headerRlt, data) + self.datasetsRlt.append(orso_data) + j += 1 + + def analyze_timeslices(self, i): + wallTime_e = np.float64(self.dataset.data.events.wallTime)/1e9 + pulseTimeS = np.float64(self.dataset.data.pulses.time)/1e9 + interval = self.reduction_config.timeSlize[0] + try: + start = self.reduction_config.timeSlize[1] + except IndexError: + start = 0 + try: + stop = self.reduction_config.timeSlize[2] + except IndexError: + stop = wallTime_e[-1] + # make overwriting log lines possible by removing newline at the end + #logging.StreamHandler.terminator = "\r" + logging.warning(f' time slizing') + logging.info(' slize time monitor') + for ti, time in enumerate(np.arange(start, stop, interval)): + slice = self.dataset.get_timeslice(time, time+interval) + self.monitor = np.sum(slice.data.pulses.monitor) + logging.info(f' {ti:<4d} {time:6.0f} {self.monitor:7.2f} {MONITOR_UNITS[self.experiment_config.monitorType]}') + + proj: LZProjection = self.project_on_lz(slice) + try: + scale = self.reduction_config.scale[i] + except IndexError: + scale = self.reduction_config.scale[-1] + proj.scale(scale) + + # projection on qz-grid + result = proj.project_on_qz() + + if self.reduction_config.autoscale: + # scale every slice the same + if ti==0: + if i==0: + atscale = result.autoscale(self.reduction_config.autoscale) + else: + atscale = result.stitch(self.last_result) + else: + result.scale(atscale) + + if self.subtract: + if len(result.Q)==len(self.sq_q): + result.subtract(self.sR_q, self.sdR_q) + else: + logging.warning( + f'backgroung file {self.sFileName} not compatible with q_z scale ({len(self.sq_q)} vs. {len(result.Q)})') + + headerRqz = self.header.orso_header( + extra_columns=[fileio.Column('time', 's', 'time relative to start of measurement series')]) + headerRqz.data_set = f'{i}_{ti}: time = {time:8.1f} s to {time+interval:8.1f} s' + orso_data = fileio.OrsoDataset(headerRqz, result.data_for_time(time)) + self.datasetsRqz.append(orso_data) + self.last_result = result + # reset normal logging behavior + #logging.StreamHandler.terminator = "\n" + logging.info(f' done {min(time+interval, pulseTimeS[-1]):5.0f}') + + def save_Rqz(self): + fname = os.path.join(self.output_config.outputPath, f'{self.output_config.outputName}.Rqz.ort') + logging.warning(f' {fname}') + theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | R(q_z)' + fileio.save_orso(self.datasetsRqz, fname, data_separator='\n', comment=theSecondLine) + + def save_Rtl(self): + fname = os.path.join(self.output_config.outputPath, f'{self.output_config.outputName}.Rlt.ort') + logging.warning(f' {fname}') + theSecondLine = f' {self.header.experiment.title} | {self.header.experiment.start_date} | sample {self.header.sample.name} | R(lambda, theta)' + fileio.save_orso(self.datasetsRlt, fname, data_separator='\n', comment=theSecondLine) + + def loadRqz(self, name): + fname = os.path.join(self.output_config.outputPath, name) + if os.path.exists(fname): + fileName = fname + elif os.path.exists(f'{fname}.Rqz.ort'): + fileName = f'{fname}.Rqz.ort' + else: + sys.exit(f'### the background file \'{fname}\' does not exist! => stopping') + + q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) + + return q_q, Sq_q, dS_q, fileName + + def create_normalisation_map(self, short_notation): + outputPath = self.output_config.outputPath + normalisation_list = self.path_resolver.expand_file_list(short_notation) + name = '_'.join(map(str, normalisation_list)) + n_path = os.path.join(outputPath, f'{name}.norm') + + self.norm = None + if os.path.exists(n_path): + logging.debug(f'trying to load matrix from file {n_path}') + try: + self.norm = LZNormalisation.from_file(n_path, self.normevent_actions.action_hash()) + except (ValueError, EOFError): + self.norm =None + else: + logging.warning(f'normalisation matrix: found and using {n_path}') + if self.norm is None: + # in case file does not exist or the action hash doesn't match, create new normalization + logging.warning(f'normalisation matrix: using the files {normalisation_list}') + normalization_files = list(map(self.path_resolver.get_path, normalisation_list)) + reference = AmorEventData(normalization_files[0]) + self.normevent_actions(reference) + for nfi in normalization_files[1:]: + toadd = AmorEventData(nfi) + self.normevent_actions(toadd) + reference.append(toadd) + self.norm = LZNormalisation(reference, self.reduction_config.normalisationMethod, self.grid) + if reference.data.events.shape[0] > 1e6: + self.norm.safe(n_path, self.normevent_actions.action_hash()) + + self.header.measurement_additional_files = self.norm.file_list + self.header.reduction.corrections.append('normalisation with \'additional files\'') + + def project_on_lz(self, dataset=None): + if dataset is None: + dataset=self.dataset + proj = LZProjection.from_dataset(dataset, self.grid, + has_offspecular=(self.experiment_config.incidentAngle!=IncidentAngle.alphaF)) + + if not self.reduction_config.is_default('thetaRangeR'): + t0 = dataset.geometry.nu - dataset.geometry.mu + # adjust range based on detector center + thetaRange = [ti+t0 for ti in self.reduction_config.thetaRangeR] + proj.apply_theta_mask(thetaRange) + elif not self.reduction_config.is_default('thetaRange'): + proj.apply_theta_mask(self.reduction_config.thetaRange) + else: + thetaRange = [dataset.geometry.nu - dataset.geometry.mu - dataset.geometry.div/2, + dataset.geometry.nu - dataset.geometry.mu + dataset.geometry.div/2] + proj.apply_theta_mask(thetaRange) + + proj.apply_lamda_mask(self.experiment_config.lambdaRange) + + proj.apply_norm_mask(self.norm) + + proj.project(dataset, self.monitor) + + if self.reduction_config.normalisationMethod == NormalisationMethod.over_illuminated: + logging.debug(' assuming an overilluminated sample and correcting for the angle of incidence') + proj.normalize_over_illuminated(self.norm) + elif self.reduction_config.normalisationMethod==NormalisationMethod.under_illuminated: + logging.debug(' assuming an underilluminated sample and ignoring the angle of incidence') + proj.normalize_no_footprint(self.norm) + elif self.reduction_config.normalisationMethod==NormalisationMethod.direct_beam: + logging.debug(' assuming direct beam for normalisation and ignoring the angle of incidence') + proj.normalize_no_footprint(self.norm) + else: + logging.error('unknown normalisation method! Use [u]nder, [o]ver or [d]irect illumination') + proj.normalize_no_footprint(self.norm) + if self.monitor<=1e-6: + logging.info(' low monitor -> nan output') + proj.data.ref *= np.nan + + return proj diff --git a/events2histogram.py b/events2histogram.py index 9ecc8ee..29a8d6c 100755 --- a/events2histogram.py +++ b/events2histogram.py @@ -1,11 +1,10 @@ -__version__ = '2024-03-15' +__version__ = '2025-06-07' -# essential changes with regard to 2022 version -# - imprved orso header -# - nexus compatible -# - new theta grid -# - new content in data_e (angleas rather than distances) -# - bug fixes: tof correction within detector +# essential changes with regard to 2024 version +# - accepts new hdf structure +# TODO: +# - output path +# - solve the confusion for negative file numbers and the connecting '-' import os import sys @@ -88,7 +87,7 @@ def analyse_ev(event_e, tof_e, yMin, yMax, thetaMin, thetaMax): # correct tof for beam size effect at chopper # t_cor = (delta / 180 deg) * tau - data_e[:,0] -= ( data_e[:,5] / 180. ) * tau + data_e[:,0] -= ( data_e[:,5] / 180. ) * tau # effective flight path length #data_e[:,6] = chopperDetectorDistance + data_e[:,6] @@ -145,22 +144,23 @@ class Meta: 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)) + try: + self.mu = float(np.take(fh['/entry1/Amor/instrument_control_parameters/mu'], 0)) + self.nu = float(np.take(fh['/entry1/Amor/instrument_control_parameters/nu'], 0)) + self.kap = float(np.take(fh['/entry1/Amor/instrument_control_parameters/kappa'], 0)) + self.kad = float(np.take(fh['/entry1/Amor/instrument_control_parameters/kappa_offset'], 0)) + self.div = float(np.take(fh['/entry1/Amor/instrument_control_parameters/div'], 0)) + chopperSpeed = float(np.take(fh['/entry1/Amor/chopper/rotation_speed'], 0)) + chopperPhase = float(np.take(fh['/entry1/Amor/chopper/phase'], 0)) + ch1TriggerPhase = float(np.take(fh['/entry1/Amor/chopper/ch1_trigger_phase'], 0)) + polarizationConfigLabel = int(np.take(fh['/entry1/Amor/polarization/configuration/value'], 0)) + #polarizationConfigLabel = 1 except (KeyError, IndexError): logging.warning(f" using parameters from nicos cache") - #cachePath = '/home/amor/nicosdata/amor/cache/' - #cachePath = '/home/nicos/amorcache/' - cachePath = '/home/amor/cache/' + cachePath = '/home/amor/nicosdata/amor/cache/' + year_date = str(datetime.today()).split(' ')[0].replace("-", "/", 1) 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] @@ -172,31 +172,23 @@ class Meta: 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 + chopperSpeed = float(value) + #value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-chopper_speed'+year_date)).split('\t')[-1] + value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-chopper_phase/'+year_date)).split('\t')[-1] + chopperPhase = float(value) + value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-ch1_trigger_phase/'+year_date)).split('\t')[-1] + ch1TriggerPhase = float(value) + value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-polarization_config_label/'+year_date)).split('\t')[-1] + polarizationConfigLabel = int(value) - if chSp: - self.tau = 30. / chSp - else: - self.tau = 0 + self.tau = 30. / chopperSpeed - 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 = '?' + self.chopperTriggerPhase = ch1TriggerPhase + chopperPhase/2 + print(f'# chopper trigger phase = {ch1TriggerPhase}') + polarizationConfigs = ['undefined', 'oo', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm'] + self.polarizationConfig = polarizationConfigs[int(polarizationConfigLabel)] + #self.polarizationConfig = 'po' # for .ort header @@ -391,19 +383,13 @@ class PlotSelection: # header / meta data - def header(self, filename, mu, nu, totalCounts, countingTime, spin): + def header(self, filename, mu, nu, totalCounts, countingTime, polarizationConfig): 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) + header = "#{} \u03bc={:>1.2f} \u03bd={:>1.2f} {} {:>10} cts {:>8.1f} s".format(number, mu+5e-3, nu+5e-3, polarizationConfig, 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) + def headline(self, numberString, totalCounts, polarizationConfig): + headLine = "#{} \u03bc={:>1.2f} \u03bd={:>1.2f} p={} {:>12,} cts {:>8.1f} s".format(numberString, mu+5e-3, nu+5e-3, polarizationConfig, totalCounts, countingTime) return headLine # grids @@ -446,7 +432,7 @@ class PlotSelection: # create PNG with several plots - def all(self, numberString, arg, data_e): + def all(self, numberString, arg, data_e, polarizationConfig): #cmap='gist_earth' cmap = mpl.cm.gnuplot(np.arange(256)) cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1]) @@ -501,7 +487,7 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=2.8, c='r') fig.colorbar(cb, ax=mlt) plt.subplots_adjust(hspace=0.6, wspace=0.1) @@ -510,7 +496,7 @@ class PlotSelection: # create PNG with one plot - def Iyz(self, numberString, arg, data_e): + def Iyz(self, numberString, arg, data_e, polarizationConfig): det = Detector() cmap = mpl.cm.gnuplot(np.arange(256)) cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1]) @@ -533,13 +519,13 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.colorbar() plt.savefig(output, format='png', dpi=150) #plt.close() - def Ilt(self, numberString, arg, data_e) : + def Ilt(self, numberString, arg, data_e, polarizationConfig) : 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]) @@ -568,13 +554,13 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.colorbar() plt.savefig(output, format='png', dpi=150) #plt.close() - def Itz(self, numberString, arg, data_e): + def Itz(self, numberString, arg, data_e, polarizationConfig): det = Detector() cmap = mpl.cm.gnuplot(np.arange(256)) cmap[:1, :] = np.array([256/256, 255/256, 236/256, 1]) @@ -604,13 +590,13 @@ class PlotSelection: plt.ylim(0,) plt.xlabel('$t ~/~ \\mathrm{s}$') plt.ylabel('$z$ pixel row') - headline = self.headline(numberString, np.shape(data_e)[0]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.colorbar() plt.savefig(output, format='png', dpi=150) #plt.close() - def Iq(self, numberString, arg, data_e): + def Iq(self, numberString, arg, data_e, polarizationConfig): 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 ), @@ -634,12 +620,12 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.savefig(output, format='png', dpi=150) #plt.close() - def Il(self, numberString, arg, data_e): + def Il(self, numberString, arg, data_e, polarizationConfig): I_l, bins_l = np.histogram(data_e[:,7], bins = self.lamda_grid()) if arg == 'file': header = '# lambda counts' @@ -653,12 +639,12 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.savefig(output, format='png', dpi=150) #plt.close() - def It(self, numberString, arg, data_e): + def It(self, numberString, arg, data_e, polarizationConfig): I_t, bins_t = np.histogram(data_e[:,8], bins = self.theta_grid()) if arg == 'file': header = '# 2theta counts' @@ -669,12 +655,12 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.savefig(output, format='png', dpi=150) #plt.close() - def tof(self, numberString, arg, data_e): + def tof(self, numberString, arg, data_e, polarizationConfig): if foldback: time_grid = np.arange(0, 1.3*tau, 0.0005) else: @@ -696,7 +682,7 @@ class PlotSelection: 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]) + headline = self.headline(numberString, np.shape(data_e)[0], polarizationConfig) plt.title(headline, loc='left', y=1.0, c='r') plt.savefig(output, format='png') @@ -774,7 +760,7 @@ def process(dataPath, ident, clas): timeMin = 0 timeMax = 0 - chopperPhase = clas.chopperPhase + chopperTriggerPhase = clas.chopperTriggerPhase tofOffset = clas.TOFOffset thetaMin = clas.thetaRange[0] thetaMax = clas.thetaRange[1] @@ -840,29 +826,27 @@ def process(dataPath, ident, clas): tau = meta.tau try: - chPh + chopperTriggerPhase except NameError: - chPh = meta.chPh - spin = meta.spin + chopperTriggerPhase = meta.chopperTriggerPhase - global countingTime, detectorDistance, chopperDetectorDistance + + global countingTime, detectorDistance, chopperDetectorDistance, polarizationConfig + polarizationConfig = meta.polarizationConfig 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 <-|') + logging.info(f' polarization config: {polarizationConfig}') 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 + tofOffset = chopperTriggerPhase * tau / 180 # mismatch of chopper pulse and time-zero 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 @@ -896,8 +880,8 @@ def process(dataPath, ident, clas): 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('chopper speed={:>4.0f} rpm, phase={:>5.3f} deg, tau={} s'.format(30/tau, chopperTriggerPhase, tau)) + logging.info('nr={}, spn={}, cnt={}, tme={}'.format(number, polarizationConfig, 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() @@ -912,7 +896,7 @@ def process(dataPath, ident, clas): #string = f"plott.{plotType} (numberString, '{arg}', data_e)" try: plotFunction = getattr(plott, plotType) - plotFunction(numberString, arg, data_e) + plotFunction(numberString, arg, data_e, polarizationConfig) plt.close() except Exception as e: logging.error(f"ERROR: '{plotType}' is no known output format!") @@ -932,6 +916,8 @@ def commandLineArgs(): help ="chopper speed in rpm") clas.add_argument("-d", "--dataPath", help ="relative path to directory with .hdf files") + clas.add_argument("-D", "--absDataPath", + help ="absolute path to directory with .hdf files") clas.add_argument("-f", "--fileIdent", default='0', help ="file number or offset (if negative)") @@ -959,7 +945,7 @@ def commandLineArgs(): default=99., type=float, help ="value of nu") - clas.add_argument("-P", "--chopperPhase", + clas.add_argument("-P", "--chopperTriggerPhase", default=-7.5, type=float, help ="chopper phase offset") @@ -1007,12 +993,14 @@ def get_dataPath(clas): dataPath = clas.dataPath + '/' if not os.path.exists(dataPath): sys.exit('# *** the directory "'+dataPath+'" does not exist ***') + if clas.absDataPath: + dataPath = clas.absDataPath + '/' 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 , default is "./raw")') + sys.exit('# *** please provide the path to the .hdf data files (-d | -D , default is "./raw")') return dataPath diff --git a/events2histogram_2025.py b/events2histogram_2025.py deleted file mode 100755 index 5a97802..0000000 --- a/events2histogram_2025.py +++ /dev/null @@ -1,1051 +0,0 @@ -__version__ = '2025-06-07' - -# essential changes with regard to 2024 version -# - accepts new hdf structure -# TODO: -# - data path for output -# - solve confusion between negative file number and ranges (check for '-') - -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 -from orsopy import fileio - -#============================================================================== -#============================================================================== -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 - # t_cor = (delta / 180 deg) * tau - 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 + np.rad2deg( np.arctan2(data_e[:,5], detectorDistance) ) - 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 - if mu > -0.25: - data_e[:,9] = 4*np.pi * np.sin( np.deg2rad( data_e[:,8] ) ) / data_e[:,7] - else: - 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/instrument_control_parameters/mu'], 0)) - self.nu = float(np.take(fh['/entry1/Amor/instrument_control_parameters/nu'], 0)) - self.kap = float(np.take(fh['/entry1/Amor/instrument_control_parameters/kap'], 0)) - self.kad = float(np.take(fh['/entry1/Amor/instrument_control_parameters/kad'], 0)) - self.div = float(np.take(fh['/entry1/Amor/instrument_control_parameters/div'], 0)) - chopperSpeed = float(np.take(fh['/entry1/Amor/chopper/rotation_speed'], 0)) - chopperPhase = float(np.take(fh['/entry1/Amor/chopper/phase'], 0)) - ch1TriggerPhase = float(np.take(fh['/entry1/Amor/chopper/ch1_trigger_phase'], 0)) - polarizationConfigLabel = float(np.take(fh['/entry1/Amor/polarization/configuration/value'], 0)) - except (KeyError, IndexError): - logging.warning(f" using parameters from nicos cache") - cachePath = '/home/amor/cache/' - 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] - chopperSpeed = float(value) - value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-chopper_speed'+year_date)).split('\t')[-1] - chopperPhase = float(value) - value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-chopper_phase'+year_date)).split('\t')[-1] - ch1TriggerPhase = float(value) - value = str(subprocess.getoutput('/usr/bin/grep "value" '+cachePath+'nicos-ch1_trigger_phase'+year_date)).split('\t')[-1] - - self.tau = 30. / chopperSpeed - - self.chopperTriggerPhase = ch1TriggerPhase + chopperPhase/2 - - polarizationConfigs = ['undefined', 'oo', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm'] - self.polarizationConfig = polarizationConfigs[int(polarizationConfigLabel)] - - # for .ort header - - self.title = str(fh['entry1/title'][0].decode('utf-8')) - - self.proposal_id = str(fh['entry1/proposal_id'][0].decode('utf-8')) - - self.userName = str(fh['entry1/user/name'][0].decode('utf-8')) - try: - self.userEmail = str(fh['entry1/user/email'][0].decode('utf-8')) - except: - self.userEmail = None - - self.sampleName = str(fh['entry1/sample/name'][0].decode('utf-8')) - self.sampleModel = str(fh['entry1/sample/model'][0].decode('utf-8')) - - #!!! self.instrumentName = str(fh['entry1/Amor/name'][0]) - self.instrumentName = 'Amor' - #!!! self.instrumentType = str(fh['entry1/Amor/type'][0]) - self.source = str(fh['entry1/Amor/source/name'][0]) - #!!! self.sourceType = str(fh['entry1/Amor/source/type'][0]) - #!!! self.sourceProbe = str(fh['entry1/Amor/source/probe'][0]) - self.sourceProbe = 'neutron' - - 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)=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, - -#============================================================================== -def ort_header(fileName): - - #print(fileName) - meta = Meta(fileName) - ''' - Build information object for ORSO file headers. - ''' - - #for fidx in self.get_flist(self.Files): - # fdata = data_vault['tmp%s'%fidx] - # fdate = datetime.strptime(fdata['timestamp'], '%Y-%m-%d %H:%M:%S') - # datafiles.append(fileio.File(file = fdata['srcname'], timestamp = fdate)) - - inst = fileio.InstrumentSettings( - incident_angle = fileio.ValueRange( - mu + meta.kap + meta.kad - meta.div/2, - mu + meta.kap + meta.kad + meta.div/2, - 'deg' - ), - wavelength = fileio.ValueRange( - lamdaMin, - lamdaMax, - 'angstrom' - ), - ) - inst.mu = fileio.Value(mu, 'deg', comment = 'sample angle to horizon') - inst.nu = fileio.Value(nu, 'deg', comment = 'detector angle to horizon') - inst.kap = fileio.Value(meta.kap, 'deg', comment = 'nominal beam inclination') - inst.kad = fileio.Value(meta.kad, 'deg', comment = 'offset of beam inclination') - inst.div = fileio.Value(meta.div, 'deg', comment = 'incoming beam divergence') - mess = fileio.Measurement( - instrument_settings = inst, - data_files = fileName.split('/')[-1], - counting_time = fileio.Value(meta.countingTime, 's'), - scheme = 'angle- and energy-dispersive') - se_dict = None - smpl = fileio.Sample( - name = meta.sampleName, - model = meta.sampleModel, - sample_parameters = se_dict - ) - experiment = fileio.Experiment( - title = meta.title, - instrument = meta.instrumentName, - #instrument_type = meta.instrumentType, - start_date = meta.startDate, - probe = meta.sourceProbe, - #facility = meta.sourceName, - #source_type = meta.sourceType, - proposalID = meta.proposal_id, - ) - ds = fileio.DataSource( - fileio.Person( - meta.userName, - None, - contact = meta.userEmail - ), - experiment, - smpl, - mess - ) - red = fileio.Reduction( - software = fileio.Software( - 'events2histogram', - version = __version__ - ), - timestamp = datetime.now(), - creator = None, - corrections = None, - computer = None, - call = 'python '+' '.join(sys.argv) - ) - cols = [ - fileio.Column('Qz', '1/angstrom'), - fileio.Column('R', comment = 'uncorrected intensity'), - fileio.Column('sR', comment = 'sigma of gaussian probability function'), - fileio.Column('sQz', '1/angstrom', comment = 'resolution based only on sigma_lambda!') - ] - - header = fileio.Orso(ds, red, cols) - - return header - -#============================================================================== -class PlotSelection: - - # header / meta data - - def header(self, filename, mu, nu, totalCounts, countingTime, polarizationConfig): - number = filename.split('n')[1].split('.')[0].lstrip('0') - header = "#{} \u03bc={:>1.2f} \u03bd={:>1.2f} {} {:>10} cts {:>8.1f} s".format(number, mu+5e-3, nu+5e-3, polarizationConfig, totalCounts, countingTime) - return header - - def headline(self, numberString, totalCounts, polarizationConfig): - headLine = "#{} \u03bc={:>1.2f} \u03bd={:>1.2f} p={} {:>12,} cts {:>8.1f} s".format(numberString, mu+5e-3, nu+5e-3, polarizationConfig, totalCounts, countingTime) - return headLine - - # grids - - def lamda_grid(self): - dldl = 0.005 # Delta lambda / lambda - if foldback: - lamda_grid = lamdaMin*(1+dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+dldl)+1)) - else: - lamda_grid = np.arange(0.01, 2.*lamdaMax-2.*lamdaMin, 0.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=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, polarizationConfig): - #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 )]) - q_lim = np.array([qMin, qMax]) - 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], polarizationConfig) - 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) - #plt.close() - - # create PNG with one plot - - def Iyz(self, numberString, arg, data_e, polarizationConfig): - 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 == 'file': - print('# y z conts') - for y in range(len(bins_y)-1): - for z in range(len(bins_z)-1): - print(" %6.3f %6.4f %10.3e" %(bins_y[y], bins_z[z], I_yz[y,z])) - print("") - return - elif arg == 'log': - vmin = 0 - vmax = max(1, np.log(np.max(I_yz)+.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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.colorbar() - plt.savefig(output, format='png', dpi=150) - #plt.close() - - def Ilt(self, numberString, arg, data_e, polarizationConfig) : - 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 == 'file': - print('# lambda theta conts') - for l in range(len(bins_l)-1): - for t in range(len(bins_t)-1): - print(" %6.3f %6.4f %10.3e" %(bins_l[l], bins_t[t], I_lt[l,t]/ bins_t[t])) - print("") - return - elif 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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.colorbar() - plt.savefig(output, format='png', dpi=150) - #plt.close() - - def Itz(self, numberString, arg, data_e, polarizationConfig): - 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]) - if foldback: - time_grid = np.arange(0, tau, 0.0005) - else: - time_grid = np.arange(0, 2.*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 == 'file': - print('# time z conts') - for t in range(len(bins_t)-1): - for z in range(len(bins_z)-1): - print(" %6.3f %6.4f %10.3e" %(bins_t[t], bins_z[z], I_tz[t,z])) - print("") - return - elif 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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.colorbar() - plt.savefig(output, format='png', dpi=150) - #plt.close() - - def Iq(self, numberString, arg, data_e, polarizationConfig): - 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 )]) - q_lim = [qMin, qMax] - print(q_lim) - if arg == 'file': - header = '# q counts' - I_q = np.vstack((bins_q[:-1], I_q, err_q)) - np.savetxt(sys.stdout, I_q.T, header=header) - logging.info(' use `-p ort` instead of `-p Iq`!') - return - elif 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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.savefig(output, format='png', dpi=150) - #plt.close() - - def Il(self, numberString, arg, data_e, polarizationConfig): - I_l, bins_l = np.histogram(data_e[:,7], bins = self.lamda_grid()) - if arg == 'file': - header = '# lambda counts' - I_l = np.vstack((bins_l[:-1], I_l)) - np.savetxt(sys.stdout, I_l.T, header=header) - return - elif 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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.savefig(output, format='png', dpi=150) - #plt.close() - - def It(self, numberString, arg, data_e, polarizationConfig): - I_t, bins_t = np.histogram(data_e[:,8], bins = self.theta_grid()) - if arg == 'file': - header = '# 2theta counts' - I_t = np.vstack((bins_t[:-1], I_t)) - np.savetxt(sys.stdout, I_t.T, header=header) - return - else: - 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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.savefig(output, format='png', dpi=150) - #plt.close() - - def tof(self, numberString, arg, data_e, polarizationConfig): - if foldback: - time_grid = np.arange(0, 1.3*tau, 0.0005) - else: - time_grid = np.arange(0, 2.*tau, 0.0005) - I_t, bins_t = np.histogram(data_e[:,0], bins = time_grid) - if arg == 'file': - header = '# time counts' - I_t = np.vstack((bins_t[:-1], I_t)) - np.savetxt(sys.stdout, I_t.T, header=header) - return - elif 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], polarizationConfig) - plt.title(headline, loc='left', y=1.0, c='r') - plt.savefig(output, format='png') - - def ort(self, numberString, arg, data_e): - I_q, bins_q = np.histogram(data_e[:,9], bins = self.q_grid()) - sI_q = np.sqrt(I_q) - sq_q = bins_q[:-1]*0.022 - I_q = np.vstack((bins_q[:-1], I_q, sI_q, sq_q)) - - datasets = [] - fileNumber = get_flist(numberString)[-1] - fileName = fileNameCreator(dataPath, fileNumber)[0] - header = ort_header(fileName) - orso_data=fileio.OrsoDataset(header, I_q.T) - datasets.append(orso_data) - fileio.save_orso(datasets, 'e2h.ort', data_separator='\n') - -#============================================================================== -#============================================================================== -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] - global foldback - foldback = not clas.noTOFCorrection - - #================================ - # 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: - chopperTriggerPhase - except NameError: - chopperTriggerPhase = meta.chopperTriggerPhase - - - global countingTime, detectorDistance, chopperDetectorDistance, polarizationConfig - polarizationConfig = meta.polarizationConfig - detectorDistance = meta.detectorDistance - chopperDetectorDistance = meta.chopperDetectorDistance - countingTime = meta.countingTime - - if verbous: - logging.info(" mu = {:>4.2f} deg, nu = {:>4.2f} deg".format(mu, nu)) - logging.info(f' polarization config: {polarizationConfig}') - - 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) - - if foldback: - tof_e = np.where(tof_etau+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, chopperTriggerPhase, tau)) - logging.info('nr={}, spn={}, cnt={}, tme={}'.format(number, polarizationConfig, 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, polarizationConfig) - 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("-b", "--noTOFCorrection", - action='store_true', - help ="do not correct tof of seond neutron pulse") - 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("-D", "--absDataPath", - help ="absolute 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=-7.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("-r", "--iDonno", - action='store_true', - help ="no idea") - clas.add_argument("-s", "--spy", - action='store_true', - help ="report a few key values, no plotting or writing") - 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("-u", "--update", - action='store_true', - help ="update output every 5 seconds") - 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 ***') - if clas.absDataPath: - dataPath = clas.absDataPath + '/' - 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 | -D , 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() - - if clas.update: - verbous = False - logging.basicConfig(level=logging.ERROR, format='# %(message)s') - logging.error('e2h: recreating "e2h_life.png" every 5 s') - delay = 5 - output = 'e2h_life.png' - oldtime = 0 - while True: - dataPath = get_directDataPath(clas) - fileName = fileNameCreator(dataPath, 0)[0] - newtime = os.path.getmtime(fileName) - if newtime-oldtime: - print('\r# processing (press ^C to stop)', end='', flush=True) - process(dataPath, '0', clas) - oldtime = newtime - for i in range(5): - print('\r# waiting'+'.'*i+' '*(5-i)+' (press ^C to stop)', end='', flush=True) - time.sleep(delay/5) - signal.signal(signal.SIGINT, endIt) - else: - 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() - - diff --git a/libeos/command_line.py b/libeos/command_line.py deleted file mode 100644 index f77d6c6..0000000 --- a/libeos/command_line.py +++ /dev/null @@ -1,220 +0,0 @@ -import argparse - -from .logconfig import update_loglevel -from .options import ReaderConfig, EOSConfig, ExperimentConfig, OutputConfig, ReductionConfig, Defaults - - -def commandLineArgs(): - """ - Process command line argument. - The type of the default values is used for conversion and validation. - """ - msg = "eos reads data from (one or several) raw file(s) of the .hdf format, \ - performs various corrections, conversations and projections and exports\ - the resulting reflectivity in an orso-compatible format." - clas = argparse.ArgumentParser(description = msg) - - input_data = clas.add_argument_group('input data') - input_data.add_argument("-f", "--fileIdentifier", - required = True, - nargs = '+', - help = "file number(s) or offset (if < 1)") - input_data.add_argument("-n", "--normalisationFileIdentifier", - default = Defaults.normalisationFileIdentifier, - nargs = '+', - help = "file number(s) of normalisation measurement") - input_data.add_argument("-rp", "--rawPath", - type = str, - default = Defaults.rawPath, - help = "ath to directory with .hdf files") - input_data.add_argument("-Y", "--year", - default = Defaults.year, - type = int, - help = "year the measurement was performed") - input_data.add_argument("-sub", "--subtract", - help = "R(q_z) curve to be subtracted (in .Rqz.ort format)") - input_data.add_argument("-nm", "--normalisationMethod", - default = Defaults.normalisationMethod, - help = "normalisation method: [o]verillumination, [u]nderillumination, [d]irect_beam") - input_data.add_argument("-mt", "--monitorType", - type = str, - default = Defaults.monitorType, - help = "one of [p]rotonCurrent, [t]ime or [n]eutronMonitor") - - output = clas.add_argument_group('output') - output.add_argument("-o", "--outputName", - default = Defaults.outputName, - help = "output file name (withot suffix)") - output.add_argument("-op", "--outputPath", - type = str, - default = Defaults.outputPath, - help = "path for output") - output.add_argument("-of", "--outputFormat", - nargs = '+', - default = Defaults.outputFormat, - help = "one of [Rqz.ort, Rlt.ort]") - output.add_argument("-ai", "--incidentAngle", - type = str, - default = Defaults.incidentAngle, - help = "calulate alpha_i from [alphaF, mu, nu]", - ) - output.add_argument("-r", "--qResolution", - default = Defaults.qResolution, - type = float, - help = "q_z resolution") - output.add_argument("-ts", "--timeSlize", - nargs = '+', - type = float, - help = "time slizing ,[ [,stop]]") - output.add_argument("-s", "--scale", - nargs = '+', - default = Defaults.scale, - type = float, - help = "scaling factor for R(q_z)") - output.add_argument("-S", "--autoscale", - nargs = 2, - type = float, - help = "scale to 1 in the given q_z range") - - masks = clas.add_argument_group('masks') - masks.add_argument("-l", "--lambdaRange", - default = Defaults.lambdaRange, - nargs = 2, - type = float, - help = "wavelength range") - masks.add_argument("-t", "--thetaRange", - default = Defaults.thetaRange, - nargs = 2, - type = float, - help = "absolute theta range") - masks.add_argument("-T", "--thetaRangeR", - default = Defaults.thetaRangeR, - nargs = 2, - type = float, - help = "relative theta range") - masks.add_argument("-y", "--yRange", - default = Defaults.yRange, - nargs = 2, - type = int, - help = "detector y range") - masks.add_argument("-q", "--qzRange", - default = Defaults.qzRange, - nargs = 2, - type = float, - help = "q_z range") - masks.add_argument("-ct", "--lowCurrentThreshold", - default = Defaults.lowCurrentThreshold, - type = float, - help = "proton current threshold for discarding neutron pulses") - - - overwrite = clas.add_argument_group('overwrite') - overwrite.add_argument("-cs", "--chopperSpeed", - default = Defaults.chopperSpeed, - type = float, - help = "chopper speed in rpm") - overwrite.add_argument("-cp", "--chopperPhase", - default = Defaults.chopperPhase, - type = float, - help = "chopper phase") - overwrite.add_argument("-co", "--chopperPhaseOffset", - default = Defaults.chopperPhaseOffset, - type = float, - help = "phase offset between chopper opening and trigger pulse") - overwrite.add_argument("-m", "--muOffset", - default = Defaults.muOffset, - type = float, - help = "mu offset") - overwrite.add_argument("-mu", "--mu", - default = Defaults.mu, - type = float, - help ="value of mu") - overwrite.add_argument("-nu", "--nu", - default = Defaults.nu, - type = float, - help = "value of nu") - overwrite.add_argument("-sm", "--sampleModel", - default = Defaults.sampleModel, - type = str, - help = "1-line orso sample model description") - - misc = clas.add_argument_group('misc') - misc.add_argument('-v', '--verbose', action='store_true') - misc.add_argument('-vv', '--debug', action='store_true') - - return clas.parse_args() - - -def expand_file_list(short_notation): - """Evaluate string entry for file number lists""" - #log().debug('Executing get_flist') - file_list=[] - for i in short_notation.split(','): - if '-' in i: - if ':' in i: - step = i.split(':', 1)[1] - file_list += range(int(i.split('-', 1)[0]), int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, int(step)) - else: - step = 1 - file_list += range(int(i.split('-', 1)[0]), int(i.split('-', 1)[1])+1, int(step)) - else: - file_list += [int(i)] - - return sorted(file_list) - - -def output_format_list(outputFormat): - format_list = [] - if 'ort' in outputFormat or 'Rqz.ort' in outputFormat or 'Rqz' in outputFormat: - format_list.append('Rqz.ort') - if 'ort' in outputFormat or 'Rlt.ort' in outputFormat or 'Rlt' in outputFormat: - format_list.append('Rlt.ort') - if 'orb' in outputFormat or 'Rqz.orb' in outputFormat or 'Rqz' in outputFormat: - format_list.append('Rqz.orb') - if 'orb' in outputFormat or 'Rlt.orb' in outputFormat or 'Rlt' in outputFormat: - format_list.append('Rlt.orb') - return sorted(format_list, reverse=True) - -def command_line_options(): - clas = commandLineArgs() - update_loglevel(clas.verbose, clas.debug) - - reader_config = ReaderConfig( - year = clas.year, - rawPath = clas.rawPath, - ) - experiment_config = ExperimentConfig( - sampleModel = clas.sampleModel, - chopperSpeed = clas.chopperSpeed, - chopperPhase = clas.chopperPhase, - chopperPhaseOffset = clas.chopperPhaseOffset, - yRange = clas.yRange, - lambdaRange = clas.lambdaRange, - qzRange = clas.qzRange, - lowCurrentThreshold = clas.lowCurrentThreshold, - incidentAngle = clas.incidentAngle, - mu = clas.mu, - nu = clas.nu, - muOffset = clas.muOffset, - monitorType = clas.monitorType, - ) - reduction_config = ReductionConfig( - qResolution = clas.qResolution, - qzRange = clas.qzRange, - autoscale = clas.autoscale, - thetaRange = clas.thetaRange, - thetaRangeR = clas.thetaRangeR, - fileIdentifier = clas.fileIdentifier, - scale = clas.scale, - subtract = clas.subtract, - normalisationFileIdentifier = clas.normalisationFileIdentifier, - normalisationMethod = clas.normalisationMethod, - timeSlize = clas.timeSlize, - ) - output_config = OutputConfig( - outputFormats = output_format_list(clas.outputFormat), - outputName = clas.outputName, - outputPath = clas.outputPath, - ) - - return EOSConfig(reader_config, experiment_config, reduction_config, output_config) diff --git a/libeos/file_reader.py b/libeos/file_reader.py deleted file mode 100644 index e96349f..0000000 --- a/libeos/file_reader.py +++ /dev/null @@ -1,503 +0,0 @@ -import logging -import os -import subprocess -import sys -from datetime import datetime, timezone -try: - import zoneinfo -except ImportError: - # for python versions < 3.9 try to use the backports version - from backports import zoneinfo -from typing import List - -import h5py -import numpy as np -from orsopy import fileio -from orsopy.fileio.model_language import SampleModel - -from . import const -from .header import Header -from .instrument import Detector -from .options import ExperimentConfig, ReaderConfig - -try: - from . import nb_helpers -except Exception: - nb_helpers = None - -# Time zone used to interpret time strings -AMOR_LOCAL_TIMEZONE = zoneinfo.ZoneInfo(key='Europe/Zurich') - -class AmorData: - """read meta-data and event streams from .hdf file(s), apply filters and conversions""" - chopperDetectorDistance: float - chopperDistance: float - chopperPhase: float - chopperSpeed: float - chopper1TriggerPhase: float - chopper2TriggerPhase: float - div: float - data_file_numbers: List[int] - delta_z: np.ndarray - detZ_e: np.ndarray - lamda_e: np.ndarray - wallTime_e: np.ndarray - kad: float - kap: float - lambdaMax: float - lambda_e: np.ndarray - #monitor: float - mu: float - nu: float - tau: float - tofCut: float - start_date: str - monitorType: str - - seriesStartTime = None - - #------------------------------------------------------------------------------------------------- - def __init__(self, header: Header, reader_config: ReaderConfig, config: ExperimentConfig, - short_notation:str, norm=False): - #self.startTime = reader_config.startTime - self.header = header - self.config = config - self.reader_config = reader_config - self.expand_file_list(short_notation) - self.read_data(norm=norm) - - #------------------------------------------------------------------------------------------------- - def read_data(self, norm=False): - self.file_list = [] - for number in self.data_file_numbers: - self.file_list.append(self.path_generator(number)) - ## read specific meta data and measurement from first file - if norm: - self.readHeaderInfo = False - else: - self.readHeaderInfo = True - - _detZ_e = [] - _lamda_e = [] - _wallTime_e = [] - #_monitor = 0 - _monitorPerPulse = [] - _pulseTimeS = [] - for file in self.file_list: - self.read_individual_data(file, norm) - _detZ_e = np.append(_detZ_e, self.detZ_e) - _lamda_e = np.append(_lamda_e, self.lamda_e) - _wallTime_e = np.append(_wallTime_e, self.wallTime_e) - _monitorPerPulse = np.append(_monitorPerPulse, self.monitorPerPulse) - _pulseTimeS = np.append(_pulseTimeS, self.pulseTimeS) - #_monitor += self.monitor - self.detZ_e = _detZ_e - self.lamda_e = _lamda_e - self.wallTime_e = _wallTime_e - #self.monitor = _monitor - self.monitorPerPulse = _monitorPerPulse - self.pulseTimeS = _pulseTimeS - - #------------------------------------------------------------------------------------------------- - def path_generator(self, number): - fileName = f'amor{self.reader_config.year}n{number:06d}.hdf' - path = '' - for rawd in self.reader_config.rawPath: - if os.path.exists(os.path.join(rawd,fileName)): - path = rawd - break - if not path: - if os.path.exists(f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}/{fileName}'): - path = f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}' - else: - sys.exit(f'# ERROR: the file {fileName} can not be found in {self.reader_config.rawPath}') - return os.path.join(path, fileName) - #------------------------------------------------------------------------------------------------- - def expand_file_list(self, short_notation): - """Evaluate string entry for file number lists""" - #log().debug('Executing get_flist') - file_list=[] - for i in short_notation.split(','): - if '-' in i: - if ':' in i: - step = i.split(':', 1)[1] - file_list += range(int(i.split('-', 1)[0]), int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, int(step)) - else: - step = 1 - file_list += range(int(i.split('-', 1)[0]), int(i.split('-', 1)[1])+1, int(step)) - else: - file_list += [int(i)] - self.data_file_numbers=sorted(file_list) - #------------------------------------------------------------------------------------------------- - def resolve_pixels(self): - """determine spatial coordinats and angles from pixel number""" - nPixel = Detector.nWires * Detector.nStripes * Detector.nBlades - pixelID = np.arange(nPixel) - (bladeNr, bPixel) = np.divmod(pixelID, Detector.nWires * Detector.nStripes) - (bZi, detYi) = np.divmod(bPixel, Detector.nStripes) # z index on blade, y index on detector - detZi = bladeNr * Detector.nWires + bZi # z index on detector - detX = bZi * Detector.dX # x position in detector - # detZ = Detector.zero - bladeNr * Detector.bladeZ - bZi * Detector.dZ # z position on detector - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*Detector.bladeZ / Detector.distance) ) - delta = (Detector.nBlades/2. - bladeNr) * bladeAngle \ - - np.rad2deg( np.arctan(bZi*Detector.dZ / ( Detector.distance + bZi * Detector.dX) ) ) - self.delta_z = delta[detYi==1] - return np.vstack((detYi.T, detZi.T, detX.T, delta.T)).T - #------------------------------------------------------------------------------------------------- - def read_individual_data(self, fileName, norm=False): - self.hdf = h5py.File(fileName, 'r', swmr=True) - - if self.readHeaderInfo: - self.read_header_info() - - logging.warning(f' from file: {fileName}') - self.read_individual_header() - - # add header content - if self.readHeaderInfo: - self.readHeaderInfo = False - self.header.measurement_instrument_settings = fileio.InstrumentSettings( - incident_angle = fileio.ValueRange(round(self.mu+self.kap+self.kad-0.5*self.div, 3), - round(self.mu+self.kap+self.kad+0.5*self.div, 3), - 'deg'), - wavelength = fileio.ValueRange(const.lamdaCut, self.config.lambdaRange[1], 'angstrom'), - #polarization = fileio.Polarization.unpolarized, - polarization = self.polarizationConfig - ) - self.header.measurement_instrument_settings.mu = fileio.Value(round(self.mu, 3), 'deg', comment='sample angle to horizon') - self.header.measurement_instrument_settings.nu = fileio.Value(round(self.nu, 3), 'deg', comment='detector angle to horizon') - self.header.measurement_instrument_settings.div = fileio.Value(round(self.div, 3), 'deg', comment='incoming beam divergence') - self.header.measurement_instrument_settings.kap = fileio.Value(round(self.kap, 3), 'deg', comment='incoming beam inclination') - if abs(self.kad)>0.02: - self.header.measurement_instrument_settings.kad = fileio.Value(round(self.kad, 3), 'deg', comment='incoming beam angular offset') - if norm: - self.header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=self.fileDate)) - else: - self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=self.fileDate)) - logging.info(f' mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kad:6.3f}') - - self.read_event_stream() - - self.correct_for_chopper_phases() - - self.read_chopper_trigger_stream() - - self.extract_walltime(norm) - - self.read_proton_current_stream() - - self.associate_pulse_with_monitor() - - # following lines: debugging output to trace the time-offset of proton current and neutron pulses - if self.config.monitorType == 'x': - cpp, t_bins = np.histogram(self.wallTime_e, self.pulseTimeS) - np.savetxt('tme.hst', np.vstack((self.pulseTimeS[:-1], cpp, self.monitorPerPulse[:-1])).T) - - #self.average_events_per_pulse() # for debugging only. VERY time consuming!!! - - self.monitor_threshold() - - self.filter_strange_times() - - self.merge_time_frames() - - self.filter_project_x() - - self.correct_for_chopper_opening() - - self.calculate_derived_properties() - - self.filter_qz_range(norm) - - logging.info(f' number of events: total = {self.totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') - - def read_event_stream(self): - self.tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 - self.pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=np.int64) - self.dataPacket_p = np.array(self.hdf['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) - self.dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.int64) - - def correct_for_chopper_phases(self): - #print(f'chopperTriggerPhase: {self.ch1TriggerPhase}') - self.tof_e += self.tau * (self.ch1TriggerPhase + self.chopperPhase/2)/180 - - def extract_walltime(self, norm): - if nb_helpers: - self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) - else: - self.wallTime_e = np.empty(np.shape(self.tof_e)[0], dtype=np.int64) - for i in range(len(self.dataPacket_p)-1): - self.wallTime_e[self.dataPacket_p[i]:self.dataPacket_p[i+1]] = self.dataPacketTime_p[i] - self.wallTime_e[self.dataPacket_p[-1]:] = self.dataPacketTime_p[-1] - self.wallTime_e -= np.int64(self.seriesStartTime) - logging.debug(f' wall time from {self.wallTime_e[0]/1e9:6.1f} s to {self.wallTime_e[-1]/1e9:6.1f} s') - - def read_chopper_trigger_stream(self): - self.chopper1TriggerTime = np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][:-2], dtype=np.int64) - #self.chopper2TriggerTime = self.chopper1TriggerTime + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time'][:-2], dtype=np.int64) - # + np.array(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][:], dtype=np.int64) - if np.shape(self.chopper1TriggerTime)[0] > 2: - self.startTime = self.chopper1TriggerTime[0] - self.stopTime = self.chopper1TriggerTime[-1] - self.pulseTimeS = self.chopper1TriggerTime - else: - logging.warn(' no chopper trigger data available, using event steram instead') - self.startTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][0], dtype=np.int64) - self.stopTime = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-2], dtype=np.int64) - self.pulseTimeS = np.arange(self.startTime, self.stopTime, self.tau*1e9) - if self.seriesStartTime is None: - self.seriesStartTime = self.startTime - logging.debug(f' series start time (epoch): {self.seriesStartTime/1e9:13.2f} s') - self.pulseTimeS -= self.seriesStartTime - logging.debug(f' epoch time from {self.startTime/1e9:13.2f} s to {self.stopTime/1e9:13.2f} s') - logging.debug(f' => counting time {self.stopTime/1e9-self.startTime/1e9:8.2f} s') - - def read_proton_current_stream(self): - self.currentTime = np.array(self.hdf['entry1/Amor/detector/proton_current/time'][:], dtype=np.int64) - self.current = np.array(self.hdf['entry1/Amor/detector/proton_current/value'][:,0], dtype=float) - - def get_current_per_pulse(self, pulseTimeS, currentTimeS, currents): - # add currents for early pulses and current time value after last pulse (j+1) - currentTimeS = np.hstack([[0], currentTimeS, [pulseTimeS[-1]+1]]) - currents = np.hstack([[0], currents]) - pulseCurrentS = np.zeros(pulseTimeS.shape[0], dtype=float) - j = 0 - for i, ti in enumerate(pulseTimeS): - while ti >= currentTimeS[j+1]: - j += 1 - pulseCurrentS[i] = currents[j] - #print(f' {i} {pulseTimeS[i]} {pulseCurrentS[i]}') - return pulseCurrentS - - def associate_pulse_with_monitor(self): - if self.config.monitorType == 'p': # protonCharge - self.currentTime -= np.int64(self.seriesStartTime) - self.monitorPerPulse = self.get_current_per_pulse(self.pulseTimeS, self.currentTime, self.current) * 2*self.tau * 1e-3 - # filter low-current pulses - self.monitorPerPulse = np.where(self.monitorPerPulse > 2*self.tau * self.config.lowCurrentThreshold * 1e-3, self.monitorPerPulse, 0) - elif self.config.monitorType == 't': # countingTime - self.monitorPerPulse = np.ones(np.shape(self.pulseTimeS)[0])*2*self.tau - else: # pulses - self.monitorPerPulse = np.ones(np.shape(self.pulseTimeS)[0]) - - def average_events_per_pulse(self): - if self.config.monitorType == 'p': - for i, time in enumerate(self.pulseTimeS): - events = np.shape(self.wallTime_e[self.wallTime_e == time])[0] - logging.info(f'pulse: {i:6.0f}, events: {events:6.0f}, monitor: {self.monitorPerPulse[i]:6.2f}') - - def monitor_threshold(self): - #if self.config.monitorType == 'p': # fix to check for file compatibility - self.totalNumber = np.shape(self.tof_e[self.tof_e<=self.stopTime])[0] - if True: - goodTimeS = self.pulseTimeS[self.monitorPerPulse!=0] - filter_e = np.where(np.isin(self.wallTime_e, goodTimeS), True, False) - self.tof_e = self.tof_e[filter_e] - self.pixelID_e = self.pixelID_e[filter_e] - self.wallTime_e = self.wallTime_e[filter_e] - logging.info(f' low-beam rejected pulses: {np.shape(self.monitorPerPulse)[0]-1-np.shape(goodTimeS)[0]} out of {np.shape(self.monitorPerPulse)[0]-1}') - logging.info(f' with {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]} events') - logging.info(f' average counts per pulse = {np.shape(self.tof_e)[0] / np.shape(goodTimeS[goodTimeS!=0])[0]:7.1f}') - - def filter_qz_range(self, norm): - if self.config.qzRange[1]<0.5 and not norm: - self.mask_e = np.logical_and(self.mask_e, - (self.config.qzRange[0]<=self.qz_e) & (self.qz_e<=self.config.qzRange[1])) - self.detZ_e = self.detZ_e[self.mask_e] - self.lamda_e = self.lamda_e[self.mask_e] - self.wallTime_e = self.wallTime_e[self.mask_e] - - def calculate_derived_properties(self): - self.lamdaMax = const.lamdaCut+1.e13*self.tau*const.hdm/(self.chopperDetectorDistance+124.) - #if nb_helpers: - if False: - self.lamda_e, self.qz_e, self.mask_e = nb_helpers.calculate_derived_properties_focussing( - self.tof_e, self.detXdist_e, self.delta_e, self.mask_e, - self.config.lambdaRange[0], self.config.lambdaRange[1], self.nu, self.mu, - self.chopperDetectorDistance, const.hdm - ) - return - # lambda - self.lamda_e = (1.e13*const.hdm)*self.tof_e/(self.chopperDetectorDistance+self.detXdist_e) - self.mask_e = np.logical_and(self.mask_e, (self.config.lambdaRange[0]<=self.lamda_e) & ( - self.lamda_e<=self.config.lambdaRange[1])) - # alpha_f - # q_z - if self.config.incidentAngle == 'alphaF': - alphaF_e = self.nu - self.mu + self.delta_e - self.qz_e = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/self.lamda_e) - # qx_e = 0. - self.header.measurement_scheme = 'angle- and energy-dispersive' - elif self.config.incidentAngle == 'nu': - alphaF_e = (self.nu + self.delta_e + self.kap + self.kad) / 2. - self.qz_e = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/self.lamda_e) - # qx_e = 0. - self.header.measurement_scheme = 'energy-dispersive' - else: - alphaF_e = self.nu - self.mu + self.delta_e - alphaI = self.kap + self.kad + self.mu - self.qz_e = 2*np.pi * ((np.sin(np.deg2rad(alphaF_e)) + np.sin(np.deg2rad(alphaI)))/self.lamda_e) - self.qx_e = 2*np.pi * ((np.cos(np.deg2rad(alphaF_e)) - np.cos(np.deg2rad(alphaI)))/self.lamda_e) - self.header.measurement_scheme = 'energy-dispersive' - - def correct_for_chopper_opening(self): - # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau - if self.config.incidentAngle == 'alphaF': - self.tof_e -= ( self.delta_e / 180. ) * self.tau - else: - # TODO: check sign of correction - self.tof_e -= ( self.kad / 180. ) * self.tau - - def filter_project_x(self): - pixelLookUp = self.resolve_pixels() - if nb_helpers: - (self.detZ_e, self.detXdist_e, self.delta_e, self.mask_e) = nb_helpers.filter_project_x( - pixelLookUp, self.pixelID_e.astype(np.int64), self.config.yRange[0], self.config.yRange[1] - ) - else: - # resolve pixel ID into y and z indicees, x position and angle - (detY_e, self.detZ_e, self.detXdist_e, self.delta_e) = pixelLookUp[np.int_(self.pixelID_e)-1, :].T - # define mask and filter y range - self.mask_e = (self.config.yRange[0]<=detY_e) & (detY_e<=self.config.yRange[1]) - - # TODO: - handle each neutron pulse individually, - associate with correct monitor also for slow neutrons - def merge_time_frames(self): - total_offset = self.tofCut + self.tau * (self.ch1TriggerPhase + self.chopperPhase/2)/180 - if nb_helpers: - self.tof_e = nb_helpers.merge_frames(self.tof_e, self.tofCut, self.tau, total_offset) - else: - self.tof_e = np.remainder(self.tof_e-(self.tofCut-self.tau), self.tau)+total_offset # tof shifted to 1 frame - - def filter_strange_times(self): - # 'strange' tof times are those with t > 2 tau (originating from the efu) - filter_e = (self.tof_e<=2*self.tau) - self.tof_e = self.tof_e[filter_e] - self.pixelID_e = self.pixelID_e[filter_e] - self.wallTime_e = self.wallTime_e[filter_e] - if np.shape(filter_e)[0]-np.shape(self.tof_e)[0]>0.5: - logging.warning(f' strange times: {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]}') - - def read_individual_header(self): - self.chopperDistance = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0)) - self.detectorDistance = float(np.take(self.hdf['entry1/Amor/detector/transformation/distance'], 0)) - self.chopperDetectorDistance = self.detectorDistance-float(np.take(self.hdf['entry1/Amor/chopper/distance'], 0)) - self.tofCut = const.lamdaCut*self.chopperDetectorDistance/const.hdm*1.e-13 - - polarizationConfigs = ['undefined', 'unpolarized', 'po', 'mo', 'op', 'pp', 'mp', 'om', 'pm', 'mm'] - try: - self.mu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/mu'], 0)) - self.nu = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/nu'], 0)) - self.kap = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kap'], 0)) - self.kad = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/kad'], 0)) - self.div = float(np.take(self.hdf['/entry1/Amor/instrument_control_parameters/div'], 0)) - self.ch1TriggerPhase = float(np.take(self.hdf['/entry1/Amor/chopper/ch1_trigger_phase'], 0)) - self.ch2TriggerPhase = float(np.take(self.hdf['/entry1/Amor/chopper/ch2_trigger_phase'], 0)) - try: - chopperTriggerTime = float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][2])\ - - float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_zero'][1]) - self.tau = int(1e-6*chopperTriggerTime/2+0.5)*(1e-3) - self.chopperSpeed = 30/self.tau - chopperTriggerTimeDiff = float(self.hdf['entry1/Amor/chopper/ch2_trigger/event_time_offset'][2]) - chopperTriggerPhase = 180e-9*chopperTriggerTimeDiff/self.tau - #TODO: check the next line - self.chopperPhase = chopperTriggerPhase + self.ch1TriggerPhase - self.ch2TriggerPhase - except(KeyError, IndexError): - logging.debug(' chopper speed and phase taken from .hdf file') - self.chopperSpeed = float(np.take(self.hdf['/entry1/Amor/chopper/rotation_speed'], 0)) - self.chopperPhase = float(np.take(self.hdf['/entry1/Amor/chopper/phase'], 0)) - self.tau = 30/self.chopperSpeed - try: - polarizationConfigLabel = int(self.hdf['/entry1/Amor/polarization/configuration/value'][0]) - except(KeyError, IndexError): - polarizationConfigLabel = 0 - self.polarizationConfig = polarizationConfigs[polarizationConfigLabel] - logging.debug(f' polarization configuration: {self.polarizationConfig} (index {polarizationConfigLabel} (index {polarizationConfigLabel}))') - except(KeyError, IndexError): - logging.warning(" using parameters from nicos cache") - year_date = str(self.start_date).replace('-', '/', 1) - # TODO: check new cache pathes - cachePath = '/home/amor/cache/' - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-mu/{year_date}')).split('\t')[-1] - self.mu = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-nu/{year_date}')).split('\t')[-1] - self.nu = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kap/{year_date}')).split('\t')[-1] - self.kap = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kad/{year_date}')).split('\t')[-1] - self.kad = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-div/{year_date}')).split('\t')[-1] - self.div = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-ch1_speed/{year_date}')).split('\t')[-1] - self.chopperSpeed = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-chopper_phase/{year_date}')).split('\t')[-1] - self.chopperPhase = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-ch1_trigger_phase/{year_date}')).split('\t')[-1] - self.ch1TriggerPhase = float(value) - value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-ch2_trigger_phase/{year_date}')).split('\t')[-1] - self.ch2TriggerPhase = float(value) - - self.tau = 30. / self.chopperSpeed - - logging.debug(f' tau = {self.tau} s') - if self.config.muOffset: - logging.debug(f' set muOffset = {self.config.muOffset}') - self.mu += self.config.muOffset - if self.config.mu: - logging.debug(f' replaced mu = {self.mu} with {self.config.mu}') - self.mu = self.config.mu - if self.config.nu: - logging.debug(f' replaced nu = {self.nu} with {self.config.nu}') - self.nu = self.config.nu - if self.config.chopperPhaseOffset: - logging.debug(f' replaced ch1TriggerPhase = {self.ch1TriggerPhase} with {self.config.chopperPhaseOffset}') - self.ch1TriggerPhase = self.config.chopperPhaseOffset - - # extract start time as unix time, adding UTC offset of 1h to time string - dz = datetime.fromisoformat(self.hdf['/entry1/start_time'][0].decode('utf-8')) - self.fileDate=dz.replace(tzinfo=AMOR_LOCAL_TIMEZONE) - #self.startTime = np.int64( (self.fileDate.timestamp() ) * 1e9 ) - #if self.seriesStartTime is None: - # self.seriesStartTime = self.startTime - - def read_header_info(self): - # read general information and first data set - logging.info(f' meta data from: {self.file_list[0]}') - self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) - title = self.hdf['entry1/title'][0].decode('utf-8') - proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') - user_name = self.hdf['entry1/user/name'][0].decode('utf-8') - user_affiliation = 'unknown' - user_email = self.hdf['entry1/user/email'][0].decode('utf-8') - user_orcid = None - sampleName = self.hdf['entry1/sample/name'][0].decode('utf-8') - model = self.hdf['entry1/sample/model'][0].decode('utf-8') - instrumentName = 'Amor' - source = self.hdf['entry1/Amor/source/name'][0].decode('utf-8') - sourceProbe = 'neutron' - start_time = self.hdf['entry1/start_time'][0].decode('utf-8') - self.start_date = start_time.split(' ')[0] - if self.config.sampleModel: - model = self.config.sampleModel - # assembling orso header information - self.header.owner = fileio.Person( - name=user_name, - affiliation=user_affiliation, - contact=user_email, - ) - if user_orcid: - self.header.owner.orcid = user_orcid - self.header.experiment = fileio.Experiment( - title=title, - instrument=instrumentName, - start_date=self.start_date, - probe=sourceProbe, - facility=source, - proposalID=proposal_id - ) - self.header.sample = fileio.Sample( - name=sampleName, - model=SampleModel(stack=model), - sample_parameters=None, - ) - self.header.measurement_scheme = 'angle- and energy-dispersive' - diff --git a/libeos/options.py b/libeos/options.py deleted file mode 100644 index a936a13..0000000 --- a/libeos/options.py +++ /dev/null @@ -1,193 +0,0 @@ -""" -Classes for stroing various configurations needed for reduction. -""" -from dataclasses import dataclass, field -from typing import Optional, Tuple -from datetime import datetime -from os import path -import numpy as np - -import logging - -class Defaults: - # fileIdentifier - outputPath = '.' - rawPath = ['.', path.join('.','raw'), path.join('..','raw'), path.join('..','..','raw')] - year = datetime.now().year - normalisationFileIdentifier = [] - normalisationMethod = 'o' - monitorType = 'auto' - # subtract - outputName = "fromEOS" - outputFormat = ['Rqz.ort'] - incidentAngle = 'alphaF' - qResolution = 0.01 - #timeSlize - scale = [1] - # autoscale - lambdaRange = [2., 15.] - thetaRange = [-12., 12.] - thetaRangeR = [-0.75, 0.75] - yRange = [11, 41] - qzRange = [0.005, 0.51] - chopperSpeed = 500 - chopperPhase = 0.0 - chopperPhaseOffset = -9.1 - muOffset = 0 - mu = 0 - nu = 0 - sampleModel = None - lowCurrentThreshold = 50 - # - - - -@dataclass -class ReaderConfig: - year: int - rawPath: Tuple[str] - startTime: Optional[float] = 0 - -@dataclass -class ExperimentConfig: - incidentAngle: str - chopperPhase: float - chopperSpeed: float - yRange: Tuple[float, float] - lambdaRange: Tuple[float, float] - qzRange: Tuple[float, float] - monitorType: str - lowCurrentThreshold: float - - sampleModel: Optional[str] = None - chopperPhaseOffset: float = 0 - mu: Optional[float] = None - nu: Optional[float] = None - muOffset: Optional[float] = None - -@dataclass -class ReductionConfig: - normalisationMethod: str - qResolution: float - qzRange: Tuple[float, float] - thetaRange: Tuple[float, float] - thetaRangeR: Tuple[float, float] - - fileIdentifier: list = field(default_factory=lambda: ["0"]) - scale: list = field(default_factory=lambda: [1]) #per file scaling; if less elements than files use the last one - - autoscale: Optional[Tuple[bool, bool]] = None - subtract: Optional[str] = None - normalisationFileIdentifier: Optional[list] = None - timeSlize: Optional[list] = None - -@dataclass -class OutputConfig: - outputFormats: list - outputName: str - outputPath: str - -@dataclass -class EOSConfig: - reader: ReaderConfig - experiment: ExperimentConfig - reduction: ReductionConfig - output: OutputConfig - - _call_string_overwrite=None - - #@property - #def call_string(self)->str: - # if self._call_string_overwrite: - # return self._call_string_overwrite - # else: - # return self.calculate_call_string() - - def call_string(self): - base = 'python eos.py' - - inpt = '' - if self.reader.year: - inpt += f' -Y {self.reader.year}' - else: - inpt += f' -Y {datetime.now().year}' - if np.shape(self.reader.rawPath)[0] == 1: - inpt += f' --rawPath {self.reader.rawPath}' - if self.reduction.subtract: - inpt += f' -subtract {self.reduction.subtract}' - if self.reduction.normalisationFileIdentifier: - inpt += f' -n {" ".join(self.reduction.normalisationFileIdentifier)}' - if self.reduction.fileIdentifier: - inpt += f' -f {" ".join(self.reduction.fileIdentifier)}' - - otpt = '' - if self.reduction.qResolution: - otpt += f' -r {self.reduction.qResolution}' - if self.output.outputPath != '.': - inpt += f' --outputdPath {self.output.outputPath}' - if self.output.outputName: - otpt += f' -o {self.output.outputName}' - if self.output.outputFormats != ['Rqz.ort']: - otpt += f' -of {" ".join(self.output.outputFormats)}' - - mask = '' - if self.experiment.yRange != Defaults.yRange: - mask += f' -y {" ".join(str(ii) for ii in self.experiment.yRange)}' - if self.experiment.lambdaRange!= Defaults.lambdaRange: - mask += f' -l {" ".join(str(ff) for ff in self.experiment.lambdaRange)}' - if self.reduction.thetaRange != Defaults.thetaRange: - mask += f' -t {" ".join(str(ff) for ff in self.reduction.thetaRange)}' - elif self.reduction.thetaRangeR != Defaults.thetaRangeR: - mask += f' -T {" ".join(str(ff) for ff in self.reduction.thetaRangeR)}' - if self.experiment.qzRange!= Defaults.qzRange: - mask += f' -q {" ".join(str(ff) for ff in self.experiment.qzRange)}' - - para = '' - if self.experiment.chopperPhase != Defaults.chopperPhase: - para += f' --chopperPhase {self.experiment.chopperPhase}' - if self.experiment.chopperPhaseOffset != Defaults.chopperPhaseOffset: - para += f' --chopperPhaseOffset {self.experiment.chopperPhaseOffset}' - if self.experiment.mu: - para += f' --mu {self.experiment.mu}' - elif self.experiment.muOffset: - para += f' --muOffset {self.experiment.muOffset}' - if self.experiment.nu: - para += f' --nu {self.experiment.nu}' - - modl = '' - if self.experiment.sampleModel: - modl += f" --sampleModel '{self.experiment.sampleModel}'" - - acts = '' - if self.reduction.autoscale: - acts += f' --autoscale {" ".join(str(ff) for ff in self.reduction.autoscale)}' - if self.reduction.scale != Defaults.scale: - acts += f' --scale {self.reduction.scale}' - if self.reduction.timeSlize: - acts += f' --timeSlize {" ".join(str(ff) for ff in self.reduction.timeSlize)}' - - mlst = base + inpt + otpt - if mask: - mlst += mask - if para: - mlst += para - if acts: - mlst += acts - if modl: - mlst += modl - - if len(mlst) > 70: - mlst = base + ' ' + inpt + ' ' + otpt - if mask: - mlst += ' ' + mask - if para: - mlst += ' ' + para - if acts: - mlst += ' ' + acts - if modl: - mlst += ' ' + modl - - logging.debug(f'Argument list build in EOSConfig.call_string: {mlst}') - return mlst - - diff --git a/libeos/reduction.py b/libeos/reduction.py deleted file mode 100644 index 32d7c1f..0000000 --- a/libeos/reduction.py +++ /dev/null @@ -1,505 +0,0 @@ -import logging -import os -import sys - -import numpy as np -from orsopy import fileio - -from .command_line import expand_file_list -from .file_reader import AmorData -from .header import Header -from .options import EOSConfig -from .instrument import Grid - -class AmorReduction: - def __init__(self, config: EOSConfig): - self.experiment_config = config.experiment - self.reader_config = config.reader - self.reduction_config = config.reduction - self.output_config = config.output - self.grid = Grid(config.reduction.qResolution, config.experiment.qzRange) - - self.header = Header() - self.header.reduction.call = config.call_string() - - self.monitorUnit = {'n': 'cnts', 'p': 'mC', 't': 's', 'auto': 'pulses'} - - def reduce(self): - if not os.path.exists(f'{self.output_config.outputPath}'): - logging.debug(f'Creating destination path {self.output_config.outputPath}') - os.system(f'mkdir {self.output_config.outputPath}') - - # load or create normalisation matrix - if self.reduction_config.normalisationFileIdentifier: - self.create_normalisation_map(self.reduction_config.normalisationFileIdentifier[0]) - else: - self.norm_lz = self.grid.lz() - self.normAngle = 1. - self.normMonitor = 1. - - logging.warning('normalisation matrix: none requested') - - # load R(q_z) curve to be subtracted: - if self.reduction_config.subtract: - self.sq_q, self.sR_q, self.sdR_q, self.sFileName = self.loadRqz(self.reduction_config.subtract) - logging.warning(f'loaded background file: {self.sFileName}') - self.header.reduction.corrections.append(f'background from \'{self.sFileName}\' subtracted') - self.subtract = True - else: - self.subtract = False - - # load measurement data and do the reduction - self.datasetsRqz = [] - self.datasetsRlt = [] - for i, short_notation in enumerate(self.reduction_config.fileIdentifier): - self.read_file_block(i, short_notation) - - # output - logging.warning('output:') - - if 'Rqz.ort' in self.output_config.outputFormats: - self.save_Rqz() - - if 'Rlt.ort' in self.output_config.outputFormats: - self.save_Rtl() - - def read_file_block(self, i, short_notation): - logging.warning('reading input:') - self.header.measurement_data_files = [] - self.file_reader = AmorData(header=self.header, - reader_config=self.reader_config, - config=self.experiment_config, - short_notation=short_notation) - if self.reduction_config.timeSlize: - self.read_timeslices(i) - else: - self.read_unsliced(i) - - def read_unsliced(self, i): - lamda_e = self.file_reader.lamda_e - detZ_e = self.file_reader.detZ_e - self.monitor = np.sum(self.file_reader.monitorPerPulse) - logging.warning(f' monitor = {self.monitor:8.2f} {self.monitorUnit[self.experiment_config.monitorType]}') - qz_lz, qx_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, self.mask_lz = self.project_on_lz( - self.file_reader, self.norm_lz, self.normAngle, lamda_e, detZ_e) - #if self.monitor>1 : - # ref_lz /= self.monitor - # err_lz /= self.monitor - try: - ref_lz *= self.reduction_config.scale[i] - err_lz *= self.reduction_config.scale[i] - except IndexError: - ref_lz *= self.reduction_config.scale[-1] - err_lz *= self.reduction_config.scale[-1] - if 'Rqz.ort' in self.output_config.outputFormats: - headerRqz = self.header.orso_header() - headerRqz.data_set = f'Nr {i} : mu = {self.file_reader.mu:6.3f} deg' - - if qz_lz[0,int(np.shape(qz_lz)[1]/2)] < 0: - # assuming a 'measurement from below' when center of detector at negative qz - qz_lz *= -1 - - # projection on qz-grid - q_q, R_q, dR_q, dq_q = self.project_on_qz(qz_lz, ref_lz, err_lz, res_lz, self.norm_lz, self.mask_lz) - - # The filtering is now done by restricting the qz-grid - #filter_q = np.where((self.experiment_config.qzRange[0]>q_q) & (q_q>self.experiment_config.qzRange[1]), - # False, True) - #q_q = q_q[filter_q] - #R_q = R_q[filter_q] - #dR_q = dR_q[filter_q] - #dq_q = dq_q[filter_q] - - if self.reduction_config.autoscale: - if i==0: - R_q, dR_q = self.autoscale(q_q, R_q, dR_q) - else: - pRq_z = self.datasetsRqz[i-1].data[:, 1] - pdRq_z = self.datasetsRqz[i-1].data[:, 2] - R_q, dR_q = self.autoscale(q_q, R_q, dR_q, pRq_z, pdRq_z) - - if self.subtract: - if len(q_q)==len(self.sq_q): - R_q -= self.sR_q - dR_q = np.sqrt(dR_q**2+self.sdR_q**2) - else: - logging.warning( - f'backgroung file {self.sFileName} not compatible with q_z scale ({len(self.sq_q)} vs. {len(q_q)})') - - data = np.array([q_q, R_q, dR_q, dq_q]).T - orso_data = fileio.OrsoDataset(headerRqz, data) - self.datasetsRqz.append(orso_data) - if 'Rlt.ort' in self.output_config.outputFormats: - columns = [ - fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), - fileio.Column('R', '', 'specular reflectivity'), - fileio.ErrorColumn(error_of='R', error_type='uncertainty', value_is='sigma'), - fileio.ErrorColumn(error_of='Qz', error_type='resolution', value_is='sigma'), - fileio.Column('lambda', 'angstrom', 'wavelength'), - fileio.Column('alpha_f', 'deg', 'final angle'), - fileio.Column('l', '', 'index of lambda-bin'), - fileio.Column('t', '', 'index of theta bin'), - fileio.Column('intensity', '', 'filtered neutron events per pixel'), - fileio.Column('norm', '', 'normalisation matrix'), - fileio.Column('mask', '', 'pixels used for calculating R(q_z)'), - fileio.Column('Qx', '1/angstrom', 'parallel momentum transfer'), - ] - # data_source = file_reader.data_source - - ts, zs = ref_lz.shape - lindex_lz = np.tile(np.arange(1, ts+1), (zs, 1)).T - tindex_lz = np.tile(np.arange(1, zs+1), (ts, 1)) - - j = 0 - for item in zip( - qz_lz.T, - ref_lz.T, - err_lz.T, - res_lz.T, - lamda_lz.T, - theta_lz.T, - lindex_lz.T, - tindex_lz.T, - int_lz.T, - self.norm_lz.T, - np.where(self.mask_lz, 1, 0).T, - qx_lz.T, - ): - data = np.array(list(item)).T - headerRlt = self.header.orso_header(columns=columns) - headerRlt.data_set = f'dataset_{i}_{j+1} : alpha_f = {theta_lz[0, j]:6.3f} deg' - orso_data = fileio.OrsoDataset(headerRlt, data) - self.datasetsRlt.append(orso_data) - j += 1 - - def read_timeslices(self, i): - wallTime_e = np.float64(self.file_reader.wallTime_e)/1e9 - pulseTimeS = np.float64(self.file_reader.pulseTimeS)/1e9 - interval = self.reduction_config.timeSlize[0] - try: - start = self.reduction_config.timeSlize[1] - except IndexError: - start = 0 - try: - stop = self.reduction_config.timeSlize[2] - except IndexError: - stop = wallTime_e[-1] - # make overwriting log lines possible by removing newline at the end - #logging.StreamHandler.terminator = "\r" - logging.warning(f' time slizing') - logging.info(' slize time monitor') - for ti, time in enumerate(np.arange(start, stop, interval)): - - filter_e = np.where((time0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) - else: - logging.warning(' automatic scaling not possible') - scale = 1. - else: - filter_q = np.where(np.isnan(pR_q*R_q), False, True) - filter_q = np.where(R_q>0, filter_q, False) - filter_q = np.where(pR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ - / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) - else: - logging.warning(' automatic scaling not possible') - scale = 1. - R_q /= scale - dR_q /= scale - logging.info(f' scaling factor = {1/scale}') - - return R_q, dR_q - - def project_on_qz(self, q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz): - q_q = self.grid.q() - mask_lzf = mask_lz.flatten() - q_lzf = q_lz.flatten()[mask_lzf] - R_lzf = R_lz.flatten()[mask_lzf] - dR_lzf = dR_lz.flatten()[mask_lzf] - dq_lzf = dq_lz.flatten()[mask_lzf] - norm_lzf = norm_lz.flatten()[mask_lzf] - - weights_lzf = norm_lzf - #weights_lzf = np.sqrt(norm_lzf) - #weights_lzf = 1 / dR_lzf - - N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - - R_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf * R_lzf )[0] - R_q = R_q / N_q - - dR_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dR_lzf)**2 )[0] - dR_q = np.sqrt( dR_q ) / N_q - - # TODO: different error propagations for dR and dq! - # this is what should work: - #dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0] - #dq_q = np.sqrt( dq_q ) / N_q - # and this actually works: - N_q = np.histogram(q_lzf, bins = q_q, weights = weights_lzf**2 )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - dq_q = np.histogram(q_lzf, bins = q_q, weights = (weights_lzf * dq_lzf)**2 )[0] - dq_q = np.sqrt( dq_q / N_q ) - - q_q = 0.5 * (q_q + np.roll(q_q, 1)) - - return q_q[1:], R_q, dR_q, dq_q - - def loadRqz(self, name): - fname = os.path.join(self.output_config.outputPath, name) - if os.path.exists(fname): - fileName = fname - elif os.path.exists(f'{fname}.Rqz.ort'): - fileName = f'{fname}.Rqz.ort' - else: - sys.exit(f'### the background file \'{fname}\' does not exist! => stopping') - - q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) - - return q_q, Sq_q, dS_q, fileName - - def create_normalisation_map(self, short_notation): - outputPath = self.output_config.outputPath - normalisation_list = expand_file_list(short_notation) - name = str(normalisation_list[0]) - for i in range(1, len(normalisation_list), 1): - name = f'{name}_{normalisation_list[i]}' - n_path = os.path.join(outputPath, f'{name}.norm') - if os.path.exists(n_path): - logging.warning(f'normalisation matrix: found and using {n_path}') - with open(n_path, 'rb') as fh: - self.normFileList = np.load(fh, allow_pickle=True) - self.normAngle = np.load(fh, allow_pickle=True) - self.norm_lz = np.load(fh, allow_pickle=True) - self.normMonitor = np.load(fh, allow_pickle=True) - for i, entry in enumerate(self.normFileList): - self.normFileList[i] = entry.split('/')[-1] - self.header.measurement_additional_files = self.normFileList - else: - logging.warning(f'normalisation matrix: using the files {normalisation_list}') - fromHDF = AmorData(header=self.header, - reader_config=self.reader_config, - config=self.experiment_config, - short_notation=short_notation, norm=True) - self.normAngle = fromHDF.nu - fromHDF.mu - lamda_e = fromHDF.lamda_e - detZ_e = fromHDF.detZ_e - self.normMonitor = np.sum(fromHDF.monitorPerPulse) - norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (self.grid.lamda(), self.grid.z())) - norm_lz = np.where(norm_lz>2, norm_lz, np.nan) - if self.reduction_config.normalisationMethod == 'd': - # direct reference => invert map vertically - self.norm_lz = np.flip(norm_lz, 1) - else: - # correct for reference sm reflectivity - lamda_l = self.grid.lamda() - theta_z = self.normAngle + fromHDF.delta_z - lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T - theta_lz = self.grid.lz()*theta_z - qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz - # TODO: introduce variable for `m` and propably for the slope - Rsm_lz = np.ones(np.shape(qz_lz)) - Rsm_lz = np.where(qz_lz>0.0217, 1-(qz_lz-0.0217)*(0.0625/0.0217), Rsm_lz) - Rsm_lz = np.where(qz_lz>0.0217*5, np.nan, Rsm_lz) - self.norm_lz = norm_lz / Rsm_lz - - if len(lamda_e) > 1e6: - with open(n_path, 'wb') as fh: - np.save(fh, np.array(fromHDF.file_list), allow_pickle=False) - np.save(fh, np.array(self.normAngle), allow_pickle=False) - np.save(fh, self.norm_lz, allow_pickle=False) - np.save(fh, self.normMonitor, allow_pickle=False) - self.normFileList = fromHDF.file_list - self.header.reduction.corrections.append('normalisation with \'additional files\'') - - def project_on_lz(self, fromHDF, norm_lz, normAngle, lamda_e, detZ_e): - # projection on lambda-z-grid - lamda_l = self.grid.lamda() - alphaF_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z - # TODO: implement various methods to obtain alpha_i. - #if self.experiment_config.incidentAngle == 'alphaF': - # # for specular reflectometry with a highly divergent beam - # alphaF_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z - #elif self.experiment_config.incidentAngle == 'nu': - # # for specular reflectometry, using kappa nad nu but ignoring mu - # alphaF_z = (fromHDF.nu + fromHDF.delta_z + fromHDF.kap + fromHDF.kad) / 2. - #else: - # # using kappa, for a collimated incoming beam - # pass - lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T - alphaF_lz = self.grid.lz()*alphaF_z - - mask_lz = np.where(np.isnan(norm_lz), False, True) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(alphaF_lz)>5e-3, True, False)) - if self.reduction_config.thetaRangeR[1]<12: - t0 = fromHDF.nu - fromHDF.mu - mask_lz = np.logical_and(mask_lz, np.where(alphaF_lz-t0 >= self.reduction_config.thetaRangeR[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(alphaF_lz-t0 <= self.reduction_config.thetaRangeR[1], True, False)) - elif self.reduction_config.thetaRange[1]<12: - mask_lz = np.logical_and(mask_lz, np.where(alphaF_lz >= self.reduction_config.thetaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(alphaF_lz <= self.reduction_config.thetaRange[1], True, False)) - else: - self.reduction_config.thetaRange = [fromHDF.nu - fromHDF.mu - fromHDF.div/2, - fromHDF.nu - fromHDF.mu + fromHDF.div/2] - mask_lz = np.logical_and(mask_lz, np.where(alphaF_lz >= self.reduction_config.thetaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(alphaF_lz <= self.reduction_config.thetaRange[1], True, False)) - if self.experiment_config.lambdaRange[1]<15: - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= self.experiment_config.lambdaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= self.experiment_config.lambdaRange[1], True, False)) - - # gravity correction - #alphaF_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) - alphaF_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) - - if self.experiment_config.incidentAngle == 'alphaF': - #alphaI_lz = alphaF_lz - qz_lz = 4.0*np.pi * np.sin(np.deg2rad(alphaF_lz)) / lamda_lz - qx_lz = self.grid.lz() * 0. - else: - alphaI_lz = self.grid.lz()*(fromHDF.mu + fromHDF.kap + fromHDF.kad) - qz_lz = 2.0*np.pi * (np.sin(np.deg2rad(alphaF_lz)) + np.sin(np.deg2rad(alphaI_lz))) / lamda_lz - qx_lz = 2.0*np.pi * (np.cos(np.deg2rad(alphaF_lz)) - np.cos(np.deg2rad(alphaI_lz))) / lamda_lz - - int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, self.grid.z())) - # cut normalisation sample horizon - int_lz = np.where(mask_lz, int_lz, np.nan) - thetaF_lz = np.where(mask_lz, alphaF_lz, np.nan) - - if self.reduction_config.normalisationMethod == 'o': - logging.debug(' assuming an overilluminated sample and correcting for the angle of incidence') - thetaN_z = fromHDF.delta_z + normAngle - thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z - thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) - ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) - elif self.reduction_config.normalisationMethod == 'u': - logging.debug(' assuming an underilluminated sample and ignoring the angle of incidence') - ref_lz = (int_lz / norm_lz) - elif self.reduction_config.normalisationMethod == 'd': - logging.debug(' assuming direct beam for normalisation and ignoring the angle of incidence') - ref_lz = (int_lz / norm_lz) - else: - logging.error('unknown normalisation method! Use [u]nder, [o]ver or [d]irect illumination') - ref_lz = (int_lz / norm_lz) - if self.monitor > 1e-6 : - ref_lz *= self.normMonitor / self.monitor - else: - logging.info(' too small monitor value for normalisation -> ignoring monitors') - err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) - - # TODO: allow for non-ideal Delta lambda / lambda (rather than 2.2%) - res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(alphaF_z)[0])) * 0.022**2 - res_lz = res_lz + (0.008/alphaF_lz)**2 - res_lz = qz_lz * np.sqrt(res_lz) - - return qz_lz, qx_lz, ref_lz, err_lz, res_lz, lamda_lz, alphaF_lz, int_lz, mask_lz - - - @staticmethod - def histogram2d_lz(lamda_e, detZ_e, bins): - """ - Perform binning operation equivalent to numpy bin2d for the sepcial case - of the second dimension using integer positions (pre-defined pixels). - Based on the devide_bin algorithm below. - """ - dimension = bins[1].shape[0]-1 - if not (np.array(bins[1])==np.arange(dimension+1)).all(): - raise ValueError("histogram2d_lz requires second bin dimension to be contigous integer range") - binning = AmorReduction.devide_bin(lamda_e, detZ_e.astype(np.int64), bins[0], dimension) - return np.array(binning), bins[0], bins[1] - - @staticmethod - def devide_bin(lambda_e, position_e, lamda_edges, dimension): - ''' - Use a divide and conquer strategy to bin the data. For the actual binning the - numpy bincount function is used, as it is much faster than histogram for - counting of integer values. - - :param lambda_e: Array of wavelength for each event - :param position_e: Array of positional indices for each event - :param lamda_edges: The edges of bins to be used for the histogram - :param dimension: position number of buckets in output arrray - - :return: 2D list of dimensions (lambda, x) of counts - ''' - if len(lambda_e)==0: - # no more events in range, return empty bins - return [np.zeros(dimension, dtype=np.int64).tolist()]*(len(lamda_edges)-1) - if len(lamda_edges)==2: - # deepest recursion reached, all items should be within the two ToF edges - return [np.bincount(position_e, minlength=dimension).tolist()] - # split all events into two time of flight regions - split_idx = len(lamda_edges)//2 - left_region = lambda_e=3.8 packages = - libeos -scripts = - eos.py + eos install_requires = numpy h5py @@ -30,3 +28,7 @@ install_requires = [project.urls] Homepage = "https://github.com/jochenstahn/amor" + +[options.entry_points] +console_scripts = + eos = eos.__main__:main diff --git a/test_data/amor2023n000608.hdf b/test_data/amor2023n000608.hdf deleted file mode 100644 index da050fa..0000000 Binary files a/test_data/amor2023n000608.hdf and /dev/null differ diff --git a/test_data/amor2023n000609.hdf b/test_data/amor2023n000609.hdf deleted file mode 100644 index a3c869c..0000000 Binary files a/test_data/amor2023n000609.hdf and /dev/null differ diff --git a/test_data/amor2023n000610.hdf b/test_data/amor2023n000610.hdf deleted file mode 100644 index 204cda8..0000000 Binary files a/test_data/amor2023n000610.hdf and /dev/null differ diff --git a/test_data/amor2023n000611.hdf b/test_data/amor2023n000611.hdf deleted file mode 100644 index 077ecdf..0000000 Binary files a/test_data/amor2023n000611.hdf and /dev/null differ diff --git a/test_data/amor2023n000612.hdf b/test_data/amor2023n000612.hdf deleted file mode 100644 index 3ef5b3c..0000000 Binary files a/test_data/amor2023n000612.hdf and /dev/null differ diff --git a/test_data/amor2023n000613.hdf b/test_data/amor2023n000613.hdf deleted file mode 100644 index 8986f99..0000000 Binary files a/test_data/amor2023n000613.hdf and /dev/null differ diff --git a/test_data/amor2025n005952.hdf b/test_data/amor2025n005952.hdf new file mode 100644 index 0000000..d1d7824 --- /dev/null +++ b/test_data/amor2025n005952.hdf @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6a6c94bb1e1c2c41290bb5c13aa593bbb031dc28662b91dd4cb5f512ca501697 +size 164227616 diff --git a/test_data/amor2025n006003.hdf b/test_data/amor2025n006003.hdf new file mode 100644 index 0000000..3c06f83 --- /dev/null +++ b/test_data/amor2025n006003.hdf @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e4d32f754edc6beda8c6b3943d7962121eb5ec36be8245efa7749776ae97208f +size 511262 diff --git a/test_data/amor2025n006004.hdf b/test_data/amor2025n006004.hdf new file mode 100644 index 0000000..42e5eb7 --- /dev/null +++ b/test_data/amor2025n006004.hdf @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:48d3bea779abb6ec415f051cb77d1d67ff04e6e57e55ed2eacfe18448a7c1056 +size 765639 diff --git a/test_data/amor2025n006005.hdf b/test_data/amor2025n006005.hdf new file mode 100644 index 0000000..adf8e79 --- /dev/null +++ b/test_data/amor2025n006005.hdf @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3e523acca89e18f66d1ab73b823425ea9d0bec5c59d3b467247ea082ecd5ff5a +size 3293620 diff --git a/tests/analyze_hdf.py b/tests/analyze_hdf.py new file mode 100644 index 0000000..7d700ab --- /dev/null +++ b/tests/analyze_hdf.py @@ -0,0 +1,31 @@ +""" +Small helper to find information about hdf datafiles for debugging +""" + +import h5py + +def rec_tree(group, min_size=128): + if hasattr(group, 'keys'): + output = {} + total_size = 0 + for key in group.keys(): + subkeys, size = rec_tree(group[key], min_size) + total_size += size + if size>min_size: + if subkeys: + output[key] = subkeys + else: + output[key] = size + return output, size + elif hasattr(group, 'size'): + return None, group.size + else: + return None, 0 + +if __name__ == "__main__": + import sys + for fi in sys.argv[1:]: + print(fi) + print(rec_tree(sys.argv[1])) + + \ No newline at end of file diff --git a/tests/mock_data.py b/tests/mock_data.py new file mode 100644 index 0000000..8b5f00a --- /dev/null +++ b/tests/mock_data.py @@ -0,0 +1,64 @@ +""" +Generate a mock dataset in memory for running unit tests. +""" + +import h5py +import numpy as np + +MOCK_METADATA = { + 'title': 'Testdata', + 'proposal_id': 'none', + 'user/name': 'test user', + 'user/email': 'test@user.de', + 'sample/name': 'test sample', + 'sample/model': 'air | Fe 12 | Si', + 'Amor/source/name': 'SINQ', + 'start_time': '2025-01-01 00:00:01', + } +MOCK_META_TYPED = { + 'Amor/chopper/pair_separation': (1000.0, np.float32), + 'Amor/detector/transformation/distance': (4000.0, np.float64), + 'Amor/instrument_control_parameters/kappa': (1000.0, np.float64), + 'Amor/instrument_control_parameters/kappa_offset': (1000.0, np.float64), + 'Amor/instrument_control_parameters/div': (1.6, np.float64), + 'Amor/chopper/ch1_trigger_phase': (-9.1, np.float64), + 'Amor/chopper/ch2_trigger_phase': (6.75, np.float64), + 'Amor/chopper/ch2_trigger/event_time_zero': ([0.0]*10, np.uint64), + 'Amor/chopper/ch2_trigger/event_time_offset': ([0.0]*10, np.uint32), + 'Amor/chopper/rotation_speed': (500.0, np.float64), + 'Amor/chopper/phase': (0.0, np.float64), + 'Amor/polarization/configuration/value': (0.0, np.float64), + } + +def mock_data(mu=1.0, nu=2.0): + hdf = h5py.File.in_memory() # requires h5py >=3.13 + ds = hdf.create_group('entry1') + for key, value in MOCK_METADATA.items(): + ds.create_dataset(key, data=np.array([value.encode('utf-8')])) + for key, (value, dtype) in MOCK_META_TYPED.items(): + if type(value) is list: + ds.create_dataset(key, data=np.array(value), dtype=dtype) + else: + ds.create_dataset(key, data=np.array([value]), dtype=dtype) + + ds.create_dataset('Amor/instrument_control_parameters/mu', np.array([mu]), dtype=np.float64) + ds.create_dataset('Amor/instrument_control_parameters/nu', np.array([nu]), dtype=np.float64) + + return hdf + +def compare_with_real_data(fname): + hdf = h5py.File(fname, 'r') + ds = hdf['entry1'] + for key, value in MOCK_METADATA.items(): + try: + ds[key][0].decode('utf-8') + except KeyError: + print(f'/entry1/{key} does not exist in file') + for key, (value, dtype) in MOCK_META_TYPED.items(): + try: + item = ds[key] + except KeyError: + print(f'/entry1/{key} does not exist in file') + else: + if item.dtype != dtype: + print(f'/entry1/{key} does not match {dtype}, dataset is {item.dtype}') diff --git a/tests/test_full_analysis.py b/tests/test_full_analysis.py index a1fc875..2105814 100644 --- a/tests/test_full_analysis.py +++ b/tests/test_full_analysis.py @@ -1,16 +1,27 @@ import os import cProfile from unittest import TestCase -from libeos import options, reduction, logconfig +from dataclasses import fields, MISSING +from eos import options, reduction, logconfig logconfig.setup_logging() -logconfig.update_loglevel(True, False) +logconfig.update_loglevel(1) # TODO: add test for new features like proton charge normalization class FullAmorTest(TestCase): @classmethod def setUpClass(cls): + # generate map for option defaults + cls._field_defaults = {} + for opt in [options.ExperimentConfig, options.ReductionConfig, options.OutputConfig]: + defaults = {} + for field in fields(opt): + if field.default not in [None, MISSING]: + defaults[field.name] = field.default + elif field.default_factory not in [None, MISSING]: + defaults[field.name] = field.default_factory() + cls._field_defaults[opt.__name__] = defaults cls.pr = cProfile.Profile() @classmethod @@ -20,49 +31,47 @@ class FullAmorTest(TestCase): def setUp(self): self.pr.enable() self.reader_config = options.ReaderConfig( - year=2023, - rawPath=(os.path.join('..', "test_data"),), + year=2025, + rawPath=["test_data"], ) def tearDown(self): self.pr.disable() - for fi in ['test.Rqz.ort', '614.norm']: + for fi in ['test_results/test.Rqz.ort', 'test_results/5952.norm']: try: - os.unlink(os.path.join(self.reader_config.rawPath[0], fi)) + os.unlink(fi) except FileNotFoundError: pass - def test_time_slicing(self): experiment_config = options.ExperimentConfig( + chopperSpeed=self._field_defaults['ExperimentConfig']['chopperSpeed'], chopperPhase=-13.5, chopperPhaseOffset=-5, - monitorType=options.Defaults.monitorType, - lowCurrentThreshold=options.Defaults.lowCurrentThreshold, - yRange=(11., 41.), - lambdaRange=(2., 15.), - qzRange=(0.005, 0.30), - incidentAngle=options.Defaults.incidentAngle, + monitorType=self._field_defaults['ExperimentConfig']['monitorType'], + lowCurrentThreshold=self._field_defaults['ExperimentConfig']['lowCurrentThreshold'], + yRange=(18, 48), + lambdaRange=(3., 11.5), + incidentAngle=self._field_defaults['ExperimentConfig']['incidentAngle'], mu=0, nu=0, muOffset=0.0, sampleModel='air | 10 H2O | D2O' ) reduction_config = options.ReductionConfig( - normalisationMethod=options.Defaults.normalisationMethod, + normalisationMethod=self._field_defaults['ReductionConfig']['normalisationMethod'], qResolution=0.01, - qzRange=options.Defaults.qzRange, - thetaRange=(-12., 12.), - thetaRangeR=(-12., 12.), - fileIdentifier=["610"], + qzRange=self._field_defaults['ReductionConfig']['qzRange'], + thetaRange=(-0.75, 0.75), + fileIdentifier=["6003-6005"], scale=[1], normalisationFileIdentifier=[], timeSlize=[300.0] ) output_config = options.OutputConfig( - outputFormats=["Rqz.ort"], + outputFormats=[options.OutputFomatOption.Rqz_ort], outputName='test', - outputPath=os.path.join('..', 'test_results'), + outputPath='test_results', ) config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config) # run three times to get similar timing to noslicing runs @@ -75,33 +84,32 @@ class FullAmorTest(TestCase): def test_noslicing(self): experiment_config = options.ExperimentConfig( + chopperSpeed=self._field_defaults['ExperimentConfig']['chopperSpeed'], chopperPhase=-13.5, chopperPhaseOffset=-5, - monitorType=options.Defaults.monitorType, - lowCurrentThreshold=options.Defaults.lowCurrentThreshold, - yRange=(11., 41.), - lambdaRange=(2., 15.), - qzRange=(0.005, 0.30), - incidentAngle=options.Defaults.incidentAngle, + monitorType=self._field_defaults['ExperimentConfig']['monitorType'], + lowCurrentThreshold=self._field_defaults['ExperimentConfig']['lowCurrentThreshold'], + yRange=(18, 48), + lambdaRange=(3., 11.5), + incidentAngle=self._field_defaults['ExperimentConfig']['incidentAngle'], mu=0, nu=0, - muOffset=0.0 + muOffset=0.0, ) reduction_config = options.ReductionConfig( + normalisationMethod=self._field_defaults['ReductionConfig']['normalisationMethod'], qResolution=0.01, - qzRange=options.Defaults.qzRange, - normalisationMethod=options.Defaults.normalisationMethod, - thetaRange=(-12., 12.), - thetaRangeR=(-12., 12.), - fileIdentifier=["610", "611", "608,612-613", "609"], + qzRange=self._field_defaults['ReductionConfig']['qzRange'], + thetaRange=(-0.75, 0.75), + fileIdentifier=["6003", "6004", "6005"], scale=[1], - normalisationFileIdentifier=["608"], - autoscale=(True, True) + normalisationFileIdentifier=["5952"], + autoscale=(0.0, 0.05), ) output_config = options.OutputConfig( - outputFormats=["Rqz.ort"], + outputFormats=[options.OutputFomatOption.Rqz_ort], outputName='test', - outputPath=os.path.join('..', 'test_results'), + outputPath='test_results', ) config=options.EOSConfig(self.reader_config, experiment_config, reduction_config, output_config) reducer = reduction.AmorReduction(config) diff --git a/windows_build.spec b/windows_build.spec index c933036..47c81c4 100644 --- a/windows_build.spec +++ b/windows_build.spec @@ -2,7 +2,7 @@ a = Analysis( - ['eos.py'], + ['eos/__main__.py'], pathex=[], binaries=[], datas=[], diff --git a/windows_folder.spec b/windows_folder.spec index 7646535..4ca0124 100644 --- a/windows_folder.spec +++ b/windows_folder.spec @@ -9,7 +9,7 @@ datas += tmp_ret[0]; binaries += tmp_ret[1]; hiddenimports += tmp_ret[2] a = Analysis( - ['eos.py'], + ['eos/__main__.py'], pathex=[], binaries=[], datas=[],