diff --git a/eos/options.py b/eos/options.py index 6ae0b7c..ec54062 100644 --- a/eos/options.py +++ b/eos/options.py @@ -593,6 +593,7 @@ class ReflectivityConfig: class E2HPlotSelection(StrEnum): All = 'all' + Raw = 'raw' YZ = 'Iyz' LT = 'Ilt' YT = 'Iyt' @@ -700,6 +701,15 @@ class E2HReductionConfig(ArgParsable): }, ) + fontsize: float = field( + default=8., + metadata={ + 'short': 'pf', + 'group': 'output', + 'help': 'font size for graphs', + }, + ) + @dataclass class E2HConfig: reader: ReaderConfig diff --git a/eos/reduction_e2h.py b/eos/reduction_e2h.py index c93ba67..71762c0 100644 --- a/eos/reduction_e2h.py +++ b/eos/reduction_e2h.py @@ -48,11 +48,16 @@ class E2HReduction: self.file_index = 0 self.plot_kwds = {} self.fig = plt.figure() + plt.rcParams.update({'font.size': self.config.reduction.fontsize}) if self.config.reduction.update: # live update implies plotting self.config.reduction.show_plot = True + if self.config.reduction.plot==E2HPlotSelection.Raw: + # Raw implies fast caculations + self.config.reduction.fast = True + if not self.config.reduction.fast or self.config.reduction.plot in NEEDS_LAMDA: from . import event_analysis as ea @@ -89,9 +94,9 @@ class E2HReduction: if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, E2HPlotSelection.Q]: self.grid = LZGrid(0.01, [0.0, 0.25]) - if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.LT, - E2HPlotSelection.YZ, E2HPlotSelection.YT, - E2HPlotSelection.TZ]: + if self.config.reduction.plot in [E2HPlotSelection.All, E2HPlotSelection.Raw, + E2HPlotSelection.LT, E2HPlotSelection.YT, + E2HPlotSelection.YZ, E2HPlotSelection.TZ]: self.plot_kwds['colorbar'] = True self.plot_kwds['cmap'] = str(self.config.reduction.plot_colormap) if self.config.reduction.plotArgs==E2HPlotArguments.Linear: @@ -193,6 +198,17 @@ class E2HReduction: self.projection = CombinedProjection(3, 2, [plz, pyz, pr], [(0, 2, 0, 1), (0, 2, 1, 2), (2, 3, 0, 2)]) + if self.config.reduction.plot==E2HPlotSelection.Raw: + del(self.plot_kwds['colorbar']) + # A combined graph that does not require longer calculations + plyt = YTProjection(tthh) + pltofz = TofZProjection(last_file_header.timing.tau, foldback=not self.config.reduction.fast) + pltof = TofProjection(last_file_header.timing.tau, foldback=not self.config.reduction.fast) + plt = TProjection(tthh) + + self.projection = CombinedProjection(3, 3, [plyt, pltofz, plt, pltof], + [(0,2, 0, 1), (0, 2, 1, 3), (2,3, 0,1),(2,3,1,3)]) + def read_data(self): fileName = self.file_list[self.file_index] self.dataset = AmorEventData(fileName, max_events=self.config.reduction.max_events) diff --git a/events2histogram.py b/events2histogram.py deleted file mode 100755 index 29a8d6c..0000000 --- a/events2histogram.py +++ /dev/null @@ -1,1056 +0,0 @@ -__version__ = '2025-06-07' - -# 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 -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 - - - # 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/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/' - 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] - 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] - 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) - - self.tau = 30. / chopperSpeed - - 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 - - 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 - - chopperTriggerPhase = clas.chopperTriggerPhase - 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 - - 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 - #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", "--chopperTriggerPhase", - 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() - -