Commit 69f3c454 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

added python script to plot q 

after collect type 'plot q', see help plot for more info
parent 7f349d3f

Makefile.version

+1 −1
Original line number Diff line number Diff line
@@ -7,7 +7,7 @@ export PROJARCH=$(PROJECT)-$(VERSION_MAJOR).$(VERSION_MINOR) 


export VERSION_MAJOR=1

export VERSION_MINOR=1

export VERSION_RELEASE=1

export VERSION_RELEASE=2



git_rev=$(shell git rev-parse --short HEAD 2> /dev/null)

ifeq "$(git_rev)" ""

python/drspine/__init__.py

+4 −3
Original line number Diff line number Diff line 
#

__all__     = ['plot_sqt', 'plot_sqt_bin', 'plot_chi2', 'plot_fits', 'plot_maps',

               'plot_fit_pixel', 'plot_fit_map', 'drspine_version']

__all__     = ['plot_sqt', 'plot_sqt_bin', 'plot_chi2', 'plot_q',

               'plot_fits', 'plot_maps', 'plot_fit_pixel', 'plot_fit_map',

               'drspine_version']



from .drspine import plot_sqt , plot_sqt_bin, plot_chi2

from .drspine import plot_sqt , plot_sqt_bin, plot_chi2, plot_q

from .fits    import plot_fits, plot_maps, plot_fit_pixel, plot_fit_map

from .version import __version__ as drspine_version

python/drspine/drspine.py

+47 −43
Original line number Diff line number Diff line
@@ -250,34 +250,38 @@ def plot_sqt(filename, **kwargs): 
        plt.title("Model: %s" % model)









def plot_sqt_bin(filename, full_bw=False, save=False):

    """

    Plot Q,tau binning

    """

# #####################################################################

# plots for _bin.dat files

#

# 0,1,2:  itau, ix, iy

# 3, 4 :  tau_mean, q_mean

# 5, 6 :  sqt, e_sqt

# 7, 8 :  udratio, e_udratio (res)

# 9,10 :  udratio, e_udratio (sig)

    #11,12 :  chisq (res), chisq (sig)



    #s = np.loadtxt(filename, delimiter=',')

#11,12 :  up,   e_up    (sig)

#13,14 :  down, e_down  (sig)

#15,16 :  chisq (res), chisq (sig)

#17    :  run_number

def _load_sqt_bin(filename, full_bw=False):

    s = load_txt(filename, delimiter=',', cleanup=['NaN', 'Infinity'])



    f = s[s[:,0]==0] # only full band width

    t = s[s[:,0]>0]  # all but full band width

    if full_bw:

        return s[s[:,0]==0] # only full band width

    return s[s[:,0]>0]      # all but full band width (TOF slices)



    #centerx = np.logical_or(t[:,1]==4, t[:,1]==5)

    #centery = np.logical_or(t[:,2]==4, t[:,2]==5)

    #x = t[np.logical_and(centerx, centery)]

    plt.figure()



def plot_sqt_bin(filename, save=False, full_bw=False):

    """

    Plot Q,tau binning

    """



    plt.figure()



    if full_bw:

        plt.plot(f[:,3], f[:,4],   'ro', label='Full BW')

    plt.plot(t[:,3], t[:,4],   'k.', label='TOF slices')

        d = _load_sqt_bin(filename, True)

        plt.plot(d[:,3], d[:,4],   'ro', label='Full BW')

    d = _load_sqt_bin(filename, False)

    plt.plot(d[:,3], d[:,4],   'k.', label='TOF slices')



    plt.xscale('log')

    plt.grid(True, which='both')
@@ -291,49 +295,49 @@ def plot_sqt_bin(filename, full_bw=False, save=False): 
def _plot_chi2_histo(bins, chi2, color='r', fmt='.', label=None):

    db   = (bins[1:]-bins[:-1])     # bin width

    b    = (bins[1:]+bins[:-1])/2.0 # bin center

    N=float(sum(chi2))*db[0]

    N=float(sum(chi2))*db

    plt.step(    b,  chi2/N, where='mid', color=color, ls='-')

    plt.errorbar(b,  chi2/N, yerr=np.sqrt(chi2)/N, fmt=color+fmt, label=label)







def _plot_chi2_theory(bins, ndf, color='k', label=None):

def _plot_chi2_theory(bins, ndf, fmt='.', color='k', label=None):

    b    = (bins[1:]+bins[:-1])/2.0 # bin center

    xchi2 =  stats.chi2.pdf(b, ndf)

    plt.plot(b, xchi2, color='k', ls=':', lw=1)



    plt.plot(b, xchi2, color=color, ls=':', lw=1, label=label)





def plot_chi2(filename, full_bw=True, chimax=100.0, nbins=101, ndf=19):

    s = np.loadtxt(filename, delimiter=',')



    f = s[s[:,0]==0] # only full band width

    t = s[s[:,0]>0]  # all but full band width

def plot_chi2(filename, chimax=100.0, nbins=101, ndf=19, full_bw=True):

    d = _load_sqt_bin(filename, full_bw)

    label = ''

    if full_bw:

        label = '(full bw)'



    norm = False

    bins = np.linspace(0, chimax, nbins)



    rchi2, _ = np.histogram(f[:,-2], bins=bins, density=norm)

    schi2, _ = np.histogram(f[:,-1], bins=bins, density=norm)

    rchi3, _ = np.histogram(t[:,-2], bins=bins, density=norm)

    schi3, _ = np.histogram(t[:,-1], bins=bins, density=norm)

    rchi2, _ = np.histogram(d[:,-3], bins=bins, density=norm)

    schi2, _ = np.histogram(d[:,-2], bins=bins, density=norm)



    plt.figure(figsize=(12,9))

    plt.figure()



    if full_bw:

        _plot_chi2_histo(bins, rchi2, color='r', fmt='o', label='resolution (full bw)')

        _plot_chi2_histo(bins, schi2, color='b', fmt='o', label='signal     (full bw)')

    else:

        _plot_chi2_histo(bins, rchi3, color='m', fmt='.', label='resolution (sel bw)')

        _plot_chi2_histo(bins, schi3, color='g', fmt='.', label='signal     (sel bw)')

    _plot_chi2_histo(bins, rchi2, color='r', fmt='o', label='resolution %s' % label)

    _plot_chi2_histo(bins, schi2, color='b', fmt='o', label='signal     %s' % label)

    _plot_chi2_theory(bins, ndf=ndf, color='k', fmt='.', label=r'$\chi^2 (NDF=%d)$' % ndf)



    _plot_chi2_theory(bins, ndf=ndf, color='k', label='NDF=19')

    plt.grid(True, which='both')

    plt.xlabel(r'$\tau$ [ns]')

    plt.ylabel(r'Q [$\AA$]')

    plt.legend(loc='best')





def plot_q(filename, bins=41,  full_bw=False):

    d = _load_sqt_bin(filename, full_bw)



    plt.figure() #figsize=(12,9))



    #plt.xscale('log')

    plt.hist(d[:,4], weights=d[:,11], bins=bins, histtype='step')

    plt.grid(True, which='both')

    plt.xlabel(r'$\tau$ [ns]')

    plt.ylabel(r'Q [$\AA$]')

    plt.legend(loc='best')
 
    plt.xlabel(r'Q [$\AA$]')

    plt.ylabel('Intensity (arb)')

sources/drspine.f90

+16 −4
Original line number Diff line number Diff line
@@ -287,6 +287,7 @@ program drspine 
  istat = add_tab_expansion('plot sqt')

  istat = add_tab_expansion('plot sqt model kww-norm')

  istat = add_tab_expansion('plot sqt_bin')

  istat = add_tab_expansion('plot q')

  istat = add_tab_expansion('plot fits')

  istat = add_tab_expansion('plot maps')

  !
@@ -650,7 +651,7 @@ program drspine 
     !> COMMAND: plot

     !-------------------------------------------------------------

     if(command('plot  ',&

        ' plot [sqt|sqt_bin|fits|maps] [batch] [save]'//LF//&

        ' plot [sqt|sqt_bin|fits|q|maps] [batch] [save]'//LF//&

        '   - interface to ipython/matplotlib plotting (if installed)'//LF//&

        '   sqt  [sub-options] - plots available sqt data'//LF//&

        '        model <name>  - fit model to the data, e.g. kww-norm'//LF//&
@@ -662,6 +663,8 @@ program drspine 
        '        npix <npix>   - plot <npix> number of detector pixel in a row/column'//LF//&

        '        ymin <ymin>   - set min value for y-axis (default is 0)'//LF//&

        '        ymax <ymax>   - set max value for y-axis (default is automatic)'//LF//&

        '   q    [sub-options] - plot q histogram'//LF//&

        '        nbin <nbin>   - plot <nbin> number of histogram bins'//LF//&

        '   maps [sub-options] - plots fit maps on the detector plane'//LF//&

        '        what <map>    - select map to plot (default is sym)'// LF//&

        '        sym           - plot symmetry phase map for each pixel'//LF//&
@@ -678,7 +681,7 @@ program drspine 
        '   save               - option: if present, saves plots as PDF files'//LF//&

        '   batch              - option: if present, system creates PDF plots in batch mode'//LF//&

        ' NOTES:'//LF//&

        '   sqt, sqt_bin - data must be created with the "collect" command (see help collect)'//LF//&

        '   sqt, sqt_bin,q  - data must be created with the "collect" command (see help collect)'//LF//&

        '   fits, maps      - data must be saved with the "dump" command (see help dump)$')) then

        !               ===============

        call cmd_plot()
@@ -2092,7 +2095,7 @@ CONTAINS 
    ! arguments to parse

    integer :: irun                          ! select run to plot (fits/maps)

    integer :: itau                          ! select tau to plot (fits/maps)

    integer :: npix                          ! pixels per row to plot (fits/maps)

    integer :: npix                          ! pixels per row to plot (fits/maps) or bins (q)

    real(kind=DBL) :: ymin, ymax             ! y/z-scale limits (fits/maps)

    character(len=MAX_LINE_LENGTH) :: cmodel ! model to fit

    character(len=MAX_LINE_LENGTH) :: cmap   ! which map to plot (default 'sym' - symmetry phase)
@@ -2151,6 +2154,15 @@ CONTAINS 
       extra_args = "npix="//trim(ctmp)//", "//trim(extra_args)

       write(ctmp,'(i0)') itau

       extra_args = "tau="//trim(ctmp)//","//trim(extra_args)

    else if ( found('q') ) then

       pyfunc='plot_q'

       filename = program_param%file_sqt_bin

       extra_args=' '

       if (found('nbins') ) then

           npix = get_named_value('nbins ', npix, inew)

           write(ctmp,'(i0)') npix

           extra_args = "bins="//trim(ctmp)//", "//trim(extra_args)

       endif

    else if ( found('maps') ) then

       pyfunc='plot_maps'

       write(filename, "(a,i0)") trim(program_param%file_fitpref), irun