Commit 2f633c0b authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

cleanup, unit testing, phases parsing

parent 12122166

pysen/__init__.py

+1 −1
Original line number Diff line number Diff line
@@ -5,5 +5,5 @@ Main PySEN package 
from .revision  import version                                             # NOQA

from .constants import *

from .physics   import (get_theta, get_q, q_binning, l_binning, t_binning, # NOQA

                        get_tau)                                           # NOQA

                        get_tau, tof2lambda, lambda2tof)                   # NOQA

from .misc      import setup_logger, dictionary_hash, DEFAULT_LOG_LEVEL    # NOQA

pysen/inout/nexus.py

+36 −43
Original line number Diff line number Diff line
@@ -25,7 +25,7 @@ import datetime as dtm 
import numpy as np

import h5py



from pysen import HMN, MICRO, ANGSTROM

from pysen import MICRO, ANGSTROM, tof2lambda



TDC_SHIFT =    10

TDC_MAXCH = 0x3FF
@@ -36,12 +36,6 @@ Y0_PIX = 518 
XWID2     =   813//2 # 410 # half-width

YWID2     =   821//2 # 410



TAU_IDX   = 1000

DN_IDX    =  100

UP_IDX    =  200



TOF2LAM = MICRO*HMN/ANGSTROM # conversion factor for TOF [us] to lambda [A]



def pixel_decompose(pixel):

    "helper function to decompose pixel into x and y coordinates"

    ypix   = (pixel >> TDC_SHIFT) & TDC_MAXCH
@@ -60,6 +54,18 @@ def timevalue_array(nexus_file, vname, prefix="/entry/DASlogs", dtype=TimeFloatV 
    v_value = nexus_file[f"{prefix}/{vname}/value"][:]

    return np.array(list(zip(v_time,v_value)), dtype=dtype)



def get_average(val, begin, end):

    """

    calculate average between two time stamps

    """

    sel = np.logical_and(begin<val['time'],   val['time']<end)

    sel = val[sel]

    if sel.size:

        avg = np.average(sel['value'])

    else:

        avg = np.average(val['value'])

    return avg



def time_format(ctime, fmt='%c'):

    "convert 'nexus' time to 'echo' time format"

    ctime = ctime.split('.') # remove nanoseconds
@@ -67,8 +73,8 @@ def time_format(ctime, fmt='%c'): 
    ptime = dtm.datetime.strptime(ctime, '%Y-%m-%dT%H:%M:%S%z')

    return  ptime.strftime(fmt) #Fri Oct 24 05:26:54 2014



def get_echo_info(nxsfile, *extra_keys):

    "get a base run info as a dictionary"

def get_run_info(nxsfile):

    "get a base SNS-NSE run info as a dictionary"

    base = os.path.basename(nxsfile.filename)

    log = logging.getLogger('nexus')
@@ -95,12 +101,12 @@ def get_echo_info(nxsfile, *extra_keys): 
                    base, key, default_value)

            value = default_value

        return value



    # FIXME: put these in some confg list

    res['mo_l']  = get_average_value('/entry/DASlogs/BL15:Mot:mo_l',  21.3)

    res['mophi'] = get_average_value('/entry/DASlogs/BL15:Mot:mophi', None)

    res['qmin' ] = get_average_value('/entry/DASlogs/BL15:Mot:Qmin', 0.1)

    res['lmin' ] = get_average_value('/entry/DASlogs/BL15:Chop:Skf1:WavelengthSet', None)

    res['lmax' ] = get_average_value('/entry/DASlogs/BL15:Chop:Skf4:WavelengthSet', None)

    res['qmin' ] = get_average_value('/entry/DASlogs/BL15:Mot:Qmin', 0.1)

    #

    res['phasesens'] = get_average_value('/entry/DASlogs/BL15:Legacy:PhaseAngleSensitivity',1585.0)

    res['sample_temp'] = get_average_value('/entry/DASlogs/BL15:SE:LS1:SAMPLE_TEMP', 303.0)
@@ -114,32 +120,14 @@ def get_echo_info(nxsfile, *extra_keys): 
        except KeyError:

            res[key] = None



    # extra keys

    for key in extra_keys:

        res[key] = get_average_value(key, None)





    return res





def get_tau_index(scan_idx):

    "return tau, phase index"

    idx   = scan_idx['value']

    return (idx // TAU_IDX, idx %  TAU_IDX)



def get_averages(values, begin, end):

    "get averages"

    res = {}

    if values is None:

        return res

    for key, val in values.items():

        sel = np.logical_and(begin<val['time'],   val['time']<end)

        sel = val[sel]

        if len(sel):

            avg = np.average(sel['value'])

        else:

            avg = np.average(val['value'])

        res[key] = avg

def print_nexus_info(runinfo):

    """print nexus info"""

    res = ""

    with np.printoptions(precision=3, threshold=3, linewidth=132, suppress=True):

        for k, v in runinfo.items():

            res = res + f"\n{k:16}: {v}"

    return res



def process_scan_point(events, pcharge, ltot, begin=None, end=None, values=None):
@@ -154,7 +142,8 @@ def process_scan_point(events, pcharge, ltot, begin=None, end=None, values=None) 
    res['pcharge'] = np.sum(pcharge[sel]['value'])



    # other stuff

    res.update(get_averages(values, begin, end))

    for key, val in values.items():

        res[key] = get_average(val, begin, end)



    # events selection

    ev_sel = np.logical_and(begin<ev_t0,ev_t0<end)
@@ -162,7 +151,7 @@ def process_scan_point(events, pcharge, ltot, begin=None, end=None, values=None) 
    id1, id2 = evid[0], evid[-1] # first and last event



    # event_time_offset is in usec

    lam = events['event_time_offset'][id1:id2]*TOF2LAM/ltot

    lam = tof2lambda(events['event_time_offset'][id1:id2]*MICRO, ltot)/ANGSTROM

    pix = pixel_decompose(events['event_id'][id1:id2])

    res['lambda']  = lam

    res['pixel']   = pix
@@ -175,20 +164,24 @@ def dummy_scanpoint_callback(*_args, **_kwargs): 
    return True



# process scan

def process_scan(filename, scan_point_callback=dummy_scanpoint_callback, **kwargs):

def process_nexus(filename, scan_point_callback=dummy_scanpoint_callback, **kwargs):

    "process nexus file (scan indices)"

    log = logging.getLogger('nexus')

    base = os.path.basename(filename)

    nxsfile = h5py.File(filename, 'r')



    # get file info

    info    = get_run_info(nxsfile)



    # get the scan indices

    scan_index     = timevalue_array(nxsfile, 'scan_index', dtype=TimeIntValueDType)

    scan_index_len = len(scan_index)

    if scan_index_len <= 1: # need at least two indices

        log.warning("%s: no subscans found", base)

        return None



    info    = get_echo_info(nxsfile)

    #

    info['scan_indices'] = scan_index

    #

    pcharge = timevalue_array(nxsfile, 'proton_charge')

    events  = nxsfile['/entry/bank1_events']

    #

pysen/inout/nxs2echo.py

+22 −9
Original line number Diff line number Diff line
@@ -3,6 +3,7 @@ 
Piotr Zolnierczuk, zolnierczukp AT ornl.gov

Oak Ridge National Lab, 2024

"""

import os.path

import logging

import argparse
@@ -11,23 +12,35 @@ from pysen.inout import EchoScan, make_echofilename 
def main():

    "the main"

    parser = argparse.ArgumentParser()

    parser.set_defaults(ntof=42, npix=32, nphases=19, outdir='.')

    parser.set_defaults(ntof=42, npix=32, step=45.0, outdir='.', loglevel=logging.INFO)

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

    parser.add_argument('--npix', type=int, dest='npix',

                        help='set number of detector pixels (default %(default)s)')

    parser.add_argument('--ntof', type=int, dest='ntof',

                        help='set number of tof channels (default %(default)s)')

    parser.add_argument('--nphases', type=int, dest='nphases',

                        help='set number of phases (default %(default)s)')

    parser.add_argument('--outdir', '-o', metavar='out', dest='outdir',

                        help='set output directory (default %(default)s)')

    parser.add_argument('--npix', '-P', type=int, dest='npix',

                        help='set number of detector pixels (default %(default)s)')

    parser.add_argument('--ntof', '-T', type=int, dest='ntof',

                        help='set number of tof channels (default %(default)s)')

    parser.add_argument('--step', '-S', type=float, dest='step',

                        help='set precession phase step (default %(default)s)')

    #

    parser.add_argument('--verbose', '-v', dest='loglevel', action='store_const',

                        const=logging.DEBUG,   help='increase verbosity level')

    parser.add_argument('--quiet'  , '-q', dest='loglevel', action='store_const',

                        const=logging.WARNING, help='suppress info messages')



    args = parser.parse_args()



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



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

    log = logging.getLogger()



    for file_name in args.filename:

        if not os.path.exists(file_name):

            log.warning("file '%s' does not exist", file_name)

            continue

        echo = EchoScan(make_echofilename(file_name))

        if not echo.read_nexus(file_name, npix=args.npix, ntof=args.ntof, nphases=args.nphases):

        if not echo.read_nexus(file_name,

            npix=args.npix, ntof=args.ntof, phase_step=args.step):

            continue

        echo.to_echo(args.outdir)

pysen/inout/writer.py

+72 −20
Original line number Diff line number Diff line 
#!/usr/bin/env python

"""Simple nexus to echo converter

"""Simple echo writer

Piotr Zolnierczuk, zolnierczukp AT ornl.gov

Oak Ridge National Lab, 2024

"""
@@ -12,29 +12,34 @@ import numpy as np 
from ..config import (  XDET, YDET, NXCHAN as NPIX, NTCHAN as NTOF,

                        DEFAULT_PROTON_ENERGY as PROTON_ENERGY )

from ..constants import PICO, NANO, MICRO, ANGSTROM, GAUSS

from .nexus   import process_scan, get_tau_index, X0_PIX, Y0_PIX, XWID2, YWID2

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

PHASE_INDEX = 100



PHASE_UP    = 2 # up phase padding

PHASE_DN    = 1 # down phase padding

PHASES_NPOL = 2+3 # 2xPC for down and 3xPC for up

#



B_SENSORS = (('sample'  , 'b_sample'),

             ('fpi21'   , 'b_encl1' ),

             ('fpi22'   , 'b_encl2' ),

             ('external', 'b_ext1'),

             ('external', 'b_ext2'),

             ('external', 'b_ext3'),

             ('external', 'b_ext4'),

             ('external', 'b_ext5'))

#             name        type for echo file

B_SENSORS = (('b_sample', 'sample'  ),

             ('b_encl1' , 'fpi21'   ),

             ('b_encl2' , 'fpi22'   ),

             ('b_ext1'  , 'external'),

             ('b_ext2'  , 'external'),

             ('b_ext3'  , 'external'),

             ('b_ext4'  , 'external'),

             ('b_ext5'  , 'external'))



MOTORS = ('mophi', 'mopsi', 'mo_z', 'moana', 'mo_l', 'moatt')



DEFAULT_PHASE_STEP = 45.0 # degrees

DEFAULT_NUM_PHASES = 19



TAU_IDX   = 1000

DN_IDX    =  100

UP_IDX    =  200





def make_echofilename(file_name):

    "convert nexus filename to .echo filename"
@@ -50,6 +55,11 @@ def get_phase_index(scan_index, offset=0): 
        return updn + offset

    return scan_index%PHASE_INDEX - offset - 1



def get_tau_index(scan_idx):

    "return tau, phase index"

    idx   = scan_idx['value']

    return (idx // TAU_IDX, idx %  TAU_IDX)





class EchoScan:

    "Convert EPICS/ADARA/NeXus file to .echo file"
@@ -78,17 +88,54 @@ class EchoScan: 
        self.tau[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"

        """read nexus file"""

        self.log.info("reading nexus file %s", nxsfile)

        #

        kwargs.setdefault('phase_step', DEFAULT_PHASE_STEP)

        kwargs.setdefault('nphases',    DEFAULT_NUM_PHASES)

        kwargs.setdefault('npix', NPIX)

        kwargs.setdefault('ntof', NTOF)

        kwargs.setdefault('phase_step', DEFAULT_PHASE_STEP)

        #

        self.info = process_scan(nxsfile, self.phasepoint_callback)

        self.info = process_nexus(nxsfile, self.phasepoint_callback)

        if self.info is None:

            return False

        if 'echo' not in self.info.get('notes','').lower():
@@ -96,14 +143,19 @@ class EchoScan: 
            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





    def scan_info(self):

        "scan info section"

        echobase = os.path.splitext(self.echofile)[0]
@@ -177,9 +229,9 @@ class EchoScan: 
        out.write(f"phase {phasedeg:g}\n")

        out.write(f"phase_current {phasecur:.6f} {phasecur:.6f}\n")

        #

        for blbl, bsens in B_SENSORS:

        for bsens, btyp in B_SENSORS:

            b = res[f"{bsens}"]*MICRO/GAUSS

            out.write(f"Bfield {blbl:8s} {bsens:8s} {b:.6f} G 0.0 G 0.0 G\n")

            out.write(f"Bfield {btyp:8s} {bsens:8s} {b:.6f} G 0.0 G 0.0 G\n")

        # proton charge

        out.write(f"proton_charge {res['pcharge']:.0f}\n")

        # counter data (ignored)
@@ -239,7 +291,7 @@ class EchoScan: 
                fd.write(f"*  i5 {i5:14.6f} {i5:14.6f}\n")

                #

                fd.write( "c bsensors (x,y,z) aimed | actual\n")

                for _, bsens in B_SENSORS:

                for bsens, _ in B_SENSORS:

                    bfld = res[bsens]*MICRO/GAUSS

                    fd.write(f"* {bsens:8s} 7.7 7.7 7.7 {bfld:14.6f} 0.0 0.0\n")

                #

pysen/physics/__init__.py

+1 −0
Original line number Diff line number Diff line
@@ -4,3 +4,4 @@ 


from .elastic  import get_theta, get_q, get_qvec, get_dspace  # NOQA

from .spinecho import get_tau,q_binning, l_binning, t_binning # NOQA

from .misc     import tof2lambda, lambda2tof                  # NOQA