Loading pysen/echo/__init__.py +3 −3 Original line number Diff line number Diff line Loading @@ -3,5 +3,5 @@ echo fit/plot module """ from .fit import fit_echo_current, Spectrum # NOQA from .reduce import get_symmetry_phase, reduce_echo # NOQA from .reduce import get_symmetry_phase, get_phase_indices # NOQA from .phase_table import make_phase_table # NOQA pysen/echo/phase_table.py +2 −4 Original line number Diff line number Diff line Loading @@ -101,10 +101,8 @@ def make_phase_table(phtab, polyfit=0, ncur=10, nphi=11, pos='p2', threshold=0.2 ph_int = brute_force_interpolation(phase_table, yphi, xi00, theta0, maincur) ph_dif = ph_int-ph_fit ph_rel = 100*ph_dif/ph_fit if abs(ph_dif) > max_abs: max_abs = abs(ph_dif) if abs(ph_rel) > max_rel: max_rel = abs(ph_rel) max_abs = max(max_abs, abs(ph_dif)) max_rel = max(max_rel, abs(ph_rel)) log.info("i00=%9.5f phi=%9.5f, i5_fit=%8.5f i5_tab=%8.5f i5_int=%8.5f diff=(%+.4fA, %+.3f%%)", maincur, theta0, ph_fit, ph_tab, ph_int, ph_dif, ph_rel) log.info("max(diff)=(%.4fA,%.3f%%)", max_abs, max_rel) Loading pysen/echo/reduce.py +189 −183 Original line number Diff line number Diff line #!/usr/bin/env python "reduce echo data" import time import logging import numpy as np from .. import config, get_q, dictionary_hash from ..constants import NANOSECOND, ANGSTROM #from .. import config, get_q, dictionary_hash #from ..constants import NANOSECOND, ANGSTROM from .fit import Spectrum, fit_echo_current def get_phase_indices(hdfile): "get phase indices" idx = {'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } n_phases = idx['dn'] return n_phases, idx def get_symmetry_phase(hdfile, iecho=None, **kwargs): "get the symmetry phase (for phase tables)" # Loading @@ -19,17 +27,14 @@ def get_symmetry_phase(hdfile, iecho=None, **kwargs): ypix1 = kwargs.pop('ypix1', 0) ypix2 = kwargs.pop('ypix2', None) #if kwargs.get('center_only'): # tbin1, tbin2 = 4, -4 # xpix1, xpix2 = 5, -5 # ypix1, ypix2 = 5, -5 log = logging.getLogger() n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, _ = get_phase_indices(hdfile) results = [] for echo in list(hdfile['/data'].values()): Loading @@ -54,22 +59,23 @@ def get_symmetry_phase(hdfile, iecho=None, **kwargs): phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) cur = phase['phase_current'][:, 0] # actual value det = phase['detector'][...] pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha0 = np.average(pcha) pcha = pcha[:,np.newaxis,np.newaxis,np.newaxis]/float(pcha0) pcha = pcha/float(np.average(pcha)) pcha = pcha[:,np.newaxis,np.newaxis,np.newaxis] # det = det/pcha wlen = np.sum(det, axis=(0,1,2)) # all but TOF axis # pha = det[:,ypix1:ypix2,xpix1:xpix2,tbin1:tbin2] epha = np.sqrt(np.sum(pha[:nph], axis=(1,2,3))) pha = np.sum(pha[:nph], axis=(1,2,3)) det = det[:,ypix1:ypix2,xpix1:xpix2,tbin1:tbin2] edet = np.sqrt(np.sum(det[:nph], axis=(1,2,3))) det = np.sum(det[:nph], axis=(1,2,3)) spectrum = Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]) (phase_ef, _), _, _ = fit_echo_current(cur[:nph], (pha, epha), spectrum, (phase_ef, _), _, _ = fit_echo_current(cur[:nph], (det, edet), Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]), lam0=lam0, phase0=phase0, phasesens=phasesens) log.info("i00=%9.5f phi=%9.5f, i5_fit=%8.5f i5_set=%8.5f", maincur, theta0, phase_ef, phase0) Loading @@ -86,170 +92,170 @@ def write_pixel(grp, name, attrs, **kwargs): pixel.create_dataset(key, data=val) def reduce_echo(hdfile, **kwargs): "reduce echo file" # reduction parameters rpar = {} rpar['npix'] = kwargs.get('npix', 4) # pix # default TOF range rpar['tbin1'] = kwargs.get('tbin1', 0) rpar['tbin2'] = kwargs.get('tbin2', None) # default center patch rpar['xpix1'] = kwargs.get('xpix1', 10) rpar['xpix2'] = kwargs.get('xpix2', 22) rpar['ypix1'] = kwargs.get('ypix1', 10) rpar['ypix2'] = kwargs.get('ypix2', 22) #def reduce_echo(hdfile, **kwargs): # "reduce echo file" # # reduction parameters # rpar = {} # rpar['npix'] = kwargs.get('npix', 4) # pix # # default TOF range # rpar['tbin1'] = kwargs.get('tbin1', 0) # rpar['tbin2'] = kwargs.get('tbin2', None) # # default center patch # rpar['xpix1'] = kwargs.get('xpix1', 10) # rpar['xpix2'] = kwargs.get('xpix2', 22) # rpar['ypix1'] = kwargs.get('ypix1', 10) # rpar['ypix2'] = kwargs.get('ypix2', 22) # # # # misc rpar['incoherent'] = kwargs.get('incoherent', False) # # misc # rpar['incoherent'] = kwargs.get('incoherent', False) # # # log = logging.getLogger() # # # 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'] # # # rpar['tbin2'] = rpar['tbin2'] or nt # if kwargs.get('whole_detector'): # rpar['xpix1'], rpar['xpix2'] = 0, nx # rpar['ypix1'], rpar['ypix2'] = 0, ny # # # tbin1, tbin2 = rpar['tbin1'], rpar['tbin2'] # xpix1, xpix2 = rpar['xpix1'], rpar['xpix2'] # ypix1, ypix2 = rpar['ypix1'], rpar['ypix2'] # log = logging.getLogger() # rhash = dictionary_hash(rpar) # reduction_parameters = hdfile.require_group("reduction") # if reduction_parameters.attrs.get('hash') == rhash: # log.info('data already reduced with the same parameters') # return hdfile # 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'] # rpar['tbin2'] = rpar['tbin2'] or nt if kwargs.get('whole_detector'): rpar['xpix1'], rpar['xpix2'] = 0, nx rpar['ypix1'], rpar['ypix2'] = 0, ny # tbin1, tbin2 = rpar['tbin1'], rpar['tbin2'] xpix1, xpix2 = rpar['xpix1'], rpar['xpix2'] ypix1, ypix2 = rpar['ypix1'], rpar['ypix2'] rhash = dictionary_hash(rpar) reduction_parameters = hdfile.require_group("reduction") if reduction_parameters.attrs.get('hash') == rhash: log.info('data already reduced with the same parameters') return hdfile for itau, echo in hdfile['/data'].items(): try: del echo['reduced'] except KeyError: pass reduced_data = echo.require_group('reduced') phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND lam0 = physics.attrs['lambda'] phase0 = float(tech.attrs['i5'][0]) theta0 = float(tech.attrs['mophi'][0]) try: lmax = params['lambdatable'] except KeyError: lmax = params['lambdaTable'] phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) lmax = np.array(lmax).flatten() dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) avlam = lmax - dlam/2 det = phase['detector'][...]*1.0 pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha = pcha/np.average(pcha) wlen = np.sum(det, axis=(0,1,2)) YY, XX, TT = np.mgrid[ypix1:ypix2, xpix1:xpix2, tbin1:tbin2] theta, _ = config.pixel_angle(XX, YY, np.radians(theta0)) # theta,phi alam = avlam[TT] q_pix = get_q(alam, theta)*ANGSTROM tau_pix = (alam/lam0)**3*tau0 flux = np.tile(wlen[tbin1:tbin2], q_pix.shape[:-1]).reshape(q_pix.shape) qave = np.average(q_pix, weights=flux) qvar = np.sqrt(np.average((q_pix-qave)**2, weights=flux)) tave = np.average(tau_pix, weights=flux) tvar = np.sqrt(np.average((tau_pix-tave)**2, weights=flux)) spectrum = Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]) y = np.sum(det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2], axis=(1,2,3)) ey = np.sqrt(y) y = y/pcha ey = ey/pcha # dn = y[nph:n_idx['up']] up = y[n_idx['up']:] (phase_ef0, phase_err) , (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens, incoherent=rpar['incoherent'], global_fit=False) atari = np.sum(det, axis=(1,2)) upave = np.average(up) dnave = np.average(dn) res['q'] = qave, qvar res['tau'] = tave, tvar res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) res['current0'] = (phase_ef0, phase_err) res['tbins'] = (tbin1, tbin2 or -1) res['pixels'] = ((xpix1, ypix1), (xpix2, ypix2)) res['R'] = (2*res['amplitude'][0]/( upave - dnave),0) # TODO error estimate write_pixel(reduced_data, 'center', res, echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), spectrum=spectrum, atari=atari) _, nx, ny, nt = det.shape npix = rpar['npix'] nyy = ny//npix nxx = nx//npix reduced_data.attrs['nxpix'] = nxx reduced_data.attrs['nypix'] = nyy reduced_data.attrs['tbins'] = (tbin1, tbin2) pha = det.reshape(-1, nyy, npix, nxx, npix, nt) pha = pha.sum(axis=(2, 4)) # pixel rebin pha = np.sum(pha[:,:,:,tbin1:tbin2],axis=-1) # tof summation for j in range(nyy): for i in range(nxx): y = pha[:,i,j] ey = np.where(y>0, np.sqrt(y), 1) # no counts has non-zero error bar y = y/pcha ey = ey/pcha # dn = y[nph:n_idx['up']] up = y[n_idx['up']:] upave = np.average(up) dnave = np.average(dn) _ = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase_ef0, phasesens=phasesens, phase_pf=phase_ef0, incoherent=rpar['incoherent']) (phase_ef, phase_err) , (xef, yef), res = _ res['current0'] = phase_ef, phase_err res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) res['tbins'] = tbin1, tbin2 or -1 res['pixels'] = i+1,j+1 res['R'] = 2*res['amplitude'][0]/( upave - dnave),0 # TODO error estimate write_pixel(reduced_data, 'pix-%d-%d' % (i+1,j+1), attrs=res, echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), spectrum=spectrum) log.debug('pixel %d-%d done', i+1,j+1) log.info('tau=%s reduced', itau) # reduction_parameters.attrs['hash'] = rhash reduction_parameters.attrs['date'] = time.ctime() for key, value in rpar.items(): reduction_parameters.attrs[key] = value return hdfile # for itau, echo in hdfile['/data'].items(): # try: # del echo['reduced'] # except KeyError: # pass # reduced_data = echo.require_group('reduced') # # phase = echo['phase'] # physics = echo['phys'] # params = echo['params'] # tech = echo['tech'] # # tau0 = physics.attrs['fouriertime']/NANOSECOND # lam0 = physics.attrs['lambda'] # phase0 = float(tech.attrs['i5'][0]) # theta0 = float(tech.attrs['mophi'][0]) # # try: # lmax = params['lambdatable'] # except KeyError: # lmax = params['lambdaTable'] # phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) # # lmax = np.array(lmax).flatten() # dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) # avlam = lmax - dlam/2 # # det = phase['detector'][...]*1.0 # pcha = phase['proton_charge'][...] # cur = phase['phase_current'][:, 0] # actual value # # pcha = pcha/np.average(pcha) # wlen = np.sum(det, axis=(0,1,2)) # # YY, XX, TT = np.mgrid[ypix1:ypix2, xpix1:xpix2, tbin1:tbin2] # theta, _ = config.pixel_angle(XX, YY, np.radians(theta0)) # theta,phi # alam = avlam[TT] # q_pix = get_q(alam, theta)*ANGSTROM # tau_pix = (alam/lam0)**3*tau0 # # flux = np.tile(wlen[tbin1:tbin2], q_pix.shape[:-1]).reshape(q_pix.shape) # qave = np.average(q_pix, weights=flux) # qvar = np.sqrt(np.average((q_pix-qave)**2, weights=flux)) # tave = np.average(tau_pix, weights=flux) # tvar = np.sqrt(np.average((tau_pix-tave)**2, weights=flux)) # # spectrum = Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]) # y = np.sum(det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2], axis=(1,2,3)) # ey = np.sqrt(y) # y = y/pcha # ey = ey/pcha # # # dn = y[nph:n_idx['up']] # up = y[n_idx['up']:] # # (phase_ef0, phase_err) , (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, # lam0=lam0, phase0=phase0, phasesens=phasesens, # incoherent=rpar['incoherent'], global_fit=False) # atari = np.sum(det, axis=(1,2)) # upave = np.average(up) # dnave = np.average(dn) # res['q'] = qave, qvar # res['tau'] = tave, tvar # res['up'] = upave, np.sqrt(np.sum(up)/len(up)) # res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) # res['current0'] = (phase_ef0, phase_err) # res['tbins'] = (tbin1, tbin2 or -1) # res['pixels'] = ((xpix1, ypix1), (xpix2, ypix2)) # res['R'] = (2*res['amplitude'][0]/( upave - dnave),0) # TODO error estimate # # write_pixel(reduced_data, 'center', res, # echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), # spectrum=spectrum, atari=atari) # # _, nx, ny, nt = det.shape # npix = rpar['npix'] # nyy = ny//npix # nxx = nx//npix # # reduced_data.attrs['nxpix'] = nxx # reduced_data.attrs['nypix'] = nyy # reduced_data.attrs['tbins'] = (tbin1, tbin2) # # pha = det.reshape(-1, nyy, npix, nxx, npix, nt) # pha = pha.sum(axis=(2, 4)) # pixel rebin # pha = np.sum(pha[:,:,:,tbin1:tbin2],axis=-1) # tof summation # # for j in range(nyy): # for i in range(nxx): # y = pha[:,i,j] # ey = np.where(y>0, np.sqrt(y), 1) # no counts has non-zero error bar # y = y/pcha # ey = ey/pcha # # # dn = y[nph:n_idx['up']] # up = y[n_idx['up']:] # # upave = np.average(up) # dnave = np.average(dn) # # _ = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, # lam0=lam0, phase0=phase_ef0, # phasesens=phasesens, phase_pf=phase_ef0, # incoherent=rpar['incoherent']) # (phase_ef, phase_err) , (xef, yef), res = _ # res['current0'] = phase_ef, phase_err # res['up'] = upave, np.sqrt(np.sum(up)/len(up)) # res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) # res['tbins'] = tbin1, tbin2 or -1 # res['pixels'] = i+1,j+1 # res['R'] = 2*res['amplitude'][0]/( upave - dnave),0 # TODO error estimate # write_pixel(reduced_data, 'pix-%d-%d' % (i+1,j+1), attrs=res, # echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), # spectrum=spectrum) # log.debug('pixel %d-%d done', i+1,j+1) # log.info('tau=%s reduced', itau) # # # # reduction_parameters.attrs['hash'] = rhash # reduction_parameters.attrs['date'] = time.ctime() # for key, value in rpar.items(): # reduction_parameters.attrs[key] = value # # return hdfile pysen/plot/echo.py +7 −21 Original line number Diff line number Diff line Loading @@ -9,6 +9,7 @@ import matplotlib.pyplot as plt from matplotlib import colors from ..constants import ANGSTROM, NANOSECOND from ..echo import get_phase_indices from .plotutil import get_nsubplots def plot_single_echo(ax, data, nptr, **kwargs): Loading Loading @@ -61,10 +62,7 @@ def plot_center_echo(hdfile, **kwargs): comment = hdfile.attrs['master_comment'] sample = comment.split()[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, n_idx = get_phase_indices(hdfile) miny, maxy = np.inf, 0 fig0 = plt.figure(figsize=(8,8)) Loading Loading @@ -159,10 +157,7 @@ def plot_echo_portrait(hdfile, iecho=0, **kwargs): comment = hdfile.attrs['master_comment'] sample = comment.split()[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, n_idx = get_phase_indices(hdfile) #miny, maxy = np.inf, 0 itau = 0 Loading Loading @@ -221,10 +216,7 @@ def atari_plot(hdfile, iecho=0, **kwargs): vmin = kwargs.pop('vmin', 0) vmax = kwargs.pop('vmax', None) # n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, n_idx = get_phase_indices(hdfile) log = logging.getLogger() log.info('plotting') Loading Loading @@ -274,17 +266,11 @@ def atari_plot(hdfile, iecho=0, **kwargs): img = np.sum(det[:,:,:,:], axis=(0,-1)) if logz: norm = colors.LogNorm(max(1,vmin), vmax) im = ax.imshow(img , aspect='equal', norm=norm) else: norm = colors.Normalize(vmin, vmax) im = ax.imshow(img , aspect='equal', norm=norm) im = ax.imshow(img , aspect='equal') fig.colorbar(im, ax=ax) #ax.axvline(xpix1, color='gray') #ax.axhline(ypix1, color='gray') #ax.axvline(xpix2-1, color='gray') #ax.axhline(ypix2-1, color='gray') ax = axes[1,1] p = ax.errorbar(cur[:nph], edat[1,:nph], yerr=edat[2,:nph], fmt='.', lw=1) Loading pysen/plot/nseplotlib.py +7 −21 Original line number Diff line number Diff line Loading @@ -12,10 +12,10 @@ from matplotlib import colors from .. import config, get_q from ..constants import ANGSTROM, NANOSECOND, GAUSS, MICRO, OMEGA_N from ..echo import fit_echo_current, Spectrum from ..echo import fit_echo_current, Spectrum, get_phase_indices from ..mathutil import linear_fit from .plotutil import get_nsubplots, norm_pix, get_pix from .echo import plot_single_echo from .echo import plot_single_echo, get_phase_indices MAX_CHI2_DEFAULT_BFIELD = 1e-2 # TODO: explain the meaning of this Loading @@ -29,13 +29,6 @@ def _get_lambda(params): lmax = params['lambdaTable'] return np.asarray(lmax).flatten() #def _get_echo_phase_indices(hfile): # " " # n_idx = { 'dn': hfile['/'].attrs['point_to_down'], # 'up': hfile['/'].attrs['point_to_up'] , # 'nphases': hfile['/'].attrs['no_of_phases'] } # nph = n_idx['dn'] # #def _get_detector_dims(hfile): # " " # nt = hfile['/detector'].attrs['no_t_channels'] Loading Loading @@ -158,14 +151,10 @@ def plot_atari(hdfile, iecho=None, **kwargs): #ntaus = hdfile.attrs['scan_length'] #print(hdfile.filename) base = os.path.splitext(os.path.basename(hdfile.filename))[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } #ctime = hdfile['/'].attrs['preset'] 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'] nph, n_idx = get_phase_indices(hdfile) comment = hdfile.attrs['master_comment'] tbin1, tbin2 = pix['tbin1'], pix['tbin2'] Loading Loading @@ -244,10 +233,10 @@ def plot_atari(hdfile, iecho=None, **kwargs): img = np.sum(det[:,:,:,tbin1:tbin2], axis=(0,-1)) if logz: norm = colors.LogNorm(max(1,vmin), vmax) im = ax.imshow(img , aspect='equal', norm=norm) else: norm = colors.Normalize(vmin, vmax) im = ax.imshow(img , aspect='equal', norm=norm) fig.colorbar(im, ax=ax) ax.axvline(xpix1, color='gray') ax.axhline(ypix1, color='gray') Loading Loading @@ -359,14 +348,11 @@ def plot_map(hdfile, iecho=None, **kwargs): #comment = hdfile.attrs['master_comment'] #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'] #miny, maxy = np.inf, 0 nph, n_idx = get_phase_indices(hdfile) log.warning('maps under development') Loading Loading
pysen/echo/__init__.py +3 −3 Original line number Diff line number Diff line Loading @@ -3,5 +3,5 @@ echo fit/plot module """ from .fit import fit_echo_current, Spectrum # NOQA from .reduce import get_symmetry_phase, reduce_echo # NOQA from .reduce import get_symmetry_phase, get_phase_indices # NOQA from .phase_table import make_phase_table # NOQA
pysen/echo/phase_table.py +2 −4 Original line number Diff line number Diff line Loading @@ -101,10 +101,8 @@ def make_phase_table(phtab, polyfit=0, ncur=10, nphi=11, pos='p2', threshold=0.2 ph_int = brute_force_interpolation(phase_table, yphi, xi00, theta0, maincur) ph_dif = ph_int-ph_fit ph_rel = 100*ph_dif/ph_fit if abs(ph_dif) > max_abs: max_abs = abs(ph_dif) if abs(ph_rel) > max_rel: max_rel = abs(ph_rel) max_abs = max(max_abs, abs(ph_dif)) max_rel = max(max_rel, abs(ph_rel)) log.info("i00=%9.5f phi=%9.5f, i5_fit=%8.5f i5_tab=%8.5f i5_int=%8.5f diff=(%+.4fA, %+.3f%%)", maincur, theta0, ph_fit, ph_tab, ph_int, ph_dif, ph_rel) log.info("max(diff)=(%.4fA,%.3f%%)", max_abs, max_rel) Loading
pysen/echo/reduce.py +189 −183 Original line number Diff line number Diff line #!/usr/bin/env python "reduce echo data" import time import logging import numpy as np from .. import config, get_q, dictionary_hash from ..constants import NANOSECOND, ANGSTROM #from .. import config, get_q, dictionary_hash #from ..constants import NANOSECOND, ANGSTROM from .fit import Spectrum, fit_echo_current def get_phase_indices(hdfile): "get phase indices" idx = {'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } n_phases = idx['dn'] return n_phases, idx def get_symmetry_phase(hdfile, iecho=None, **kwargs): "get the symmetry phase (for phase tables)" # Loading @@ -19,17 +27,14 @@ def get_symmetry_phase(hdfile, iecho=None, **kwargs): ypix1 = kwargs.pop('ypix1', 0) ypix2 = kwargs.pop('ypix2', None) #if kwargs.get('center_only'): # tbin1, tbin2 = 4, -4 # xpix1, xpix2 = 5, -5 # ypix1, ypix2 = 5, -5 log = logging.getLogger() n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, _ = get_phase_indices(hdfile) results = [] for echo in list(hdfile['/data'].values()): Loading @@ -54,22 +59,23 @@ def get_symmetry_phase(hdfile, iecho=None, **kwargs): phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) cur = phase['phase_current'][:, 0] # actual value det = phase['detector'][...] pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha0 = np.average(pcha) pcha = pcha[:,np.newaxis,np.newaxis,np.newaxis]/float(pcha0) pcha = pcha/float(np.average(pcha)) pcha = pcha[:,np.newaxis,np.newaxis,np.newaxis] # det = det/pcha wlen = np.sum(det, axis=(0,1,2)) # all but TOF axis # pha = det[:,ypix1:ypix2,xpix1:xpix2,tbin1:tbin2] epha = np.sqrt(np.sum(pha[:nph], axis=(1,2,3))) pha = np.sum(pha[:nph], axis=(1,2,3)) det = det[:,ypix1:ypix2,xpix1:xpix2,tbin1:tbin2] edet = np.sqrt(np.sum(det[:nph], axis=(1,2,3))) det = np.sum(det[:nph], axis=(1,2,3)) spectrum = Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]) (phase_ef, _), _, _ = fit_echo_current(cur[:nph], (pha, epha), spectrum, (phase_ef, _), _, _ = fit_echo_current(cur[:nph], (det, edet), Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]), lam0=lam0, phase0=phase0, phasesens=phasesens) log.info("i00=%9.5f phi=%9.5f, i5_fit=%8.5f i5_set=%8.5f", maincur, theta0, phase_ef, phase0) Loading @@ -86,170 +92,170 @@ def write_pixel(grp, name, attrs, **kwargs): pixel.create_dataset(key, data=val) def reduce_echo(hdfile, **kwargs): "reduce echo file" # reduction parameters rpar = {} rpar['npix'] = kwargs.get('npix', 4) # pix # default TOF range rpar['tbin1'] = kwargs.get('tbin1', 0) rpar['tbin2'] = kwargs.get('tbin2', None) # default center patch rpar['xpix1'] = kwargs.get('xpix1', 10) rpar['xpix2'] = kwargs.get('xpix2', 22) rpar['ypix1'] = kwargs.get('ypix1', 10) rpar['ypix2'] = kwargs.get('ypix2', 22) #def reduce_echo(hdfile, **kwargs): # "reduce echo file" # # reduction parameters # rpar = {} # rpar['npix'] = kwargs.get('npix', 4) # pix # # default TOF range # rpar['tbin1'] = kwargs.get('tbin1', 0) # rpar['tbin2'] = kwargs.get('tbin2', None) # # default center patch # rpar['xpix1'] = kwargs.get('xpix1', 10) # rpar['xpix2'] = kwargs.get('xpix2', 22) # rpar['ypix1'] = kwargs.get('ypix1', 10) # rpar['ypix2'] = kwargs.get('ypix2', 22) # # # # misc rpar['incoherent'] = kwargs.get('incoherent', False) # # misc # rpar['incoherent'] = kwargs.get('incoherent', False) # # # log = logging.getLogger() # # # 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'] # # # rpar['tbin2'] = rpar['tbin2'] or nt # if kwargs.get('whole_detector'): # rpar['xpix1'], rpar['xpix2'] = 0, nx # rpar['ypix1'], rpar['ypix2'] = 0, ny # # # tbin1, tbin2 = rpar['tbin1'], rpar['tbin2'] # xpix1, xpix2 = rpar['xpix1'], rpar['xpix2'] # ypix1, ypix2 = rpar['ypix1'], rpar['ypix2'] # log = logging.getLogger() # rhash = dictionary_hash(rpar) # reduction_parameters = hdfile.require_group("reduction") # if reduction_parameters.attrs.get('hash') == rhash: # log.info('data already reduced with the same parameters') # return hdfile # 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'] # rpar['tbin2'] = rpar['tbin2'] or nt if kwargs.get('whole_detector'): rpar['xpix1'], rpar['xpix2'] = 0, nx rpar['ypix1'], rpar['ypix2'] = 0, ny # tbin1, tbin2 = rpar['tbin1'], rpar['tbin2'] xpix1, xpix2 = rpar['xpix1'], rpar['xpix2'] ypix1, ypix2 = rpar['ypix1'], rpar['ypix2'] rhash = dictionary_hash(rpar) reduction_parameters = hdfile.require_group("reduction") if reduction_parameters.attrs.get('hash') == rhash: log.info('data already reduced with the same parameters') return hdfile for itau, echo in hdfile['/data'].items(): try: del echo['reduced'] except KeyError: pass reduced_data = echo.require_group('reduced') phase = echo['phase'] physics = echo['phys'] params = echo['params'] tech = echo['tech'] tau0 = physics.attrs['fouriertime']/NANOSECOND lam0 = physics.attrs['lambda'] phase0 = float(tech.attrs['i5'][0]) theta0 = float(tech.attrs['mophi'][0]) try: lmax = params['lambdatable'] except KeyError: lmax = params['lambdaTable'] phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) lmax = np.array(lmax).flatten() dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) avlam = lmax - dlam/2 det = phase['detector'][...]*1.0 pcha = phase['proton_charge'][...] cur = phase['phase_current'][:, 0] # actual value pcha = pcha/np.average(pcha) wlen = np.sum(det, axis=(0,1,2)) YY, XX, TT = np.mgrid[ypix1:ypix2, xpix1:xpix2, tbin1:tbin2] theta, _ = config.pixel_angle(XX, YY, np.radians(theta0)) # theta,phi alam = avlam[TT] q_pix = get_q(alam, theta)*ANGSTROM tau_pix = (alam/lam0)**3*tau0 flux = np.tile(wlen[tbin1:tbin2], q_pix.shape[:-1]).reshape(q_pix.shape) qave = np.average(q_pix, weights=flux) qvar = np.sqrt(np.average((q_pix-qave)**2, weights=flux)) tave = np.average(tau_pix, weights=flux) tvar = np.sqrt(np.average((tau_pix-tave)**2, weights=flux)) spectrum = Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]) y = np.sum(det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2], axis=(1,2,3)) ey = np.sqrt(y) y = y/pcha ey = ey/pcha # dn = y[nph:n_idx['up']] up = y[n_idx['up']:] (phase_ef0, phase_err) , (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens, incoherent=rpar['incoherent'], global_fit=False) atari = np.sum(det, axis=(1,2)) upave = np.average(up) dnave = np.average(dn) res['q'] = qave, qvar res['tau'] = tave, tvar res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) res['current0'] = (phase_ef0, phase_err) res['tbins'] = (tbin1, tbin2 or -1) res['pixels'] = ((xpix1, ypix1), (xpix2, ypix2)) res['R'] = (2*res['amplitude'][0]/( upave - dnave),0) # TODO error estimate write_pixel(reduced_data, 'center', res, echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), spectrum=spectrum, atari=atari) _, nx, ny, nt = det.shape npix = rpar['npix'] nyy = ny//npix nxx = nx//npix reduced_data.attrs['nxpix'] = nxx reduced_data.attrs['nypix'] = nyy reduced_data.attrs['tbins'] = (tbin1, tbin2) pha = det.reshape(-1, nyy, npix, nxx, npix, nt) pha = pha.sum(axis=(2, 4)) # pixel rebin pha = np.sum(pha[:,:,:,tbin1:tbin2],axis=-1) # tof summation for j in range(nyy): for i in range(nxx): y = pha[:,i,j] ey = np.where(y>0, np.sqrt(y), 1) # no counts has non-zero error bar y = y/pcha ey = ey/pcha # dn = y[nph:n_idx['up']] up = y[n_idx['up']:] upave = np.average(up) dnave = np.average(dn) _ = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, lam0=lam0, phase0=phase_ef0, phasesens=phasesens, phase_pf=phase_ef0, incoherent=rpar['incoherent']) (phase_ef, phase_err) , (xef, yef), res = _ res['current0'] = phase_ef, phase_err res['up'] = upave, np.sqrt(np.sum(up)/len(up)) res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) res['tbins'] = tbin1, tbin2 or -1 res['pixels'] = i+1,j+1 res['R'] = 2*res['amplitude'][0]/( upave - dnave),0 # TODO error estimate write_pixel(reduced_data, 'pix-%d-%d' % (i+1,j+1), attrs=res, echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), spectrum=spectrum) log.debug('pixel %d-%d done', i+1,j+1) log.info('tau=%s reduced', itau) # reduction_parameters.attrs['hash'] = rhash reduction_parameters.attrs['date'] = time.ctime() for key, value in rpar.items(): reduction_parameters.attrs[key] = value return hdfile # for itau, echo in hdfile['/data'].items(): # try: # del echo['reduced'] # except KeyError: # pass # reduced_data = echo.require_group('reduced') # # phase = echo['phase'] # physics = echo['phys'] # params = echo['params'] # tech = echo['tech'] # # tau0 = physics.attrs['fouriertime']/NANOSECOND # lam0 = physics.attrs['lambda'] # phase0 = float(tech.attrs['i5'][0]) # theta0 = float(tech.attrs['mophi'][0]) # # try: # lmax = params['lambdatable'] # except KeyError: # lmax = params['lambdaTable'] # phasesens = float(params['phaseangle'].attrs['sensitivities'][0]) # # lmax = np.array(lmax).flatten() # dlam = np.ones_like(lmax)*(lmax[1]-lmax[0]) # avlam = lmax - dlam/2 # # det = phase['detector'][...]*1.0 # pcha = phase['proton_charge'][...] # cur = phase['phase_current'][:, 0] # actual value # # pcha = pcha/np.average(pcha) # wlen = np.sum(det, axis=(0,1,2)) # # YY, XX, TT = np.mgrid[ypix1:ypix2, xpix1:xpix2, tbin1:tbin2] # theta, _ = config.pixel_angle(XX, YY, np.radians(theta0)) # theta,phi # alam = avlam[TT] # q_pix = get_q(alam, theta)*ANGSTROM # tau_pix = (alam/lam0)**3*tau0 # # flux = np.tile(wlen[tbin1:tbin2], q_pix.shape[:-1]).reshape(q_pix.shape) # qave = np.average(q_pix, weights=flux) # qvar = np.sqrt(np.average((q_pix-qave)**2, weights=flux)) # tave = np.average(tau_pix, weights=flux) # tvar = np.sqrt(np.average((tau_pix-tave)**2, weights=flux)) # # spectrum = Spectrum(lam=lmax[tbin1:tbin2], dlam=dlam[tbin1:tbin2], flux=wlen[tbin1:tbin2]) # y = np.sum(det[:,xpix1:xpix2,ypix1:ypix2,tbin1:tbin2], axis=(1,2,3)) # ey = np.sqrt(y) # y = y/pcha # ey = ey/pcha # # # dn = y[nph:n_idx['up']] # up = y[n_idx['up']:] # # (phase_ef0, phase_err) , (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, # lam0=lam0, phase0=phase0, phasesens=phasesens, # incoherent=rpar['incoherent'], global_fit=False) # atari = np.sum(det, axis=(1,2)) # upave = np.average(up) # dnave = np.average(dn) # res['q'] = qave, qvar # res['tau'] = tave, tvar # res['up'] = upave, np.sqrt(np.sum(up)/len(up)) # res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) # res['current0'] = (phase_ef0, phase_err) # res['tbins'] = (tbin1, tbin2 or -1) # res['pixels'] = ((xpix1, ypix1), (xpix2, ypix2)) # res['R'] = (2*res['amplitude'][0]/( upave - dnave),0) # TODO error estimate # # write_pixel(reduced_data, 'center', res, # echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), # spectrum=spectrum, atari=atari) # # _, nx, ny, nt = det.shape # npix = rpar['npix'] # nyy = ny//npix # nxx = nx//npix # # reduced_data.attrs['nxpix'] = nxx # reduced_data.attrs['nypix'] = nyy # reduced_data.attrs['tbins'] = (tbin1, tbin2) # # pha = det.reshape(-1, nyy, npix, nxx, npix, nt) # pha = pha.sum(axis=(2, 4)) # pixel rebin # pha = np.sum(pha[:,:,:,tbin1:tbin2],axis=-1) # tof summation # # for j in range(nyy): # for i in range(nxx): # y = pha[:,i,j] # ey = np.where(y>0, np.sqrt(y), 1) # no counts has non-zero error bar # y = y/pcha # ey = ey/pcha # # # dn = y[nph:n_idx['up']] # up = y[n_idx['up']:] # # upave = np.average(up) # dnave = np.average(dn) # # _ = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum, # lam0=lam0, phase0=phase_ef0, # phasesens=phasesens, phase_pf=phase_ef0, # incoherent=rpar['incoherent']) # (phase_ef, phase_err) , (xef, yef), res = _ # res['current0'] = phase_ef, phase_err # res['up'] = upave, np.sqrt(np.sum(up)/len(up)) # res['dn'] = dnave, np.sqrt(np.sum(dn)/len(dn)) # res['tbins'] = tbin1, tbin2 or -1 # res['pixels'] = i+1,j+1 # res['R'] = 2*res['amplitude'][0]/( upave - dnave),0 # TODO error estimate # write_pixel(reduced_data, 'pix-%d-%d' % (i+1,j+1), attrs=res, # echo=np.vstack((cur, y, ey)), fit=np.vstack((xef, yef)), # spectrum=spectrum) # log.debug('pixel %d-%d done', i+1,j+1) # log.info('tau=%s reduced', itau) # # # # reduction_parameters.attrs['hash'] = rhash # reduction_parameters.attrs['date'] = time.ctime() # for key, value in rpar.items(): # reduction_parameters.attrs[key] = value # # return hdfile
pysen/plot/echo.py +7 −21 Original line number Diff line number Diff line Loading @@ -9,6 +9,7 @@ import matplotlib.pyplot as plt from matplotlib import colors from ..constants import ANGSTROM, NANOSECOND from ..echo import get_phase_indices from .plotutil import get_nsubplots def plot_single_echo(ax, data, nptr, **kwargs): Loading Loading @@ -61,10 +62,7 @@ def plot_center_echo(hdfile, **kwargs): comment = hdfile.attrs['master_comment'] sample = comment.split()[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, n_idx = get_phase_indices(hdfile) miny, maxy = np.inf, 0 fig0 = plt.figure(figsize=(8,8)) Loading Loading @@ -159,10 +157,7 @@ def plot_echo_portrait(hdfile, iecho=0, **kwargs): comment = hdfile.attrs['master_comment'] sample = comment.split()[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, n_idx = get_phase_indices(hdfile) #miny, maxy = np.inf, 0 itau = 0 Loading Loading @@ -221,10 +216,7 @@ def atari_plot(hdfile, iecho=0, **kwargs): vmin = kwargs.pop('vmin', 0) vmax = kwargs.pop('vmax', None) # n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } nph = n_idx['dn'] nph, n_idx = get_phase_indices(hdfile) log = logging.getLogger() log.info('plotting') Loading Loading @@ -274,17 +266,11 @@ def atari_plot(hdfile, iecho=0, **kwargs): img = np.sum(det[:,:,:,:], axis=(0,-1)) if logz: norm = colors.LogNorm(max(1,vmin), vmax) im = ax.imshow(img , aspect='equal', norm=norm) else: norm = colors.Normalize(vmin, vmax) im = ax.imshow(img , aspect='equal', norm=norm) im = ax.imshow(img , aspect='equal') fig.colorbar(im, ax=ax) #ax.axvline(xpix1, color='gray') #ax.axhline(ypix1, color='gray') #ax.axvline(xpix2-1, color='gray') #ax.axhline(ypix2-1, color='gray') ax = axes[1,1] p = ax.errorbar(cur[:nph], edat[1,:nph], yerr=edat[2,:nph], fmt='.', lw=1) Loading
pysen/plot/nseplotlib.py +7 −21 Original line number Diff line number Diff line Loading @@ -12,10 +12,10 @@ from matplotlib import colors from .. import config, get_q from ..constants import ANGSTROM, NANOSECOND, GAUSS, MICRO, OMEGA_N from ..echo import fit_echo_current, Spectrum from ..echo import fit_echo_current, Spectrum, get_phase_indices from ..mathutil import linear_fit from .plotutil import get_nsubplots, norm_pix, get_pix from .echo import plot_single_echo from .echo import plot_single_echo, get_phase_indices MAX_CHI2_DEFAULT_BFIELD = 1e-2 # TODO: explain the meaning of this Loading @@ -29,13 +29,6 @@ def _get_lambda(params): lmax = params['lambdaTable'] return np.asarray(lmax).flatten() #def _get_echo_phase_indices(hfile): # " " # n_idx = { 'dn': hfile['/'].attrs['point_to_down'], # 'up': hfile['/'].attrs['point_to_up'] , # 'nphases': hfile['/'].attrs['no_of_phases'] } # nph = n_idx['dn'] # #def _get_detector_dims(hfile): # " " # nt = hfile['/detector'].attrs['no_t_channels'] Loading Loading @@ -158,14 +151,10 @@ def plot_atari(hdfile, iecho=None, **kwargs): #ntaus = hdfile.attrs['scan_length'] #print(hdfile.filename) base = os.path.splitext(os.path.basename(hdfile.filename))[0] n_idx = { 'dn': hdfile['/'].attrs['point_to_down'], 'up': hdfile['/'].attrs['point_to_up'] , 'nphases': hdfile['/'].attrs['no_of_phases'] } #ctime = hdfile['/'].attrs['preset'] 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'] nph, n_idx = get_phase_indices(hdfile) comment = hdfile.attrs['master_comment'] tbin1, tbin2 = pix['tbin1'], pix['tbin2'] Loading Loading @@ -244,10 +233,10 @@ def plot_atari(hdfile, iecho=None, **kwargs): img = np.sum(det[:,:,:,tbin1:tbin2], axis=(0,-1)) if logz: norm = colors.LogNorm(max(1,vmin), vmax) im = ax.imshow(img , aspect='equal', norm=norm) else: norm = colors.Normalize(vmin, vmax) im = ax.imshow(img , aspect='equal', norm=norm) fig.colorbar(im, ax=ax) ax.axvline(xpix1, color='gray') ax.axhline(ypix1, color='gray') Loading Loading @@ -359,14 +348,11 @@ def plot_map(hdfile, iecho=None, **kwargs): #comment = hdfile.attrs['master_comment'] #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'] #miny, maxy = np.inf, 0 nph, n_idx = get_phase_indices(hdfile) log.warning('maps under development') Loading