Loading pysen/config.py +0 −1 Original line number Diff line number Diff line Loading @@ -325,7 +325,6 @@ def coverage(lmax, qmin, **kwargs): for j,y in enumerate(ylims): theta1 = get_theta_pix(x,y,L2,sa,ca) theta2, _ = pixel_angle(i,j,theta0) print(theta1, theta2) theta = theta2 lq = get_q(lmax, theta)*lmax q = linspace(get_q(lmax, theta), get_q(lmax-dlam, theta), len(tbins)+1) Loading pysen/revision.py +4 −4 Original line number Diff line number Diff line Loading @@ -2,9 +2,9 @@ PySEN revision module """ import sys __version__ = "0.7.3" __version__ = "0.7.4" __release__ = "dev1" __date__ = "Feb 13, 2019" __date__ = "Feb 14, 2019" # VERSION = __version__ RELEASE = __release__ Loading @@ -28,8 +28,8 @@ def version(full=False): "get pysen version number" if full: py_ver = "Python %s.%s" % (sys.version_info.major, sys.version_info.minor) return "%s.%s (%s) %s" % (__version__, __release__, __date__, py_ver) return "%s.%s" % (__version__, __release__) return "%s-%s (%s) %s" % (__version__, __release__, __date__, py_ver) return "%s-%s" % (__version__, __release__) # EOF scripts/nseplot.py 0 → 100755 +476 −0 Original line number Diff line number Diff line #!/usr/bin/env python "Echo plot" import os.path import logging import numpy as np import h5py import matplotlib.pyplot as plt from pysen import ANGSTROM, NANOSECOND, GAUSS, MICRO, OMEGA_N, version from pysen.plotutil import get_nsubplots from pysen.io import echo_to_hdf from pysen.atarilib import fit_echo_current, Spectrum MAX_CHI2_DEFAULT = 1e-2 def plot_field(ax,t,v,**kwargs): "plot mag field" label = kwargs.get('label' , '') max_chi2 = kwargs.get('max_chi2', MAX_CHI2_DEFAULT) v0 = np.sqrt(np.sum(v**2,axis=-1)) # a,b = np.polyfit(t,v0,deg=1) v1 = np.polyval([a,b],t) chi2 = np.sum((v0-v1)**2) cov = np.cov(t,v0) diag = cov[0,0]*cov[1,1] r = cov[1,0]/np.sqrt(diag) if diag>0 else np.nan #r = np.corrcoef(t,v0)[1,0] #maxd = max(abs(v0-v1)) descr1 = r'a=%.2e b=%.3g' % (a, b) descr2 = r'$\rho$=%.3g $\chi^2$=%.1e' % (r, chi2) #descr2 = r'$\rho$=%.3g $\chi^2$=%.1e $(\Delta B)_{max}$=%.3g' % (r, chi2, maxd) # va = np.average(v0) vd = max(abs(v0-va)) vmin = va-max(vd,0.5) vmax = va+max(vd,0.5) ax.plot(t,v0, '.--', lw=1, label=label) ax.plot(t,v1, '-') ax.set_ylim(bottom=vmin, top=vmax) ax.text(0.05,0.20,descr1, transform=ax.transAxes, fontsize=7) ax.text(0.05,0.10,descr2, transform=ax.transAxes, fontsize=7) ax.grid(True, lw=0.5) ax.legend(loc='upper left') if chi2>max_chi2: ax.set_facecolor('yellow') def magnetic_fields_plot(hdfile, iecho=0, **kwargs): "plot magnetic fields" # max_chi2 = kwargs.get('max_chi2', MAX_CHI2_DEFAULT) # base = os.path.splitext(os.path.basename(hdfile.filename))[0] comment = hdfile.attrs['master_comment'] comment = comment.decode("utf-8") sample = comment.split()[0] nph = hdfile['/'].attrs['point_to_down'] for echo in list(hdfile['/data'].values()): if iecho and echo.attrs['id'] != iecho: continue phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND q0 = physics.attrs['hkl'][0] lam0 = physics.attrs['lambda'] cur0 = float(tech.attrs['i5'][0]) phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) bfield = phase['bfield'] cur = phase['phase_current'][:, 0] # actual value dj = np.radians(phasesens*(cur-cur0))/OMEGA_N/lam0 nxplot, nyplot = get_nsubplots(len(bfield)) fig, axes = plt.subplots(nxplot,nyplot, figsize=(9,7)) fig.subplots_adjust(hspace=0.4, wspace=0.4) fig.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ $\tau$=%.3fns' % (sample, base, lam0/ANGSTROM, q0, tau0)) for iplt, bsensor in enumerate(bfield): bmag = bfield[bsensor] label = bsensor.split('.')[-1] ax = axes[iplt//nxplot,iplt%nxplot] plot_field(ax, dj[:nph]/MICRO, bmag[:nph,:]*GAUSS/MICRO, label=label, max_chi2=max_chi2) for iplt in range(iplt+1, nxplot*nyplot): ax = axes[iplt//nxplot,iplt%nxplot] ax.axis('off') def atari_plot(hdfile, iecho=0, **kwargs): "atari plot" # only_echo = kwargs.pop('only_echo', False) # tbin1 = kwargs.pop('tbin1',0) # TOF bins tbin2 = kwargs.pop('tbin2',None) # FIXME: HARDCODED xpix1 = kwargs.pop('xpix1', 10) xpix2 = kwargs.pop('xpix2', 22) ypix1 = kwargs.pop('ypix1', 10) ypix2 = kwargs.pop('ypix2', 22) #ntaus = hdfile.attrs['scan_length'] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nt = hdfile['/detector'].attrs['no_t_channels'] ny = hdfile['/detector'].attrs['no_y_channels'] nx = hdfile['/detector'].attrs['no_x_channels'] nph = n_idx['dn'] comment = hdfile.attrs['master_comment'] for echo in list(hdfile['/data'].values()): if iecho and echo.attrs['id'] != iecho: continue phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND q0 = physics.attrs['hkl'][0] lam0 = physics.attrs['lambda'] phase0 = float(tech.attrs['i5'][0]) try: lmax = params['lambdatable'] except KeyError: lmax = params['lambdaTable'] lmax = np.array(lmax).flatten() dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) det = phase['detector'][...] pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha0 = np.average(pcha) pcha = np.tile(pcha, (ny, nx, nt)).reshape(n_idx['nphases'],ny,nx,nt) pcha = pcha/float(pcha0) det = det/pcha # plt.figure(figsize=(8,8)) # plt.subplot(2,2,1) wlen = np.sum(det, axis=(0,1,2)) plt.step(wlen[tbin1:tbin2], lmax[tbin1:tbin2]/ANGSTROM) plt.ylabel(r'$\lambda$ [$\AA$]') plt.xlim(left=0) plt.grid(True) plt.subplot(2,2,2) xatari = np.sum(det, axis=(1,2)).T extent = (cur[0], cur[nph-1]) if tbin2 is None: extent = extent + (lmax[tbin1]/ANGSTROM, max(lmax)/ANGSTROM) else: extent = extent + (lmax[tbin1]/ANGSTROM, lmax[tbin2]/ANGSTROM) plt.imshow(xatari[tbin1:tbin2,:nph], origin='lower', extent=extent, aspect='auto') plt.subplot(2,2,3) img = np.sum(det[:,:,:,tbin1:tbin2], axis=(0,-1)) plt.imshow(img , aspect='auto') plt.axvline(xpix1, color='k') plt.axvline(xpix2-1, color='k') plt.axhline(ypix1, color='k') plt.axhline(ypix2-1, color='k') plt.subplot(2,2,4) pha = det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2] dn = np.sum(pha[nph:n_idx['up']], axis=(1,2,3)) up = np.sum(pha[n_idx['up']:], axis=(1,2,3)) epha = np.sqrt(np.sum(pha[:nph]**2, axis=(1,2,3))) pha = np.sum(pha[:nph], axis=(1,2,3)) spectrum = Spectrum(lam=lmax, dlam=dlam, flux=wlen) phase_ef, (xef, yef), res = fit_echo_current(cur[:nph], (pha, epha), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) res['up'] = np.average(up),np.sqrt(np.sum(up)/len(up)) res['dn'] = np.average(dn),np.sqrt(np.sum(dn)/len(dn)) #for key in res: # print(key, res[key]) p = plt.errorbar(cur[:nph], pha, yerr=epha, fmt='.', lw=1) dcol = p[0].get_color() plt.plot(xef, yef, '--') plt.axvline(phase0, color='blue', lw=1, ls='--', label='%.3fA' % phase0 ) plt.axvline(phase_ef, color='red', lw=2, ls='-', label='%.3fA' % phase_ef) ax = plt.gca() if not only_echo: R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] ) plt.errorbar(cur[nph:n_idx['up']], dn, fmt='.', color=dcol) plt.errorbar(cur[n_idx['up']:], up, fmt='.', color=dcol) plt.text(0.01, 1.01, r'$R(Q,t)$=%.3g' % R, transform=ax.transAxes) plt.axhline(np.average(up)) plt.axhline(np.average(dn)) plt.xlabel(r'phase current [A]') plt.legend(loc='best') plt.grid(True) plt.suptitle(r'%s $Q$=%.3f$\AA^{-1}$ $\tau$=%.3gns' % (comment.decode("utf-8"), q0,tau0)) def echo_plot(hdfile, iecho=0, **kwargs): "echo plot" # center_only = kwargs.pop('center_only', False) # npix = kwargs.pop('npix', 2) # pix # tbin1 = kwargs.pop('tbin1',0) # TOF bins tbin2 = kwargs.pop('tbin2',None) # min_counts = kwargs.pop('min_counts', 10.0) max_chi2 = kwargs.pop('max_chi2' ,500.0) min_amp = kwargs.pop('min_amp' , 0.20) # FIXME: HARDCODED #xpix1 = kwargs.pop('xpix1', 10) #xpix2 = kwargs.pop('xpix2', 22) #ypix1 = kwargs.pop('ypix1', 10) #ypix2 = kwargs.pop('ypix2', 22) base = os.path.splitext(os.path.basename(hdfile.filename))[0] comment = hdfile.attrs['master_comment'] comment = comment.decode("utf-8") sample = comment.split()[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nt = hdfile['/detector'].attrs['no_t_channels'] ny = hdfile['/detector'].attrs['no_y_channels'] nx = hdfile['/detector'].attrs['no_x_channels'] nph = n_idx['dn'] if center_only: ntaus = 1 for echo in list(hdfile['/data'].values()): if not iecho or echo.attrs['id'] == iecho: ntaus = ntaus + 1 nxtau, nytau = get_nsubplots(ntaus) fig0 = plt.figure(figsize=(8,8)) fig0.suptitle(r'%s | %s' % (sample, base)) itau = 0 resolution = [] for echo in list(hdfile['/data'].values()): if iecho and echo.attrs['id'] != iecho: continue itau = itau + 1 phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND q0 = physics.attrs['hkl'][0] lam0 = physics.attrs['lambda'] phase0 = float(tech.attrs['i5'][0]) try: lmax = params['lambdatable'] except KeyError: lmax = params['lambdaTable'] lmax = np.array(lmax).flatten() dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) det = phase['detector'][...] pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha0 = np.average(pcha) pcha = np.tile(pcha, (ny, nx, nt)).reshape(n_idx['nphases'],ny,nx,nt) pcha = pcha/float(pcha0) det = det/pcha wlen = np.sum(det, axis=(0,1,2)) spectrum = Spectrum(lam=lmax, dlam=dlam, flux=wlen) y = np.sum(det[:,:,:,tbin1:tbin2], axis=(1,2,3)) ey = np.sqrt(y) dn = y[nph:n_idx['up']] up = y[n_idx['up']:] phase_ef0, (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) if center_only: ax = fig0.add_subplot(nxtau, nytau, itau) fig0.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ ' % (sample, base, lam0/ANGSTROM, q0)) ax.errorbar(cur[:nph], y[:nph], yerr=ey[:nph], fmt='k.', lw=2, ms=3) ax.plot(xef, yef, 'r-', label=r'$\tau$=%.3gns' % tau0) ax.axvline(phase_ef0, color='blue', lw=1, ls='--') ax.errorbar(cur[nph:n_idx['up']], dn, fmt='k.', ms=3) ax.errorbar(cur[n_idx['up']: ], up, fmt='k.', ms=3) upave = np.average(up) dnave = np.average(dn) udave = (upave+dnave)/2 res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] ) resolution.append((tau0, R)) ax.legend() ax.grid() #ax.set_xticklabels([]) if itau%nytau!=1: ax.set_yticklabels([]) continue npx, nx, ny, nt = det.shape nyy = ny//npix nxx = nx//npix pha = det.reshape(npx, nyy, npix, nxx, npix, nt) pha = pha.sum(axis=(2, 4)) # pixel rebin pha = np.sum(pha[:,:,:,tbin1:tbin2],axis=-1) # tof summation fig, axes = plt.subplots(nxx,nyy, figsize=(8,8)) fig.subplots_adjust(hspace=0.05, wspace=0.05) ymax = np.amax(pha) ymin = np.amin(pha) for j in range(nyy): for i in range(nxx): ax = axes[i,j] ax.set_xticklabels([]) ax.set_yticklabels([]) y = pha[:,i,j] ey = np.sqrt(pha[:,i,j]) dn = y[nph:n_idx['up']] up = y[n_idx['up']:] # upave = np.average(up) dnave = np.average(dn) udave = (upave+dnave)/2 ave = np.average(y[:nph]) std = np.std(y[:nph]) if udave<min_counts or ave<min_counts: continue if udave<min_counts and std/ave<min_wiggle: continue res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens, phase_pf=phase_ef0) phase_ef, (xef, yef), res = res # if res['chi2']>max_chi2: continue res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] ) if abs(R)<min_amp: continue ax.errorbar(cur[:nph], y[:nph], yerr=ey[:nph], fmt='k.', lw=1, ms=1, label='(%s,%s)' % (i,j)) ax.set_ylim(bottom=ymin, top=ymax) ax.plot(xef, yef, 'r-') ax.axvline(phase_ef, color='blue', lw=1, ls='--') ax.errorbar(cur[nph:n_idx['up']], dn, fmt='k.', ms=1) ax.errorbar(cur[n_idx['up']: ], up, fmt='k.', ms=1) if R<0: ax.set_facecolor('yellow') #ax.legend() # plt.text(0.01, 1.01, r'$R(Q,t)$=%.3g' % R, transform=ax.transAxes) #ax.axhline(np.average(up)) #ax.axhline(np.average(dn)) #plt.xlabel(r'ephase current [A]') #plt.legend(loc='best') #plt.grid(True) fig.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ $\tau$=%.3gns' % (sample, base, lam0/ANGSTROM, q0,tau0)) if center_only: # figr = plt.figure(figsize=(8,8)) # figr.suptitle(r'%s | %s' % (sample, base)) ax = fig0.add_subplot(nxtau,nytau,ntaus) res = np.asarray(resolution) print(res) ax.plot(res[:,0],res[:,1], '.--') ax.set_xscale('log') ax.set_yscale('log') ax.set_ylim(top=1.1) ax.grid(True) def main(): "the main" import argparse plot_action = dict(echo=echo_plot,atari=atari_plot, bfield=magnetic_fields_plot) parser = argparse.ArgumentParser(description='echo plot') parser.set_defaults(loglevel=logging.INFO, outdir='.', plot='echo', tbin1=0, tbin2=None, tau=0) parser.add_argument('file', metavar='filename', help='file to process', nargs='+') parser.add_argument('--plot', '-p', dest='plot', choices=['echo', 'atari', 'bfield']) parser.add_argument('--atari', dest='plot', action='store_const', const='atari') parser.add_argument('--bfield', dest='plot', action='store_const', const='bfield') parser.add_argument('--echo' , dest='plot', action='store_const', const='echo') grp_echo = parser.add_argument_group('echo', 'echo plotting') grp_echo.add_argument('--center-only', '-C', dest='center_only', action='store_true', default=False, help='show only center patch') grp_echo.add_argument('--only-echo', dest='only_echo', action='store_true', default=False, help='show only echo (no up/down)') parser.add_argument('--tau', '-t', dest='tau', type=int, help='set tau to display (default=%(default)s)') parser.add_argument('--t1' , '-b', dest='tbin1', type=int, help='set min TOF bin (default=%(default)s)') parser.add_argument('--t2' , '-B', dest='tbin2', type=int, help='set max TOF bin (default=%(default)s)') parser.add_argument('--save-file', '-s', dest='savefile', help='save figure to a file (do not show)') parser.add_argument('--verbose', '-v', dest='loglevel', action='store_const', const=logging.DEBUG, help='verbose output') parser.add_argument('--quiet', '-q', dest='loglevel', action='store_const', const=logging.WARNING, help='verbose output') parser.add_argument('--version', '-V', action='version', version='%(prog)s pysen={version}'.format(version=version(full=True))) args = parser.parse_args() if args.loglevel>=logging.INFO: log_format=r'%(message)s' else: log_format=r'%(levelname)s: %(funcName)s: %(message)s' logging.basicConfig(level=args.loglevel, format=log_format) log = logging.getLogger() if args.savefile: plt.switch_backend('Agg') kwargs = vars(args) for filename in args.file: basename , ext = os.path.splitext(os.path.basename(filename)) if ext==".echo": log.info('converting %s to HDF', basename) hfile = echo_to_hdf(filename, args.outdir) else: hfile = filename with h5py.File(hfile, 'r') as hdf5file: plot_action[args.plot](hdf5file, iecho=args.tau, **kwargs) if args.savefile: plt.savefig(basename+'-'+args.savefile) if not args.savefile: plt.show() if __name__ == "__main__": main() scripts/qtau.py +13 −13 Original line number Diff line number Diff line Loading @@ -52,37 +52,37 @@ def main(): group1 = parser.add_argument_group('script', 'input data are read from the script') group1.add_argument('--file', '-f', dest='file', help='read input data from a file') group2 = parser.add_argument_group('command line options', group2 = parser.add_argument_group('shell', 'input data are taken from command line options') group2.add_argument('taus', metavar='tau', type=float, nargs='*', help='select list of taus') group2.add_argument('--lmax', '-l', dest='lmax', type=float, help='set max wavelength, lambda_max (default %(default)s)') group2.add_argument('--qmin', '-q', dest='qmin', type=float, help='set min q (default %(default)s)') group2.add_argument('--tbin', '-t', dest='tbins', type=int, action='append', help='set tbin (default %(default)s)') help='set TOF binning') # (default %(default)s)') group2.add_argument('--pos' , '-p', dest='pos' , help='set instrument position: p1, p2, p3 or p4 (default %(default)s)') group2.add_argument('--mode', '-m', dest='mode', help=('set instrument mode: standard, shorty_2' ' or mixed (default %(default)s)')) group2.add_argument('--full', '-F', dest='full', action='store_true', help='plot full coverage (default it for the center of the detector)') group2.add_argument('--lint', dest='logscalex', action='store_false', parser.add_argument('--full', '-F', dest='full', action='store_true', help='full detecor coverage (default: only for the center)') parser.add_argument('--lint', dest='logscalex', action='store_false', help='plot t in linear scale') group2.add_argument('--logt', dest='logscalex', action='store_true', parser.add_argument('--logt', dest='logscalex', action='store_true', help='plot t in log scale') group2.add_argument('--linq', dest='logscaley', action='store_false', parser.add_argument('--linq', dest='logscaley', action='store_false', help='plot q in linear scale') group2.add_argument('--logq', dest='logscaley', action='store_true', parser.add_argument('--logq', dest='logscaley', action='store_true', help='plot q in log scale') group2.add_argument('--no-errors', '-e', dest='errors', action='store_false', parser.add_argument('--no-errors', '-e', dest='errors', action='store_false', help='do not plot errors on q/tau averages') group2.add_argument('--no-legend', dest='legend', action='store_false', parser.add_argument('--no-legend', dest='legend', action='store_false', help='do not show plot legend') group2.add_argument('--color-map', '-c', dest='cmap', parser.add_argument('--color-map', '-c', dest='cmap', help='set colormap to use (default=%(default)s)') group2.add_argument('taus', metavar='tau', type=float, nargs='*', help='select list of taus') parser.add_argument('--verbose', '-v', dest='verbose', action='store_true', help='verbose output') parser.add_argument('--version', action='version', Loading scripts/the_cube.py +57 −53 File changed.Preview size limit exceeded, changes collapsed. Show changes Loading
pysen/config.py +0 −1 Original line number Diff line number Diff line Loading @@ -325,7 +325,6 @@ def coverage(lmax, qmin, **kwargs): for j,y in enumerate(ylims): theta1 = get_theta_pix(x,y,L2,sa,ca) theta2, _ = pixel_angle(i,j,theta0) print(theta1, theta2) theta = theta2 lq = get_q(lmax, theta)*lmax q = linspace(get_q(lmax, theta), get_q(lmax-dlam, theta), len(tbins)+1) Loading
pysen/revision.py +4 −4 Original line number Diff line number Diff line Loading @@ -2,9 +2,9 @@ PySEN revision module """ import sys __version__ = "0.7.3" __version__ = "0.7.4" __release__ = "dev1" __date__ = "Feb 13, 2019" __date__ = "Feb 14, 2019" # VERSION = __version__ RELEASE = __release__ Loading @@ -28,8 +28,8 @@ def version(full=False): "get pysen version number" if full: py_ver = "Python %s.%s" % (sys.version_info.major, sys.version_info.minor) return "%s.%s (%s) %s" % (__version__, __release__, __date__, py_ver) return "%s.%s" % (__version__, __release__) return "%s-%s (%s) %s" % (__version__, __release__, __date__, py_ver) return "%s-%s" % (__version__, __release__) # EOF
scripts/nseplot.py 0 → 100755 +476 −0 Original line number Diff line number Diff line #!/usr/bin/env python "Echo plot" import os.path import logging import numpy as np import h5py import matplotlib.pyplot as plt from pysen import ANGSTROM, NANOSECOND, GAUSS, MICRO, OMEGA_N, version from pysen.plotutil import get_nsubplots from pysen.io import echo_to_hdf from pysen.atarilib import fit_echo_current, Spectrum MAX_CHI2_DEFAULT = 1e-2 def plot_field(ax,t,v,**kwargs): "plot mag field" label = kwargs.get('label' , '') max_chi2 = kwargs.get('max_chi2', MAX_CHI2_DEFAULT) v0 = np.sqrt(np.sum(v**2,axis=-1)) # a,b = np.polyfit(t,v0,deg=1) v1 = np.polyval([a,b],t) chi2 = np.sum((v0-v1)**2) cov = np.cov(t,v0) diag = cov[0,0]*cov[1,1] r = cov[1,0]/np.sqrt(diag) if diag>0 else np.nan #r = np.corrcoef(t,v0)[1,0] #maxd = max(abs(v0-v1)) descr1 = r'a=%.2e b=%.3g' % (a, b) descr2 = r'$\rho$=%.3g $\chi^2$=%.1e' % (r, chi2) #descr2 = r'$\rho$=%.3g $\chi^2$=%.1e $(\Delta B)_{max}$=%.3g' % (r, chi2, maxd) # va = np.average(v0) vd = max(abs(v0-va)) vmin = va-max(vd,0.5) vmax = va+max(vd,0.5) ax.plot(t,v0, '.--', lw=1, label=label) ax.plot(t,v1, '-') ax.set_ylim(bottom=vmin, top=vmax) ax.text(0.05,0.20,descr1, transform=ax.transAxes, fontsize=7) ax.text(0.05,0.10,descr2, transform=ax.transAxes, fontsize=7) ax.grid(True, lw=0.5) ax.legend(loc='upper left') if chi2>max_chi2: ax.set_facecolor('yellow') def magnetic_fields_plot(hdfile, iecho=0, **kwargs): "plot magnetic fields" # max_chi2 = kwargs.get('max_chi2', MAX_CHI2_DEFAULT) # base = os.path.splitext(os.path.basename(hdfile.filename))[0] comment = hdfile.attrs['master_comment'] comment = comment.decode("utf-8") sample = comment.split()[0] nph = hdfile['/'].attrs['point_to_down'] for echo in list(hdfile['/data'].values()): if iecho and echo.attrs['id'] != iecho: continue phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND q0 = physics.attrs['hkl'][0] lam0 = physics.attrs['lambda'] cur0 = float(tech.attrs['i5'][0]) phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) bfield = phase['bfield'] cur = phase['phase_current'][:, 0] # actual value dj = np.radians(phasesens*(cur-cur0))/OMEGA_N/lam0 nxplot, nyplot = get_nsubplots(len(bfield)) fig, axes = plt.subplots(nxplot,nyplot, figsize=(9,7)) fig.subplots_adjust(hspace=0.4, wspace=0.4) fig.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ $\tau$=%.3fns' % (sample, base, lam0/ANGSTROM, q0, tau0)) for iplt, bsensor in enumerate(bfield): bmag = bfield[bsensor] label = bsensor.split('.')[-1] ax = axes[iplt//nxplot,iplt%nxplot] plot_field(ax, dj[:nph]/MICRO, bmag[:nph,:]*GAUSS/MICRO, label=label, max_chi2=max_chi2) for iplt in range(iplt+1, nxplot*nyplot): ax = axes[iplt//nxplot,iplt%nxplot] ax.axis('off') def atari_plot(hdfile, iecho=0, **kwargs): "atari plot" # only_echo = kwargs.pop('only_echo', False) # tbin1 = kwargs.pop('tbin1',0) # TOF bins tbin2 = kwargs.pop('tbin2',None) # FIXME: HARDCODED xpix1 = kwargs.pop('xpix1', 10) xpix2 = kwargs.pop('xpix2', 22) ypix1 = kwargs.pop('ypix1', 10) ypix2 = kwargs.pop('ypix2', 22) #ntaus = hdfile.attrs['scan_length'] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nt = hdfile['/detector'].attrs['no_t_channels'] ny = hdfile['/detector'].attrs['no_y_channels'] nx = hdfile['/detector'].attrs['no_x_channels'] nph = n_idx['dn'] comment = hdfile.attrs['master_comment'] for echo in list(hdfile['/data'].values()): if iecho and echo.attrs['id'] != iecho: continue phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND q0 = physics.attrs['hkl'][0] lam0 = physics.attrs['lambda'] phase0 = float(tech.attrs['i5'][0]) try: lmax = params['lambdatable'] except KeyError: lmax = params['lambdaTable'] lmax = np.array(lmax).flatten() dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) det = phase['detector'][...] pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha0 = np.average(pcha) pcha = np.tile(pcha, (ny, nx, nt)).reshape(n_idx['nphases'],ny,nx,nt) pcha = pcha/float(pcha0) det = det/pcha # plt.figure(figsize=(8,8)) # plt.subplot(2,2,1) wlen = np.sum(det, axis=(0,1,2)) plt.step(wlen[tbin1:tbin2], lmax[tbin1:tbin2]/ANGSTROM) plt.ylabel(r'$\lambda$ [$\AA$]') plt.xlim(left=0) plt.grid(True) plt.subplot(2,2,2) xatari = np.sum(det, axis=(1,2)).T extent = (cur[0], cur[nph-1]) if tbin2 is None: extent = extent + (lmax[tbin1]/ANGSTROM, max(lmax)/ANGSTROM) else: extent = extent + (lmax[tbin1]/ANGSTROM, lmax[tbin2]/ANGSTROM) plt.imshow(xatari[tbin1:tbin2,:nph], origin='lower', extent=extent, aspect='auto') plt.subplot(2,2,3) img = np.sum(det[:,:,:,tbin1:tbin2], axis=(0,-1)) plt.imshow(img , aspect='auto') plt.axvline(xpix1, color='k') plt.axvline(xpix2-1, color='k') plt.axhline(ypix1, color='k') plt.axhline(ypix2-1, color='k') plt.subplot(2,2,4) pha = det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2] dn = np.sum(pha[nph:n_idx['up']], axis=(1,2,3)) up = np.sum(pha[n_idx['up']:], axis=(1,2,3)) epha = np.sqrt(np.sum(pha[:nph]**2, axis=(1,2,3))) pha = np.sum(pha[:nph], axis=(1,2,3)) spectrum = Spectrum(lam=lmax, dlam=dlam, flux=wlen) phase_ef, (xef, yef), res = fit_echo_current(cur[:nph], (pha, epha), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) res['up'] = np.average(up),np.sqrt(np.sum(up)/len(up)) res['dn'] = np.average(dn),np.sqrt(np.sum(dn)/len(dn)) #for key in res: # print(key, res[key]) p = plt.errorbar(cur[:nph], pha, yerr=epha, fmt='.', lw=1) dcol = p[0].get_color() plt.plot(xef, yef, '--') plt.axvline(phase0, color='blue', lw=1, ls='--', label='%.3fA' % phase0 ) plt.axvline(phase_ef, color='red', lw=2, ls='-', label='%.3fA' % phase_ef) ax = plt.gca() if not only_echo: R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] ) plt.errorbar(cur[nph:n_idx['up']], dn, fmt='.', color=dcol) plt.errorbar(cur[n_idx['up']:], up, fmt='.', color=dcol) plt.text(0.01, 1.01, r'$R(Q,t)$=%.3g' % R, transform=ax.transAxes) plt.axhline(np.average(up)) plt.axhline(np.average(dn)) plt.xlabel(r'phase current [A]') plt.legend(loc='best') plt.grid(True) plt.suptitle(r'%s $Q$=%.3f$\AA^{-1}$ $\tau$=%.3gns' % (comment.decode("utf-8"), q0,tau0)) def echo_plot(hdfile, iecho=0, **kwargs): "echo plot" # center_only = kwargs.pop('center_only', False) # npix = kwargs.pop('npix', 2) # pix # tbin1 = kwargs.pop('tbin1',0) # TOF bins tbin2 = kwargs.pop('tbin2',None) # min_counts = kwargs.pop('min_counts', 10.0) max_chi2 = kwargs.pop('max_chi2' ,500.0) min_amp = kwargs.pop('min_amp' , 0.20) # FIXME: HARDCODED #xpix1 = kwargs.pop('xpix1', 10) #xpix2 = kwargs.pop('xpix2', 22) #ypix1 = kwargs.pop('ypix1', 10) #ypix2 = kwargs.pop('ypix2', 22) base = os.path.splitext(os.path.basename(hdfile.filename))[0] comment = hdfile.attrs['master_comment'] comment = comment.decode("utf-8") sample = comment.split()[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nt = hdfile['/detector'].attrs['no_t_channels'] ny = hdfile['/detector'].attrs['no_y_channels'] nx = hdfile['/detector'].attrs['no_x_channels'] nph = n_idx['dn'] if center_only: ntaus = 1 for echo in list(hdfile['/data'].values()): if not iecho or echo.attrs['id'] == iecho: ntaus = ntaus + 1 nxtau, nytau = get_nsubplots(ntaus) fig0 = plt.figure(figsize=(8,8)) fig0.suptitle(r'%s | %s' % (sample, base)) itau = 0 resolution = [] for echo in list(hdfile['/data'].values()): if iecho and echo.attrs['id'] != iecho: continue itau = itau + 1 phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND q0 = physics.attrs['hkl'][0] lam0 = physics.attrs['lambda'] phase0 = float(tech.attrs['i5'][0]) try: lmax = params['lambdatable'] except KeyError: lmax = params['lambdaTable'] lmax = np.array(lmax).flatten() dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) det = phase['detector'][...] pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha0 = np.average(pcha) pcha = np.tile(pcha, (ny, nx, nt)).reshape(n_idx['nphases'],ny,nx,nt) pcha = pcha/float(pcha0) det = det/pcha wlen = np.sum(det, axis=(0,1,2)) spectrum = Spectrum(lam=lmax, dlam=dlam, flux=wlen) y = np.sum(det[:,:,:,tbin1:tbin2], axis=(1,2,3)) ey = np.sqrt(y) dn = y[nph:n_idx['up']] up = y[n_idx['up']:] phase_ef0, (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) if center_only: ax = fig0.add_subplot(nxtau, nytau, itau) fig0.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ ' % (sample, base, lam0/ANGSTROM, q0)) ax.errorbar(cur[:nph], y[:nph], yerr=ey[:nph], fmt='k.', lw=2, ms=3) ax.plot(xef, yef, 'r-', label=r'$\tau$=%.3gns' % tau0) ax.axvline(phase_ef0, color='blue', lw=1, ls='--') ax.errorbar(cur[nph:n_idx['up']], dn, fmt='k.', ms=3) ax.errorbar(cur[n_idx['up']: ], up, fmt='k.', ms=3) upave = np.average(up) dnave = np.average(dn) udave = (upave+dnave)/2 res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] ) resolution.append((tau0, R)) ax.legend() ax.grid() #ax.set_xticklabels([]) if itau%nytau!=1: ax.set_yticklabels([]) continue npx, nx, ny, nt = det.shape nyy = ny//npix nxx = nx//npix pha = det.reshape(npx, nyy, npix, nxx, npix, nt) pha = pha.sum(axis=(2, 4)) # pixel rebin pha = np.sum(pha[:,:,:,tbin1:tbin2],axis=-1) # tof summation fig, axes = plt.subplots(nxx,nyy, figsize=(8,8)) fig.subplots_adjust(hspace=0.05, wspace=0.05) ymax = np.amax(pha) ymin = np.amin(pha) for j in range(nyy): for i in range(nxx): ax = axes[i,j] ax.set_xticklabels([]) ax.set_yticklabels([]) y = pha[:,i,j] ey = np.sqrt(pha[:,i,j]) dn = y[nph:n_idx['up']] up = y[n_idx['up']:] # upave = np.average(up) dnave = np.average(dn) udave = (upave+dnave)/2 ave = np.average(y[:nph]) std = np.std(y[:nph]) if udave<min_counts or ave<min_counts: continue if udave<min_counts and std/ave<min_wiggle: continue res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens, phase_pf=phase_ef0) phase_ef, (xef, yef), res = res # if res['chi2']>max_chi2: continue res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) R = 2*res['amplitude'][0]/( res['up'][0] - res['dn'][0] ) if abs(R)<min_amp: continue ax.errorbar(cur[:nph], y[:nph], yerr=ey[:nph], fmt='k.', lw=1, ms=1, label='(%s,%s)' % (i,j)) ax.set_ylim(bottom=ymin, top=ymax) ax.plot(xef, yef, 'r-') ax.axvline(phase_ef, color='blue', lw=1, ls='--') ax.errorbar(cur[nph:n_idx['up']], dn, fmt='k.', ms=1) ax.errorbar(cur[n_idx['up']: ], up, fmt='k.', ms=1) if R<0: ax.set_facecolor('yellow') #ax.legend() # plt.text(0.01, 1.01, r'$R(Q,t)$=%.3g' % R, transform=ax.transAxes) #ax.axhline(np.average(up)) #ax.axhline(np.average(dn)) #plt.xlabel(r'ephase current [A]') #plt.legend(loc='best') #plt.grid(True) fig.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ $\tau$=%.3gns' % (sample, base, lam0/ANGSTROM, q0,tau0)) if center_only: # figr = plt.figure(figsize=(8,8)) # figr.suptitle(r'%s | %s' % (sample, base)) ax = fig0.add_subplot(nxtau,nytau,ntaus) res = np.asarray(resolution) print(res) ax.plot(res[:,0],res[:,1], '.--') ax.set_xscale('log') ax.set_yscale('log') ax.set_ylim(top=1.1) ax.grid(True) def main(): "the main" import argparse plot_action = dict(echo=echo_plot,atari=atari_plot, bfield=magnetic_fields_plot) parser = argparse.ArgumentParser(description='echo plot') parser.set_defaults(loglevel=logging.INFO, outdir='.', plot='echo', tbin1=0, tbin2=None, tau=0) parser.add_argument('file', metavar='filename', help='file to process', nargs='+') parser.add_argument('--plot', '-p', dest='plot', choices=['echo', 'atari', 'bfield']) parser.add_argument('--atari', dest='plot', action='store_const', const='atari') parser.add_argument('--bfield', dest='plot', action='store_const', const='bfield') parser.add_argument('--echo' , dest='plot', action='store_const', const='echo') grp_echo = parser.add_argument_group('echo', 'echo plotting') grp_echo.add_argument('--center-only', '-C', dest='center_only', action='store_true', default=False, help='show only center patch') grp_echo.add_argument('--only-echo', dest='only_echo', action='store_true', default=False, help='show only echo (no up/down)') parser.add_argument('--tau', '-t', dest='tau', type=int, help='set tau to display (default=%(default)s)') parser.add_argument('--t1' , '-b', dest='tbin1', type=int, help='set min TOF bin (default=%(default)s)') parser.add_argument('--t2' , '-B', dest='tbin2', type=int, help='set max TOF bin (default=%(default)s)') parser.add_argument('--save-file', '-s', dest='savefile', help='save figure to a file (do not show)') parser.add_argument('--verbose', '-v', dest='loglevel', action='store_const', const=logging.DEBUG, help='verbose output') parser.add_argument('--quiet', '-q', dest='loglevel', action='store_const', const=logging.WARNING, help='verbose output') parser.add_argument('--version', '-V', action='version', version='%(prog)s pysen={version}'.format(version=version(full=True))) args = parser.parse_args() if args.loglevel>=logging.INFO: log_format=r'%(message)s' else: log_format=r'%(levelname)s: %(funcName)s: %(message)s' logging.basicConfig(level=args.loglevel, format=log_format) log = logging.getLogger() if args.savefile: plt.switch_backend('Agg') kwargs = vars(args) for filename in args.file: basename , ext = os.path.splitext(os.path.basename(filename)) if ext==".echo": log.info('converting %s to HDF', basename) hfile = echo_to_hdf(filename, args.outdir) else: hfile = filename with h5py.File(hfile, 'r') as hdf5file: plot_action[args.plot](hdf5file, iecho=args.tau, **kwargs) if args.savefile: plt.savefig(basename+'-'+args.savefile) if not args.savefile: plt.show() if __name__ == "__main__": main()
scripts/qtau.py +13 −13 Original line number Diff line number Diff line Loading @@ -52,37 +52,37 @@ def main(): group1 = parser.add_argument_group('script', 'input data are read from the script') group1.add_argument('--file', '-f', dest='file', help='read input data from a file') group2 = parser.add_argument_group('command line options', group2 = parser.add_argument_group('shell', 'input data are taken from command line options') group2.add_argument('taus', metavar='tau', type=float, nargs='*', help='select list of taus') group2.add_argument('--lmax', '-l', dest='lmax', type=float, help='set max wavelength, lambda_max (default %(default)s)') group2.add_argument('--qmin', '-q', dest='qmin', type=float, help='set min q (default %(default)s)') group2.add_argument('--tbin', '-t', dest='tbins', type=int, action='append', help='set tbin (default %(default)s)') help='set TOF binning') # (default %(default)s)') group2.add_argument('--pos' , '-p', dest='pos' , help='set instrument position: p1, p2, p3 or p4 (default %(default)s)') group2.add_argument('--mode', '-m', dest='mode', help=('set instrument mode: standard, shorty_2' ' or mixed (default %(default)s)')) group2.add_argument('--full', '-F', dest='full', action='store_true', help='plot full coverage (default it for the center of the detector)') group2.add_argument('--lint', dest='logscalex', action='store_false', parser.add_argument('--full', '-F', dest='full', action='store_true', help='full detecor coverage (default: only for the center)') parser.add_argument('--lint', dest='logscalex', action='store_false', help='plot t in linear scale') group2.add_argument('--logt', dest='logscalex', action='store_true', parser.add_argument('--logt', dest='logscalex', action='store_true', help='plot t in log scale') group2.add_argument('--linq', dest='logscaley', action='store_false', parser.add_argument('--linq', dest='logscaley', action='store_false', help='plot q in linear scale') group2.add_argument('--logq', dest='logscaley', action='store_true', parser.add_argument('--logq', dest='logscaley', action='store_true', help='plot q in log scale') group2.add_argument('--no-errors', '-e', dest='errors', action='store_false', parser.add_argument('--no-errors', '-e', dest='errors', action='store_false', help='do not plot errors on q/tau averages') group2.add_argument('--no-legend', dest='legend', action='store_false', parser.add_argument('--no-legend', dest='legend', action='store_false', help='do not show plot legend') group2.add_argument('--color-map', '-c', dest='cmap', parser.add_argument('--color-map', '-c', dest='cmap', help='set colormap to use (default=%(default)s)') group2.add_argument('taus', metavar='tau', type=float, nargs='*', help='select list of taus') parser.add_argument('--verbose', '-v', dest='verbose', action='store_true', help='verbose output') parser.add_argument('--version', action='version', Loading
scripts/the_cube.py +57 −53 File changed.Preview size limit exceeded, changes collapsed. Show changes