Source code for sas.sascalc.fit.qsmearing

"""
    Handle Q smearing
"""
#####################################################################
#This software was developed by the University of Tennessee as part of the
#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
#project funded by the US National Science Foundation.
#See the license text in license.txt
#copyright 2008, University of Tennessee
######################################################################
import math
import logging
import sys

import numpy as np  # type: ignore
from numpy import pi, exp # type:ignore

from sasmodels.resolution import Slit1D, Pinhole1D
from sasmodels.sesans import SesansTransform
from sasmodels.resolution2d import Pinhole2D

from sas.sascalc.data_util.nxsunit import Converter

[docs]def smear_selection(data, model = None): """ Creates the right type of smearer according to the data. The canSAS format has a rule that either slit smearing data OR resolution smearing data is available. For the present purpose, we choose the one that has none-zero data. If both slit and resolution smearing arrays are filled with good data (which should not happen), then we choose the resolution smearing data. :param data: Data1D object :param model: sas.model instance """ # Sanity check. If we are not dealing with a SAS Data1D # object, just return None # This checks for 2D data (does not throw exception because fail is common) if data.__class__.__name__ not in ['Data1D', 'Theory1D']: if data is None: return None elif data.dqx_data is None or data.dqy_data is None: return None return PySmear2D(data) # This checks for 1D data with smearing info in the data itself (again, fail is likely; no exceptions) if not hasattr(data, "dx") and not hasattr(data, "dxl")\ and not hasattr(data, "dxw"): return None # Look for resolution smearing data # This is the code that checks for SESANS data; it looks for the file loader # TODO: change other sanity checks to check for file loader instead of data structure? _found_sesans = False #if data.dx is not None and data.meta_data['loader']=='SESANS': if data.dx is not None and data.isSesans: #if data.dx[0] > 0.0: if np.size(data.dx[data.dx <= 0]) == 0: _found_sesans = True # if data.dx[0] <= 0.0: if np.size(data.dx[data.dx <= 0]) > 0: raise ValueError('one or more of your dx values are negative, please check the data file!') if _found_sesans: # Pre-compute the Hankel matrix (H) SElength = Converter(data._xunit)(data.x, "A") theta_max = Converter("radians")(data.sample.zacceptance)[0] q_max = 2 * np.pi / np.max(data.source.wavelength) * np.sin(theta_max) zaccept = Converter("1/A")(q_max, "1/" + data.source.wavelength_unit), Rmax = 10000000 hankel = SesansTransform(data.x, SElength, data.source.wavelength, zaccept, Rmax) # Then return the actual transform, as if it were a smearing function return PySmear(hankel, model, offset=0) _found_resolution = False if data.dx is not None and len(data.dx) == len(data.x): # Check that we have non-zero data if data.dx[0] > 0.0: _found_resolution = True #print "_found_resolution",_found_resolution #print "data1D.dx[0]",data1D.dx[0],data1D.dxl[0] # If we found resolution smearing data, return a QSmearer if _found_resolution: return pinhole_smear(data, model) # Look for slit smearing data _found_slit = False if data.dxl is not None and len(data.dxl) == len(data.x) \ and data.dxw is not None and len(data.dxw) == len(data.x): # Check that we have non-zero data if data.dxl[0] > 0.0 or data.dxw[0] > 0.0: _found_slit = True # Sanity check: all data should be the same as a function of Q for item in data.dxl: if data.dxl[0] != item: _found_resolution = False break for item in data.dxw: if data.dxw[0] != item: _found_resolution = False break # If we found slit smearing data, return a slit smearer if _found_slit: return slit_smear(data, model) return None
[docs]class PySmear(object): """ Wrapper for pure python sasmodels resolution functions. """ def __init__(self, resolution, model, offset=None): self.model = model self.resolution = resolution if offset is None: offset = np.searchsorted(self.resolution.q_calc, self.resolution.q[0]) self.offset = offset
[docs] def apply(self, iq_in, first_bin=0, last_bin=None): """ Apply the resolution function to the data. Note that this is called with iq_in matching data.x, but with iq_in[first_bin:last_bin] set to theory values for these bins, and the remainder left undefined. The first_bin, last_bin values should be those returned from get_bin_range. The returned value is of the same length as iq_in, with the range first_bin:last_bin set to the resolution smeared values. """ if last_bin is None: last_bin = len(iq_in) start, end = first_bin + self.offset, last_bin + self.offset q_calc = self.resolution.q_calc iq_calc = np.empty_like(q_calc) if start > 0: iq_calc[:start] = self.model.evalDistribution(q_calc[:start]) if end+1 < len(q_calc): iq_calc[end+1:] = self.model.evalDistribution(q_calc[end+1:]) iq_calc[start:end+1] = iq_in[first_bin:last_bin+1] smeared = self.resolution.apply(iq_calc) return smeared
__call__ = apply
[docs] def get_bin_range(self, q_min=None, q_max=None): """ For a given q_min, q_max, find the corresponding indices in the data. Returns first, last. Note that these are indexes into q from the data, not the q_calc needed by the resolution function. Note also that these are the indices, not the range limits. That is, the complete range will be q[first:last+1]. """ q = self.resolution.q first = np.searchsorted(q, q_min) last = np.searchsorted(q, q_max) return first, min(last,len(q)-1)
[docs]def slit_smear(data, model=None): q = data.x width = data.dxw if data.dxw is not None else 0 height = data.dxl if data.dxl is not None else 0 # TODO: width and height seem to be reversed return PySmear(Slit1D(q, height, width), model)
[docs]def pinhole_smear(data, model=None): q = data.x width = data.dx if data.dx is not None else 0 return PySmear(Pinhole1D(q, width), model)
[docs]class PySmear2D(object): """ Q smearing class for SAS 2d pinhole data """ def __init__(self, data=None, model=None): self.data = data self.model = model self.accuracy = 'Low' self.limit = 3.0 self.index = None self.coords = 'polar' self.smearer = True
[docs] def set_accuracy(self, accuracy='Low'): """ Set accuracy. :param accuracy: string """ self.accuracy = accuracy
[docs] def set_smearer(self, smearer=True): """ Set whether or not smearer will be used :param smearer: smear object """ self.smearer = smearer
[docs] def set_data(self, data=None): """ Set data. :param data: DataLoader.Data_info type """ self.data = data
[docs] def set_model(self, model=None): """ Set model. :param model: sas.models instance """ self.model = model
[docs] def set_index(self, index=None): """ Set index. :param index: 1d arrays """ self.index = index
[docs] def get_value(self): """ Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, then find smeared intensity """ if self.smearer: res = Pinhole2D(data=self.data, index=self.index, nsigma=3.0, accuracy=self.accuracy, coords=self.coords) val = self.model.evalDistribution(res.q_calc) return res.apply(val) else: index = self.index if self.index is not None else slice(None) qx_data = self.data.qx_data[index] qy_data = self.data.qy_data[index] q_calc = [qx_data, qy_data] val = self.model.evalDistribution(q_calc) return val