echoplot reintroduced

#!/usr/bin/env python

"ATARI Plot utils"

import warnings

from functools import partial

from collections import namedtuple

from pysen import OMEGA_N

from .echo.fit import flux_weighted_eshape, echo_fit, echo_curve

Spectrum = namedtuple('Spectrum', ['lam','dlam','flux'])

def fit_poly_current(curr, counts, iphase, deg=4, dix=3):

    imx2 = min(len(cnts), imx+dix)

    xf = curr[imx1:imx2+1]

    yf = cnts[imx1:imx2+1]

    #warnings.simplefilter('ignore', np.RankWarning)

    pf = np.polyfit(xf, yf, deg=deg)

    xf = np.linspace(min(xf), max(xf), 101)

    yf = np.polyval(pf, xf)

    imx = np.argmax(yf)

    return xf[imx], (xf,yf)

def fit_echo_current(current, counts, spectrum, lam0, phase0, phasesens=1.0):

def fit_echo_current(current, counts, spectrum, lam0, phase0, phasesens=1.0, phase_pf=None):

    "fit echo for atari plot and phase tables"

    # get maximum estimate

    if phase_pf is None:

        phase_pf, _ = fit_poly_current(current, counts, None)

    # convert currents to dJ

    dj  = np.radians(phasesens*(current  - phase0))/OMEGA_N/lam0

#!/usr/bin/env python

"ATARI Plot"

"the future echo plot"

from   os.path import basename, splitext

import numpy as np

import h5py

import matplotlib.pyplot as plt

from pysen.io.reader import read_echo

#result = read_echo('./test/s3890.echo')

result = read_echo('./test/s4749.echo')

print(list(result.keys()))

from pysen import ANGSTROM, NANOSECOND

from pysen.atarilib import fit_echo_current, Spectrum

#itau = 0

#ipha = 0

#for key, val in result['data'][itau]['data'][ipha].items():

#    print key, val

#print

##

#for i, line in enumerate(

#    print i, line

def echo_plot(hdfile, iecho=0, **kwargs):

    "echo plot"

    #

    npix    = kwargs.pop('npix', 2)      # pix

    #

    tbin1   = kwargs.pop('tbin1',0)      # TOF bins

    tbin2   = kwargs.pop('tbin2',None)

    #

    min_counts = kwargs.pop('min_counts', 50.0)

    #min_wiggle = kwargs.pop('min_wiggle',  0.1)

    max_chi2   = kwargs.pop('max_chi2'  ,150.0)

    min_amp    = kwargs.pop('min_amp'   ,  0.2)

    # FIXME: HARDCODED

    xpix1   = kwargs.pop('xpix1', 10)

    xpix2   = kwargs.pop('xpix2', 22)

    ypix1   = kwargs.pop('ypix1', 10)

    ypix2   = kwargs.pop('ypix2', 22)

    base    = splitext(basename(hdfile.filename))[0]

    comment = hdfile.attrs['master_comment']

    comment = comment.decode("utf-8")

    #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']

    for echo in list(hdfile['/data'].values()):

        if iecho is not None and echo.attrs['id'] != iecho:

            print("skipping  ", echo.name)

            continue

        else:

            print("processing", echo.name)

            pass

        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])

        lmax    = np.array(params['lambdaTable']).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)

        phase_ef0, _, _ = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum,

                                        lam0=lam0, phase0=phase0, phasesens=phasesens)

        #print(phase_ef0)

        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))

        plt.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]

                y  = pha[:,i,j]

                ey = np.sqrt(pha[:,i,j])

                dn   = y[nph:n_idx['up']]

                up   = y[n_idx['up']:]

                #

                ax.set_xticklabels([])

                ax.set_yticklabels([])

                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

                phase_ef, (xef, yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum,

                                            lam0=lam0, phase0=phase0, phasesens=phasesens, phase_pf=phase_ef0)

                #

                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 (R<min_amp):

                    continue

                #print(i,j,R,res['amplitude'], res['up'], res['dn'])

                #if R<0.2:

                #    continue

                #ax.errorbar(cur, y, yerr=ey, fmt='k.')

                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)

                #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'phase current [A]')

        #plt.legend(loc='best')

        #plt.grid(True)

        plt.suptitle(r'%s | %s | $Q$=%.2f$\AA^{-1}$ $\tau$=%.2gns' % (comment, base, q0,tau0))

def main():

    "the main"

    import argparse

    parser = argparse.ArgumentParser(description='atari plot')

    parser.add_argument('file', metavar='filename', help='file to process', nargs='+')

    parser.add_argument('--echo',  '-e',  dest='echo', type=int,   default=None,

                        help='set echo to display (default=%(default)s)')

    parser.add_argument('--t1',  '-t',  dest='tbin1', type=int,   default=0,

                        help='set min TOF bin (default=%(default)s)')

    parser.add_argument('--t2',  '-T',  dest='tbin2', type=int,   default=None,

                        help='set max TOF bin (default=%(default)s)')

    parser.add_argument('--only-echo', '-N', dest='only_echo', action='store_true',   default=False,

                        help='show only echo (no up/down)')

    args = parser.parse_args()

    kwargs = vars(args)

    for filename in args.file:

        with h5py.File(filename, 'r') as hdf5file:

            echo_plot(hdf5file, iecho=args.echo, **kwargs)

    plt.show()

if __name__ == "__main__":

    main()

        'scripts/qbin.py',

        'scripts/qtau.py',

        'scripts/atariplot.py',

        #'scripts/echoplot.py',

        'scripts/echoplot.py',

        #'scripts/stapler.py',

        #'scripts/the_cube.py',

        #'scripts/webstatus.py',