Commit 1c751479 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

diffrun cleanup

parent 320f0f3a
Loading
Loading
Loading
Loading
+18 −197
Original line number Diff line number Diff line
@@ -259,196 +259,6 @@ def _get_q(q, wavelength, t1=0, t2=-1):
    return q, dq


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

    what   = kwargs.pop('what',   None)
    up     = kwargs.pop('up',     True)
    down   = kwargs.pop('down',   True)
    t1     = kwargs.pop('tmin',   None)
    t2     = kwargs.pop('tmax',   None)
    logscale  = kwargs.pop('log_scale', False)
    output    = kwargs.pop('output', None)
    power     = kwargs.pop('power',  1.4)
    use_att   = kwargs.pop('att_table', False)

    log = logging.getLogger()

    selection = []
    title  = 'Diffraction Run Plot'
    out    = None
    lw     = 1
    cmap   = cm.get_cmap('hsv')

    if output:
        out = open(output, 'wb') #pylint: disable=consider-using-with

    if what:
        up=down=False
        ylabel = {  'flip_ratio'        :   'Flip Ratio',
                    'coherent_ratio'    :   'Coherent/Incoherent Ratio',
        }.get(what, 'Rates [arb.]')
    else:
        ylabel = 'Count Rate [counts/s] at %.1f MW' % power
        if up:
            selection.append((0,'up'))
        if down:
            selection.append((1,'down'))

    noqerr = t1 is None or t2 is None
    nfiles = len(filenames)

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

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

        log.info(filename)

        for off, lbl in selection: # zip((0, 1), ('up', 'down')):
            q     = result[off::2, 0] # q_min
            pc    = result[off::2, 1] # proton charge
            if noqerr:
                cnt   = result[off::2, 2] # detector sum
            else:
                cnt   = np.sum(det[off::2,0,t1:t2], axis=-1)
            #win  = result[off::2, 3] # detector win
            att = result[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

            if noqerr:
                dq   = np.zeros_like(q)
            else:
                q, dq = _get_q(q, wlen, t1, t2)

            data  = np.vstack( (q, dq, cnt*fac, sqrt(cnt*1.0)*fac) )
            data  = data[:, data[0].argsort()]
            sym   = {'up':'^', 'down':'v'}.get(lbl,'*')
            #
            axis.errorbar(data[0], data[2], yerr=data[3], #xerr=data[1],
                        marker=sym, ls='--',
                        color=acolor, lw=lw, label="%s (%s)" % (os.path.basename(filename), lbl))
            if out:
                np.savetxt(out, data.T, fmt='%.3g', footer=' ',
                           header = "%s '%s'\nQ counts err" % (filename, lbl))

        if what is not None:
            ndat   = len(result)//2

            q = result[:, 0] #
            if noqerr:
                q0  = q[0::2]
                q1  = q[1::2]
                cnt = result[::, 2] # detector sum
            else:
                q0, _ = _get_q(q[0::2], wlen, t1, t2)
                q1, _ = _get_q(q[1::2], wlen, t1, t2)
                cnt = np.sum(det[:,0,t1:t2], axis=-1)

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

            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
            #
            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 = 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    = sqrt(exup**2 + 1/4*exdn**2)
            eincoherent  = 3/2*exdn


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

            if what=='flip_ratio':
                fratio = np.where(fratio>1.0, fratio, 1.0/fratio)
                data = np.vstack( (qmean, fratio, eratio) )
                data = data[:, data[0].argsort()]
                axis.errorbar(data[0], data[1], yerr=data[2],
                             fmt='o--', lw=lw, color=acolor,
                             label="%s (%s)" % (os.path.basename(filename), 'FR'))
                if out:
                    np.savetxt(out, data.T, fmt='%.3g', footer=' ',
                               header="%s '%s'\nQ dQ FR err" % (filename, 'flip ratio'))

            if what=='coherent_ratio':
                fratio = (2*fratio - 1)/3
                eratio = (2*eratio)/3
                data = np.vstack( (qmean, fratio, eratio) )
                data = data[:, data[0].argsort()]
                axis.errorbar(data[0], data[1], yerr=data[2],
                             fmt='s--', lw=lw, color=acolor,
                             label="%s (%s)" % (os.path.basename(filename), 'C/I'))
                if out:
                    np.savetxt(out, data.T, fmt='%.3g', footer=' ',
                               header="%s '%s'\nQ dQ C/I err" % (filename, 'coherent/incoherent'))

            if what=='average':
                average  = (coherent+incoherent)/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],
                             fmt='s--', lw=lw, color=acolor,
                             label="%s (%s)" % (os.path.basename(filename), '<A>'))
                if out:
                    np.savetxt(out, data.T, fmt='%.3g', footer=' ',
                               header="%s '%s'\nQ dQ A err" % (filename, 'average'))

            if what in ('coherent','incoherent','both'):
                if what in ('coherent','both'):
                    data = np.vstack( (qmean, coherent, ecoherent) )
                    data = data[:, data[0].argsort()]
                    axis.errorbar(data[0], data[1], yerr=data[2],
                             fmt='^--', lw=lw, color=acolor,
                             label="%s (%s)" % (os.path.basename(filename), 'C'))
                    if out:
                        np.savetxt(out, data.T, fmt='%.3g', footer=' ',
                               header="%s '%s'\nQ dQ C err" % (filename, 'coherent'))
                if what in ('incoherent','both'):
                    data = np.vstack( (qmean, incoherent, eincoherent) )
                    data = data[:, data[0].argsort()]
                    axis.errorbar(data[0], data[1], yerr=data[2],
                             fmt='v:', lw=lw, color=acolor,
                             label="%s (%s)" % (os.path.basename(filename), 'I'))
                    if out:
                        np.savetxt(out, data.T, fmt='%.3g', footer=' ',
                               header="%s '%s'\nQ dQ I err" % (filename, 'incoherent'))

    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}$')
    axis.set_ylabel(ylabel)
    axis.set_title(title)
    axis.legend(loc='best')
    axis.grid(True)

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

@@ -472,7 +282,7 @@ def plot_diffrun(axis, filenames, **kwargs):

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

    for _i, filename in enumerate(filenames):
@@ -541,16 +351,21 @@ def plot_diffrun(axis, filenames, **kwargs):

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

        # ratios
        if 'flip_ratio' in what:
@@ -561,18 +376,24 @@ def plot_diffrun(axis, filenames, **kwargs):
            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')
            _make_plot(axis, (qmean, cratio, ecratio), marker='s', lbl='C/I')

    if axis0 is not None:
        axis0.set_yscale('log' if logscale else 'linear')
        axis0.legend(loc='best')
        axis0.grid(True, which='both')
        axis0.set_ylabel('Count Rate [counts/s] at %.1f MW' % power)
    if axis1 is not None:
        axis1.set_yscale('log' if logscale else 'linear')
        axis1.legend(loc='best')
        axis1.grid(True, which='both')
        axis1.set_ylabel('Ratio')

    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}$')
    
    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
@@ -3,8 +3,8 @@ PySEN revision module
"""
import sys
__version__  = "0.72"
__release__  = "dev1"
__date__     = "Aug 2, 2022"
__release__  = "dev2"
__date__     = "Aug 3, 2022"

def version(full=False):
    "get pysen version number"
+0 −57
Original line number Diff line number Diff line
@@ -11,62 +11,6 @@ from pysen import version, setup_logger
from pysen.config import DEFAULT_ACCELERATOR_POWER
from pysen.plot   import plot_diffrun


def main_old():
    "main entry"
    parser = argparse.ArgumentParser()
    parser.set_defaults(up=True, down=True, what=None,
                        log_scale=False, tmin=None, tmax=-1, att_table=None,
                        power=DEFAULT_ACCELERATOR_POWER,
                        plotfile=None, output=None)
    parser.add_argument('file', metavar='filename',
                        help='file to process', nargs='+')
    group = parser.add_mutually_exclusive_group()
    group.add_argument('--up', '-u',    dest='down', action='store_false',
                       help="plot 'up' counts only")
    group.add_argument('--down', '-d',  dest='up', action='store_false',
                       help="plot 'down' counts only")

    group.add_argument('--ratio',    '-r', dest='what', action='store_const', const='flip_ratio',
                       help="plot flip ratio")
    group.add_argument('--coherent', '-c', dest='what', action='store_const', const='coherent',
                       help="plot coherent and incoherent deconvoluted")
    group.add_argument('--average', '-a',  dest='what', action='store_const', const='average',
                       help="plot (coherent+incoherent)/2 sum")
    group.add_argument('--coherent-ratio', '-C', dest='what', action='store_const', const='coherent_ratio',
                       help="plot coherent/incoherent ratio")

    parser.add_argument('--tmin', '-t', dest='tmin', type=int,
                        help='')
    parser.add_argument('--tmax', '-T', dest='tmax', type=int,
                        help='')
    parser.add_argument('--logscale', '-l', dest='log_scale', action='store_true',
                        help="plot in semilogy")
    parser.add_argument('--att', '-A', dest='att_table', action='store_true',
                        help="use attenuation table")

    parser.add_argument('--save', '-s', dest='plotfile')
    parser.add_argument('--out', '-o', dest='output')
    parser.add_argument('--version', action='version',
                        version='%(prog)s pysen={version}'.format(version=version()))


    setup_logger()
    args = parser.parse_args()
    kwargs = vars(args)

    _, axis = plt.subplots(figsize=(8,8))
    try:
        plot_diffrun(axis, args.file, **kwargs)

        if args.plotfile:
            plt.savefig(args.plotfile)
        else:
            plt.show()
    except (OSError, RuntimeError) as e:
        print(e)


def main():
    "main entry"
    parser = argparse.ArgumentParser()
@@ -111,7 +55,6 @@ def main():
    setup_logger()
    args = parser.parse_args()
    kwargs = vars(args)
    print(kwargs)

    _, axis = plt.subplots(figsize=(8,8))
    try: