Commit 65fb972c authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

getting ready for "direct nexus access"

parent 1e360b56
Loading
Loading
Loading
Loading

misc/nxs_xyz.py

0 → 100755
+87 −0
Original line number Diff line number Diff line
#!/usr/bin/env python
"Simple Atari Plot for NeXus SNS-NSE Files"
import logging
import argparse

import matplotlib.pyplot as plt

from pysen.inout.scans import XYZScan
from pysen.plot.xyzplotlib import xyz_decomposition, plot_xyz_nexus

class XYZ_Plot(XYZScan):
    "XYZ PLOT"
    def xyz_analysis(self):
        """ xyz analysis"""
        rates = {}

        for key, data in self.xyz.items():
            rates[key]=0.0
            if data['pcharge']>0:
                rates[key] = data['selection_counts']/data['pcharge']*1e12
                self.log.info("%s: %-6.6s = %10.2f counts/C", self.info['base'], key.upper(), rates[key])
        #
        res = xyz_decomposition(rates)
        #
        stot = res['sum_ave']
        self.log.info("%s:", self.info['base'])
        for key in ('n_coh', 'i_inc', 'm_mag'):
            self.log.info("%s: %-6.6s = %10.2f counts/C  (%.2f%%)",
                self.info['base'], key.upper(), res[key], res[key]/stot*100)
        for key in ('m_up/2', 'm_dn/2',):
            self.log.info("%s: %-6.6s = %10.2f counts/C",
                self.info['base'], key.upper(), res[key])
        self.log.info("%s: TOT    = %10.2f counts/C", self.info['base'], stot)
        #

    def plot(self, **kwargs):
        "plot xyz"
        plot_xyz_nexus(self.xyz, self.info, **kwargs)
        self.xyz_analysis()


def main():
    "the main"
    parser = argparse.ArgumentParser()
    parser.set_defaults(ntof=42, npix=32,
                        xpix1=10,xpix2=-10,ypix1=10,ypix2=-10, whole_detector=False,
                        show=True, save=False)
    parser.add_argument('filename', metavar='file', nargs='+', help='filename to plot')

    parser.add_argument('--x1' , dest='xpix1', type=int,
                     help='set min. X pix (default=%(default)s)')
    parser.add_argument('--x2' , dest='xpix2', type=int,
                     help='set max. X pix (default=%(default)s)')
    parser.add_argument('--y1' , dest='ypix1', type=int,
                     help='set min. Y pix (default=%(default)s)')
    parser.add_argument('--y2' , dest='ypix2', type=int,
                     help='set max. Y pix (default=%(default)s)')
    parser.add_argument('--whole-detector', '-W', dest='whole_detector', action='store_true',
                        help='use the whole detector')
    #
    parser.add_argument('--npix' , dest='npix', type=int,
                     help='set number of pixels (default=%(default)s)')
    parser.add_argument('--ntof' , dest='ntof', type=int,
                     help='set number of pixels (default=%(default)s)')
    #
    parser.add_argument('--no-show', '-N', dest='show',  action='store_false',
                        help='do not show plots')
    parser.add_argument('--save',    '-S', dest='save',  action='store_true',
                        help='save figure to a file')
    args = parser.parse_args()


    logging.basicConfig(format='%(levelname)s %(message)s', level=logging.INFO)

    kwargs = vars(args)
    filenames = kwargs.pop('filename')
    for file_name in filenames:
        scan = XYZ_Plot()
        scan.read_nexus(file_name, **kwargs)
        scan.plot(**kwargs) #vmax=100)

    if args.show:
        plt.show()


if __name__ == "__main__":
    main()
+1 −1
Original line number Diff line number Diff line
@@ -5,6 +5,6 @@ from .hdf import HdfConverter, convert_to_hdf # NOQA
from .reader  import (  read_echo, read_magnetic, read_xyz,         # NOQA
                        read_detimage, read_datfile, read_diffrun,  # NOQA
                        read_datreat , read_transmission )          # NOQA
from .writer  import EchoScan, make_echofilename
from .writer  import EchoWriter, make_echofilename
# from .legacy_writer  import write_csv, generate_echo                       # NOQA
# from .legacy_history import extract_data, list_variables, parse_date       # NOQA
+9 −5
Original line number Diff line number Diff line
@@ -158,25 +158,29 @@ def process_scan_point(events, pcharge, ltot, begin=None, end=None, values=None)
    res['neutron'] = np.vstack((lam,pix))
    return res

def dummy_scanpoint_callback(*_args, **_kwargs):
def _dummy_callback(*_args, **_kwargs):
    """dummy process scan point (for testing)
    a caller to process scan should override it"""
    return True

# process scan
def process_nexus(filename, scan_point_callback=dummy_scanpoint_callback, **kwargs):
def process_nexus(filename, **kwargs):
    "process nexus file (scan indices)"
    scan_type    = kwargs.pop('scan_type', None)
    scanpoint_cb = kwargs.pop('scanpoint_cb', _dummy_callback)
    #
    log = logging.getLogger('nexus')
    base = os.path.basename(filename)
    nxsfile = h5py.File(filename, 'r')

    # get file info
    info    = get_run_info(nxsfile)
    if scan_type and scan_type not in info.get('notes','').lower():
        log.warning('not an %s scan (%s)', scan_type, info['notes'])
        return None
    if info['proton_charge']<=0.0:
        log.warning("%s: empty run (pcharge=%g)", base, info['proton_charge'])
        return None
        

    # get the scan indices
    scan_index     = timevalue_array(nxsfile, 'scan_index', dtype=TimeIntValueDType)
    scan_index_len = len(scan_index)
@@ -214,7 +218,7 @@ def process_nexus(filename, scan_point_callback=dummy_scanpoint_callback, **kwar
                    it =  info[key]['time']<=scan_index[i]['time'] # tau set before scan
                    res[key] = info[key][it]['value'][-1]

                if not scan_point_callback(res, filename=base,
                if not scanpoint_cb(res, filename=base,
                                begin=scan_index[i], end=scan_index[k], **kwargs):
                    break
        else:
+2 −2
Original line number Diff line number Diff line
@@ -7,7 +7,7 @@ import os.path
import logging
import argparse

from pysen.inout import EchoScan, make_echofilename
from pysen.inout import EchoWriter, make_echofilename

def main():
    "the main"
@@ -38,7 +38,7 @@ def main():
        if not os.path.exists(file_name):
            log.warning("file '%s' does not exist", file_name)
            continue
        echo = EchoScan(make_echofilename(file_name))
        echo = EchoWriter(make_echofilename(file_name))
        if not echo.read_nexus(file_name,
            npix=args.npix, ntof=args.ntof, phase_step=args.step):
            continue

pysen/inout/scans.py

0 → 100644
+163 −0
Original line number Diff line number Diff line
#!/usr/bin/env python
"""Simple echo writer
Piotr Zolnierczuk, zolnierczukp AT ornl.gov
Oak Ridge National Lab, 2024
"""
import logging

import numpy as np

from ..config import NXCHAN as NPIX, NTCHAN as NTOF
from .nexus   import process_nexus, print_nexus_info, X0_PIX, Y0_PIX, XWID2, YWID2

# BUNCH OF HARD CODED CONSTANTS
# TO DO: some of this is repeated in nexus.py file (merge!)
DEFAULT_PHASE_STEP = 45.0 # degrees

TAU_IDX   = 1000
DN_IDX    =  100
UP_IDX    =  200


def get_tau_index(scan_idx):
    "return tau, phase index"
    idx   = scan_idx['value']
    return (idx // TAU_IDX, idx %  TAU_IDX)


class EchoScan:
    """Read SNS-NSE EPICS/NeXus echo scan file"""
    def __init__(self):
        "the constructor"
        self.log  = logging.getLogger()
        self.data = {}
        self.info = None

    def phasepoint_callback(self, data, **kwargs):
        "phase point callback"
        filename = kwargs.get('filename', '<unset>')
        begin = kwargs.get('begin')
        end   = kwargs.get('end')
        itau, iphase = get_tau_index(begin)
        self.data.setdefault(itau, [])
        if not self.data[itau]:
            self.log.info("%s: processing tau=%2d", filename, itau)
        self.log.debug("tau:phase=%2d:%3d index=%4d (%10.3f - %10.3f)   %10.3eC %10d",
                        itau, iphase, begin['value'], begin['time'], end['time'],
                        data['pcharge'], data['phase'])
        data['iphase'] = iphase
        self.data[itau].append(data)
        return True

    def get_taus_and_phases(self):
        """Get number of taus and phases"""
        indices = self.info['scan_indices']
        ctau = 0
        taus = []
        res  = []
        for idx in sorted(set(indices[indices['value']>0]['value'])):
            itau = idx // TAU_IDX
            ipha = idx %  TAU_IDX
            if itau != ctau:
                if res:
                    taus.append((ctau, res))
                ctau = itau
                res = []
            res.append(ipha)

        if res:
            taus.append((ctau, res))

        nphases = 0
        ntaus   = len(taus)
        for t,p in taus:
            has_dn, has_up = False, False
            if DN_IDX in p:
                has_dn = True
                p.remove(DN_IDX)
            if UP_IDX in p:
                has_up = True
                p.remove(UP_IDX)
            if not (has_dn and has_up):
                self.log.warning('%s: up and/or down phase missing for tau=%d', self.info['base'], t)
            if not nphases:
                nphases=len(p)
            if nphases!=len(p):
                self.log.warning('%s: tau %s has different number of phases %d (previously %d)',
                    self.info['base'], t, len(p), nphases)
        return ntaus, nphases


    def read_nexus(self, nxsfile, **kwargs):
        """read nexus file"""
        self.log.info("reading nexus file %s", nxsfile)
        #
        kwargs.setdefault('npix', NPIX)
        kwargs.setdefault('ntof', NTOF)
        kwargs.setdefault('phase_step', DEFAULT_PHASE_STEP)
        #
        self.info = process_nexus(nxsfile, scanpoint_cb=self.phasepoint_callback, scan_type='echo')
        if self.info is None:
            return False
        self.info.update(**kwargs)

        ntaus, nphases = self.get_taus_and_phases()
        self.log.info('%s: found %s taus, with %d phases in each tau', self.info['base'], ntaus, nphases)

        lmin = self.info.get('lmin', 5.0)
        lmax = self.info.get('lmax', 8.0)
        self.info['xbin'] = np.linspace(X0_PIX-XWID2, X0_PIX+XWID2, self.info['npix']+1)
        self.info['ybin'] = np.linspace(Y0_PIX-YWID2, Y0_PIX+YWID2, self.info['npix']+1)
        self.info['tbin'] = np.linspace(lmin, lmax, self.info['ntof']+1)
        self.info['nphases'] = nphases
        # debug
        self.log.debug("nexus info: %s", print_nexus_info(self.info))
        return True

class XYZScan:
    """Read SNS-NSE EPICS/NeXus xyz scan file"""
    def __init__(self):
        "the constructor"
        self.log  = logging.getLogger()
        self.xyz  =  {}
        self.info = None

    def phasepoint_callback(self, data, **kwargs):
        "phase point callback"
        filename = kwargs.get('filename', '<unset>')
        begin = kwargs.get('begin')
        xyz   = { 10: 'z_up',
                  11: 'z_dn',
                  20: 'y_up',
                  21: 'y_dn',
                  30: 'x_up',
                  31: 'x_dn',}.get(begin['value'])

        rate=0.0
        if data['pcharge']>0:
            rate = data['neutron'].shape[1]/data['pcharge']*1e12
        self.xyz.setdefault(xyz, data)
        self.log.debug("%s: processing xyz=%s, pc=%.5f, rate=%10.1f", filename, xyz, data['pcharge']*1e-12, rate)
        return True

    def read_nexus(self, nxsfile, **kwargs):
        """read nexus file"""
        self.log.info("reading nexus file %s", nxsfile)
        #
        kwargs.setdefault('npix', NPIX)
        kwargs.setdefault('ntof', NTOF)
        #
        self.info = process_nexus(nxsfile, scanpoint_cb=self.phasepoint_callback, scan_type='xyz')
        if self.info is None:
            return False
        self.info.update(**kwargs)
        #

        lmin = self.info.get('lmin', 5.0)
        lmax = self.info.get('lmax', 8.0)
        self.info['xbin'] = np.linspace(X0_PIX-XWID2, X0_PIX+XWID2, self.info['npix']+1)
        self.info['ybin'] = np.linspace(Y0_PIX-YWID2, Y0_PIX+YWID2, self.info['npix']+1)
        self.info['tbin'] = np.linspace(lmin, lmax, self.info['ntof']+1)
        self.log.debug("nexus info: %s", print_nexus_info(self.info))
        return True
# EOF
Loading