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