Commit cc5343e2 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

mag field reading

parent e2ec5235
Loading
Loading
Loading
Loading
+46 −13
Original line number Diff line number Diff line
@@ -17,16 +17,19 @@ SCHEMA_DEFAULT = {
    'sns-nse': {
        'array.plain'        : ('phase', 'phase_current', 'proton_charge',),
        'array.compressible' : ('counter', 'counter_ccard', 'detector'),
        'array.phase'        : ('Bfield','Temperature',), # 'Pressure'),
        'attribute.drop'     : ('start', 'end', ),
    },
    'sns-nse-detimage': {
        'array.plain'        : (),
        'array.compressible' : (),
        'attribute.drop'     : (),
        'array.phase'        : (),
    },
    'j-nse'  : {
        'array.plain'        : ('phase', ),
        'array.compressible' : ('counters', 'detector'),
        'array.phase'        : (),
        'attribute.drop'     : (),
    },
}
@@ -93,43 +96,73 @@ class HdfConverter:

        arrays_plain        = self._schema.get('array.plain', [])
        arrays_compressible = self._schema.get('array.compressible', [])
        arrays_phase        = self._schema.get('array.phase',[])

        for name, data in echo.items():
            subgroup = grpecho.create_group(name)
            self._log.debug('writing HDF5 echo point %s %s', grpecho.name, name)
            if name == 'phase': # list
                #
                for vname in arrays_plain:
                    value = np.asarray([phasept[vname] for phasept in data])
                    subgroup.create_dataset(vname, data=value)

                #
                for vname in arrays_compressible:
                    value = np.asarray([phasept[vname] for phasept in data])
                    subgroup.create_dataset(vname, data=value,
                                            compression=self.compression, shuffle=shuffle)
                #
                self.write_phaseinfo(data, subgroup, arrays_phase)

            else:
                self.write_pars(subgroup, data)


    def write_pars(self, group, data):
    def write_pars(self, group, data, **kwargs):
        """write "regular" group"""
        attribute_drop = self._schema.get('attribute.drop', [])
        drop_key    = kwargs.get('drop', self._schema.get('attribute.drop', []))
        for key, value in data.items():
            if key in attribute_drop:
            if key in drop_key:
                continue
            keys = key.split(".")
            subgroup, subkey = self.make_subgroup(key, group)
            if isinstance(value, np.ndarray):
                subgroup.create_dataset(subkey, data=value)
            elif key:
                value = encode(value)
                subgroup.attrs[subkey] = value

    def write_phaseinfo(self, data, group, names):
        "write phase info for names"
        phasedict = dict()
        for phasept in data:
            for key in phasept:
                if key.startswith(names):
                    phasedict.setdefault(key,[])
                    phasedict[key].append(phasept[key])

        for vname in phasedict:
            phasegrp, phasekey = self.make_subgroup(vname, group)
            varray = np.asarray(phasedict[vname])
            nph,_  = varray.shape
            varray = varray.reshape(nph,-1,2)
            varrd  = np.asarray(varray[:,:,0], dtype=float)
            varru  = np.asarray(varray[:,:,1])
            dset = phasegrp.create_dataset(phasekey, data=varrd)
            dset.attrs['unit'] = "%s" % varru[0,0]

    @staticmethod
    def make_subgroup(key, group):
        "create subgroup if necessary"
        key  = key.lower()
        keys = key.split(".")
        if len(keys)>1:
            subgroup = group.require_group(keys[0])
            subkey   = ".".join(keys[1:])
        else:
            subgroup = group
            subkey   = key
        return subgroup, subkey

            if isinstance(value, np.ndarray):
                subgroup.create_dataset(subkey, data=value)
            elif key:
                value = encode(value)
                subgroup.attrs[subkey] = value


def echo_to_hdf(filename, outdir, **kwargs):
+2 −2
Original line number Diff line number Diff line
@@ -2,9 +2,9 @@
PySEN revision module
"""
import sys
__version__  = "0.6.1"
__version__  = "0.7.0"
__release__  = "dev1"
__date__     = "Feb 11, 2019"
__date__     = "Feb 13, 2019"
#
VERSION  = __version__
RELEASE  = __release__
+6 −1
Original line number Diff line number Diff line
@@ -52,7 +52,12 @@ def atari_plot(hdfile, iecho=0, **kwargs):
        lam0    = physics.attrs['lambda']
        phase0  = float(tech.attrs['i5'][0])

        lmax    = np.array(params['lambdaTable']).flatten()
        try:
            lmax = params['lambdatable']
        except:
            lmax = params['lambdaTable']
        lmax = np.array(lmax).flatten()

        dlam    = np.ones_like(lmax)*(lmax[1]-lmax[0])

        phasesens = float(params['phaseangle'].attrs['sensitivities'][0])
+6 −1
Original line number Diff line number Diff line
@@ -76,7 +76,12 @@ def echo_plot(hdfile, iecho=0, **kwargs):
        lam0    = physics.attrs['lambda']
        phase0  = float(tech.attrs['i5'][0])

        lmax    = np.array(params['lambdaTable']).flatten()
        try:
            lmax = params['lambdatable']
        except:
            lmax = params['lambdaTable']
        lmax = np.array(lmax).flatten()

        dlam    = np.ones_like(lmax)*(lmax[1]-lmax[0])

        phasesens = float(params['phaseangle'].attrs['sensitivities'][0])

scripts/mag_field.py

0 → 100755
+147 −0
Original line number Diff line number Diff line
#!/usr/bin/env python
"ATARI Plot"

import os.path 
import logging

import numpy as np
import h5py
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

from pysen import ANGSTROM, NANOSECOND, OMEGA_N, MICROTM, version
from pysen.io import echo_to_hdf
from pysen.atarilib import fit_echo_current, Spectrum

GAUSS2UT = 1e2 # Gauss -> micro-Tesla

def bfield_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', 10.0)
    max_chi2   = kwargs.pop('max_chi2'  ,500.0)
    min_amp    = kwargs.pop('min_amp'   ,  0.0)
    # FIXME: HARDCODED
    xpix1   = kwargs.pop('xpix1', 10)
    xpix2   = kwargs.pop('xpix2', 22)
    ypix1   = kwargs.pop('ypix1', 10)
    ypix2   = kwargs.pop('ypix2', 22)

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

    comment = hdfile.attrs['master_comment']
    comment = comment.decode("utf-8")
    sample  = comment.split()[0]

    ntaus   = hdfile.attrs['scan_length']
    nph     = hdfile['/'].attrs['point_to_down']

    results = {}

    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']
        cur0    = float(tech.attrs['i5'][0])
        phasesens = float(params['phaseangle'].attrs['sensitivities'][0])


        bfield = phase['bfield']
        cur    = phase['phase_current'][:, 0] # actual value
        dj     = np.radians(phasesens*(cur-cur0))/OMEGA_N/lam0

        nplots = int(np.ceil(np.sqrt(len(bfield))))

        fig, axes = plt.subplots(nplots,nplots, figsize=(9,6))
        fig.subplots_adjust(hspace=0.3, wspace=0.3)
        fig.suptitle(r'%s | %s | $\lambda$=%.2g$\AA$ $Q$=%.3f$\AA^{-1}$ $\tau$=%.3fns' % (sample, base, lam0/ANGSTROM, q0, tau0))

        for iplt, bsensor in enumerate(bfield):
            bmag = bfield[bsensor]
            ax  = axes[iplt//nplots,iplt%nplots]
            t    = dj[:nph]/MICROTM
            label = bsensor.split('.')[-1]
            v0 = np.sqrt(np.sum(bmag[:nph,:]**2,axis=-1))*GAUSS2UT
            #
            a,b  = np.polyfit(t,v0,deg=1)
            v1   = np.polyval([a,b],t)
            r    = np.corrcoef(t,v0)
            r    = r[1,0]
            chi2 = np.sum((v0-v1)**2)
            maxd = max(abs(v0-v1))
            #
            va   = np.average(v0)
            vd   = max(abs(v0-va))
            vmin = va-max(vd,0.5)
            vmax = va+max(vd,0.5)
            ax.plot(t,v0, '.', label=label)
            ax.plot(t,v1, '-')
            ax.set_ylim(bottom=vmin, top=vmax)
            ax.legend(loc='best')
            if chi2>1e-2:
                ax.set_facecolor('yellow')
        for iplt in range(iplt+1, nplots*nplots):
            ax  = axes[iplt//nplots,iplt%nplots]
            ax.axis('off')


def main():
    "the main"
    import argparse

    parser = argparse.ArgumentParser(description='echo plot')
    parser.set_defaults(loglevel=logging.INFO, outdir='.')
    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('--verbose', '-v', dest='loglevel', action='store_const',
                        const=logging.DEBUG,  help='verbose output')
    parser.add_argument('--quiet', '-q', dest='loglevel', action='store_const',
                        const=logging.WARNING,  help='verbose output')
    parser.add_argument('--version', action='version',
                        version='%(prog)s pysen={version}'.format(version=version(full=True)))
    args = parser.parse_args()

    if args.loglevel>=logging.INFO:
        log_format=r'%(message)s'
    else:
        log_format=r'%(levelname)s: %(funcName)s: %(message)s'

    logging.basicConfig(level=args.loglevel, format=log_format)
    log = logging.getLogger()

    kwargs = vars(args)
    for filename in args.file:
        basename , ext = os.path.splitext(os.path.basename(filename))
        if ext==".echo":
            log.info('converting %s to HDF', basename)
            hfile = echo_to_hdf(filename, args.outdir)
        else:
            hfile = filename
        with h5py.File(hfile, 'r') as hdf5file:
            bfield_plot(hdf5file, iecho=args.echo, **kwargs)
    plt.show()

if __name__ == "__main__":
    main()
Loading