Commit 320f0f3a authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

fixing plot_diffrun

parent 6ded22f9
Loading
Loading
Loading
Loading
+103 −45
Original line number Diff line number Diff line
@@ -12,6 +12,7 @@ import numpy as np
from matplotlib import cm
from matplotlib.patches import Polygon
from matplotlib.colors import LogNorm
from numpy import sqrt
#pylint: disable=no-name-in-module
from scipy.spatial import ConvexHull

@@ -129,7 +130,7 @@ def _plot_counts(ax, x, cnts, label, fmt='ko-', normalize=None):
    "plot counts as a function of phase point"
    if x is None:
        x = range(1,len(cnts)+1)
    ecnts = np.sqrt(cnts*1.0)
    ecnts = sqrt(cnts*1.0)
    if normalize is not None:
        cnts  = cnts *1.0/normalize
        ecnts = ecnts*1.0/normalize
@@ -254,7 +255,7 @@ def _get_q(q, wavelength, t1=0, t2=-1):
    lave   = np.average(wavelength[:,t1:t2])
    dlam   = np.abs(wavelength[:,t1]-wavelength[:,t2])/2.0
    q      = q*wavelength[:,-1]/lave
    dq     = np.sqrt((dlam/lave)**2) #+ (0.5*dtheta/np.tan(theta/2.0))**2)
    dq     = sqrt((dlam/lave)**2) #+ (0.5*dtheta/np.tan(theta/2.0))**2)
    return q, dq


@@ -327,7 +328,7 @@ def plot_diffrun_old(axis, filenames, **kwargs):
            else:
                q, dq = _get_q(q, wlen, t1, t2)

            data  = np.vstack( (q, dq, cnt*fac, np.sqrt(cnt*1.0)*fac) )
            data  = np.vstack( (q, dq, cnt*fac, sqrt(cnt*1.0)*fac) )
            data  = data[:, data[0].argsort()]
            sym   = {'up':'^', 'down':'v'}.get(lbl,'*')
            #
@@ -355,7 +356,6 @@ def plot_diffrun_old(axis, filenames, **kwargs):
            pc0    = result[0::2, 1]
            att0   = result[0::2, 4]*1.0 if use_att else np.ones_like(pc0)
            cnt0   = cnt[0::2]

            #down count
            pc1    = result[1::2, 1]
            att1   = result[1::2, 4]*1.0 if use_att else np.ones_like(pc1)
@@ -364,24 +364,24 @@ def plot_diffrun_old(axis, filenames, **kwargs):
            fac    = (att0[:ndat]*pc1[:ndat])/(att1[:ndat]*pc0[:ndat])
            fratio =  cnt0[:ndat]/cnt1[:ndat]
            eratio = 1.0/cnt0[:ndat] + 1.0/cnt1[:ndat]
            eratio = np.sqrt(eratio)*fratio*fac
            eratio = sqrt(eratio)*fratio*fac
            #
            fup  = att0[:ndat]/pc0[:ndat]*PCHARGE_PER_SECOND
            fdn  = att1[:ndat]/pc1[:ndat]*PCHARGE_PER_SECOND
            xup  = cnt0[:ndat]*fup
            xdn  = cnt1[:ndat]*fdn
            exup = np.sqrt(cnt0[:ndat])*fup
            exdn = np.sqrt(cnt1[:ndat])*fdn
            exup = sqrt(cnt0[:ndat])*fup
            exdn = sqrt(cnt1[:ndat])*fdn
            # up = c + 1/3i
            # dn = 2/3i
            coherent     = xup - 1/2*xdn
            incoherent   = 3/2*xdn
            ecoherent    = np.sqrt(exup**2 + 1/4*exdn**2)
            ecoherent    = sqrt(exup**2 + 1/4*exdn**2)
            eincoherent  = 3/2*exdn


            qmean  = 0.5*(q0[:ndat]+q1[:ndat])
            #dqmean = np.sqrt(dq0[:ndat]**2+dq1[:ndat]**2)
            #dqmean = sqrt(dq0[:ndat]**2+dq1[:ndat]**2)

            if what=='flip_ratio':
                fratio = np.where(fratio>1.0, fratio, 1.0/fratio)
@@ -408,7 +408,7 @@ def plot_diffrun_old(axis, filenames, **kwargs):

            if what=='average':
                average  = (coherent+incoherent)/2
                eaverage = np.sqrt(ecoherent**2+eincoherent**2)
                eaverage = sqrt(ecoherent**2+eincoherent**2)
                data = np.vstack( (qmean, average, eaverage) )
                data = data[:, data[0].argsort()]
                axis.errorbar(data[0], data[1], yerr=data[2],
@@ -449,11 +449,10 @@ def plot_diffrun_old(axis, filenames, **kwargs):
    axis.legend(loc='best')
    axis.grid(True)


def plot_diffrun_new(axis, filenames, **kwargs):
def plot_diffrun(axis, filenames, **kwargs):
    "process files"

    what   = kwargs.pop('what',   ('up', 'down'))
    what   = kwargs.pop('what',   None) or ('up', 'down')
    t1     = kwargs.pop('tmin',   None)
    t2     = kwargs.pop('tmax',   None)
    logscale  = kwargs.pop('log_scale', False)
@@ -473,48 +472,107 @@ def plot_diffrun_new(axis, filenames, **kwargs):

    noqerr = t1 is None or t2 is None
    nfiles = len(filenames)
    ylabel = 'Count Rate [counts/s] at %.1f MW' % power
    axis1  = None

    for _i, filename in enumerate(filenames):
        res  = read_diffrun(filename) #
        meta  = res.get('summary')
        wlen = res.get('wavelength')
        wlen = res.get('wavelength')/ANGSTROM
        det  = res.get('detector')
        #mon1   = res.get('monitor1')
        #mon2   = res.get('monitor2')

        acolor = cmap((_i)/nfiles)
        wlen = wlen/ANGSTROM

        updn = {}

        log.info(filename)
        result = meta
        acolor = cmap((_i)/nfiles)
        ndat   = len(result)//2

        for off, lbl in zip((0, 1), ('up', 'down')):
            q     = meta[off::2, 0] # q_min
            pc    = meta[off::2, 1] # proton charge
        q = result[:, 0] #
        if noqerr:
                cnt   = meta[off::2, 2] # detector sum
            q0, dq0 = q[0::2], None
            q1, dq1 = q[1::2], None
            cnt = result[::, 2] # detector sum
        else:
                cnt   = np.sum(det[off::2,0,t1:t2], axis=-1)
            #win  = result[off::2, 3] # detector win
            att = meta[off::2, 4]*1.0 if use_att else np.ones_like(pc) # attenuator table
            #lmax = result[off::2, 5]/ANGSTROM
            fac   = (att/pc)*PCHARGE_PER_SECOND
            q0, dq0 = _get_q(q[0::2], wlen, t1, t2)
            q1, dq1 = _get_q(q[1::2], wlen, t1, t2)
            cnt = np.sum(det[:,0,t1:t2], axis=-1)

            if noqerr:
                dq   = np.zeros_like(q)
            else:
                q, dq = _get_q(q, wlen, t1, t2)
        # up count
        pc0    = result[0::2, 1]
        att0   = result[0::2, 4]*1.0 if use_att else np.ones_like(pc0)
        cnt0   = cnt[0::2]
        # down count
        pc1    = result[1::2, 1]
        att1   = result[1::2, 4]*1.0 if use_att else np.ones_like(pc1)
        cnt1   = cnt[1::2]
        #
        fup  = att0/pc0*PCHARGE_PER_SECOND
        fdn  = att1/pc1*PCHARGE_PER_SECOND
        cup  = cnt0*fup
        cdn  = cnt1*fdn
        ecup = sqrt(cnt0)*fup
        ecdn = sqrt(cnt1)*fdn

        # flip ratio
        fac    = (att0[:ndat]*pc1[:ndat])/(att1[:ndat]*pc0[:ndat])
        fratio =  cnt0[:ndat]/cnt1[:ndat]
        eratio = 1.0/cnt0[:ndat] + 1.0/cnt1[:ndat]
        eratio = sqrt(eratio)*fratio*fac
        # up = c + 1/3i
        # dn = 2/3i
        coherent     = cup[:ndat] - 1/2*cdn[:ndat]
        incoherent   = 3/2*cdn[:ndat]
        ecoherent    = sqrt(ecup[:ndat]**2 + 1/4*ecdn[:ndat]**2)
        eincoherent  = 3/2*ecdn[:ndat]

        qmean  = 0.5*(q0[:ndat]+q1[:ndat])

            data       = np.vstack( (q, dq, cnt*fac, np.sqrt(cnt*1.0)*fac) )
            updn[off]  = data[:, data[0].argsort()]
            #if out:
            #    np.savetxt(out, data.T, fmt='%.3g', footer=' ',
            #               header = "%s '%s'\nQ counts err" % (filename, lbl))
        def _make_plot(axis, data, **kwargs):
            kwargs.setdefault('ls', '--')
            kwargs.setdefault('lw', 1)
            kwargs.setdefault('color', acolor)
            data = np.vstack( data )
            data = data[:, data[0].argsort()]
            label="%s (%s)" % (os.path.basename(filename), kwargs.pop('lbl'))
            axis.errorbar(data[0], data[1], yerr=data[2], label=label, **kwargs)
            if out:
                np.savetxt(out, data.T, fmt='%.3g', header = "%s\nQ counts err" % label)

        if 'up' in what:
            _make_plot(axis, (q0, cup, ecup ), marker='^', lbl='up')
        if 'down' in what:
            _make_plot(axis, (q1, cdn, ecdn ), marker='v', lbl='down')
        if 'coherent' in what:
            _make_plot(axis, (qmean, coherent, ecoherent ),     marker='^', lbl='C')
        if 'incoherent' in what:
            _make_plot(axis, (qmean, incoherent, eincoherent ), marker='v', lbl='I')
        if 'average' in what:
            average  = (coherent+incoherent)/2
            eaverage = sqrt(ecoherent**2+eincoherent**2)
            _make_plot(axis, (qmean, average, eaverage), marker='s', lbl='<A>')

        # ratios
        if 'flip_ratio' in what:
            #fratio = np.where(fratio>1.0, fratio, 1.0/fratio)
            if axis1 is None: axis1 = axis.twinx()
            _make_plot(axis1, (qmean, fratio, eratio ), marker='o', lbl='FR')
        if 'coherent_ratio' in what:
            cratio  = (2*fratio - 1)/3
            ecratio = (2*eratio)/3
            if axis1 is None: axis1 = axis.twinx()
            _make_plot(axis1, (qmean, cratio, ecratio), marker='s', lbl='C/I')

        #print(data.keys())
    axis.set_yscale('log' if logscale else 'linear')
    if noqerr:
        axis.set_xlabel(r'$Q_{min}$ $\AA^{-1}$')
    else:
        axis.set_xlabel(r'$<Q>$ $\AA^{-1}$')
    
plot_diffrun = plot_diffrun_new
    axis.set_ylabel(ylabel)
    axis.set_title(title)
    axis.legend(loc='best')
    axis.grid(True)


def plot_transmission(ax0, ax1, filenames, **kwargs):
+2 −2
Original line number Diff line number Diff line
@@ -2,8 +2,8 @@
PySEN revision module
"""
import sys
__version__  = "0.71"
__release__  = "dev2"
__version__  = "0.72"
__release__  = "dev1"
__date__     = "Aug 2, 2022"

def version(full=False):
+5 −2
Original line number Diff line number Diff line
@@ -80,10 +80,13 @@ def main():
                       help="plot 'up' counts")
    parser.add_argument('--down',    '-d', dest='what', action='append_const', const='down',
                       help="plot 'down' counts")
    parser.add_argument('--ratio',    '-r', dest='what', action='append_const',  const='flip_ratio',
                       help="plot flip ratio")
    parser.add_argument('--coherent', '-c', dest='what', action='append_const', const='coherent',
                       help="plot coherent and incoherent deconvoluted")
    parser.add_argument('--incoherent', '-i', dest='what', action='append_const', const='incoherent',
                       help="plot coherent and incoherent deconvoluted")
    #
    parser.add_argument('--ratio',    '-r', dest='what', action='append_const',  const='flip_ratio',
                       help="plot flip ratio")
    parser.add_argument('--average',  '-a', dest='what', action='append_const', const='average',
                       help="plot (up+down)/2 sum")
    parser.add_argument('--coherent-ratio', '-C', dest='what', action='append_const', const='coherent_ratio',