Commit 491e3506 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

added error estimates to xyz

parent 454aae6e
Loading
Loading
Loading
Loading
+12 −8
Original line number Diff line number Diff line
@@ -3,6 +3,7 @@
import logging
import argparse

import numpy as np
import matplotlib.pyplot as plt

from pysen.inout.scans import XYZScan
@@ -21,20 +22,23 @@ def plot_xyz(filename, **kwargs):
    for key, data in scan.data.items():
        rates[key]=0.0
        if data['pcharge']>0:
            rates[key] = data['selection_counts']/data['pcharge']*1e12
            log.info("%s: %-6.6s = %10.2f counts/C", base, key.upper(), rates[key])
            counts = data['selection_counts']
            scale  = data['pcharge']/1e12
            rates[key]        = counts/scale
            rates[key+'_sig'] = np.sqrt(counts)/scale
            log.info("%s: %-6.6s = (%9.1f +/- %9.1f) counts/C", base, key.upper(), rates[key], rates[key+'_sig'])
        #
    res = xyz_decomposition(rates)
    res = xyz_decomposition(rates, errors=True)
    #
    stot = res['sum_ave']
    log.info("%s:", base)
    for key in ('n_coh', 'i_inc', 'm_mag'):
        log.info("%s: %-6.6s = %10.2f counts/C  (%.2f%%)",
            base, key.upper(), res[key], res[key]/stot*100)
        log.info("%s: %-6.6s = (%9.1f +/- %9.1f) counts/C  (%5.2f +/- %4.2f)%%",
            base, key.upper(), res[key], res[key+'_sig'], res[key]/stot*100, res[key+'_sig']/stot*100)
    for key in ('m_up/2', 'm_dn/2',):
        log.info("%s: %-6.6s = %10.2f counts/C",
            base, key.upper(), res[key])
    log.info("%s: TOT    = %10.2f counts/C", base, stot)
        log.info("%s: %-6.6s = (%9.1f +/- %9.1f) counts/C",
            base, key.upper(), res[key], res[key+'_sig'])
    log.info("%s: TOT    = %9.1f counts/C", base, stot)


def main():
+39 −8
Original line number Diff line number Diff line
@@ -15,7 +15,7 @@ from .. import config, get_q, ANGSTROM
from .plotutil import norm_pix, get_pix


def xyz_decomposition(xyzdata):
def xyz_decomposition(xyzdata, errors=False):
    """Decompose {x,y,z}_{up,down} measurements
    into nuclear coherent, incoherent and magnetic components

@@ -27,6 +27,9 @@ def xyz_decomposition(xyzdata):
    sX  = xyzdata['x_up'] + xyzdata['x_dn']
    sY  = xyzdata['y_up'] + xyzdata['y_dn']
    sZ  = xyzdata['z_up'] + xyzdata['z_dn']
    #



    # diff up-down for x, y and z measurements
    dX  = xyzdata['x_up'] - xyzdata['x_dn']
@@ -46,7 +49,32 @@ def xyz_decomposition(xyzdata):
    sM  =  2*dZ - (dX + dY)  # magnetic "cross section"
    sN  = (2*sUP - sDN)/6    # nuclear coherent
    #sI  = sTOT - sN - sM    # incoherent
    sI  = sUP/2 - sM
    sI  = sDN/2 - sM

    if errors:
        dX_sig     = np.sqrt(xyzdata['x_up_sig']**2 +  xyzdata['x_dn_sig']**2)
        dY_sig     = np.sqrt(xyzdata['y_up_sig']**2 +  xyzdata['y_dn_sig']**2)
        dZ_sig     = np.sqrt(xyzdata['z_up_sig']**2 +  xyzdata['z_dn_sig']**2)
        #
        sUP_sig    = np.sqrt(xyzdata['z_up_sig']**2 + xyzdata['x_up_sig']**2 + xyzdata['y_up_sig']**2)
        sDN_sig    = np.sqrt(xyzdata['z_dn_sig']**2 + xyzdata['x_dn_sig']**2 + xyzdata['y_dn_sig']**2)

        sM_UP_sig  = np.sqrt(4*xyzdata['z_up_sig']**2 + xyzdata['x_up_sig']**2 + xyzdata['y_up_sig']**2)
        sM_DN_sig  = np.sqrt(4*xyzdata['z_dn_sig']**2 + xyzdata['x_dn_sig']**2 + xyzdata['y_dn_sig']**2)
        

        sTOT_sig   = np.sqrt(sUP_sig**2 + sDN_sig**2)/3
        sM_sig     = np.sqrt(4*dZ_sig**2 + dY_sig**2 + dX_sig**2)
        sN_sig     = np.sqrt(4*sUP_sig**2 + sDN_sig**2)/6
        sI_sig     = np.sqrt(sDN_sig**2/4 + sM_sig**2)

        #
        res['sum_ave_sig'] = sTOT_sig
        res['n_coh_sig']   = sN_sig
        res['i_inc_sig']   = sI_sig
        res['m_mag_sig']   = sM_sig
        res['m_up/2_sig']  = sM_UP_sig
        res['m_dn/2_sig']  = sM_DN_sig

    #
    res['n_coh'] = sN
@@ -60,6 +88,9 @@ def xyz_decomposition(xyzdata):
    res['m_ave'  ] = sM_UP + sM_DN
    res['m_up/2']  = sM_UP
    res['m_dn/2']  = sM_DN



    return res


+2 −2
Original line number Diff line number Diff line
@@ -3,8 +3,8 @@ PySEN revision module
"""
import sys
__version__  = "2.0"
__release__  = "a6"
__date__     = "Oct 23, 2024"
__release__  = "a7"
__date__     = "Nov 13, 2024"

def version(full=False):
    "get pysen version number"