Changes in sasmodels/resolution.py [bbb887d:2d81cfe] in sasmodels


Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/resolution.py

    rbbb887d r2d81cfe  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
    22 PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC 
    2322 
    2423class Resolution(object): 
    … …  
    6665    *q_calc* is the list of points to calculate, or None if this should 
    6766    be estimated from the *q* and *q_width*. 
    68  
    69     *nsigma* is the width of the resolution function.  Should be 2.5. 
    70     See :func:`pinhole_resolution` for details. 
    71     """ 
    72     def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA): 
     67    """ 
     68    def __init__(self, q, q_width, q_calc=None, nsigma=3): 
    7369        #*min_step* is the minimum point spacing to use when computing the 
    7470        #underlying model.  It should be on the order of 
    … …  
    8682 
    8783        # Protect against models which are not defined for very low q.  Limit 
    88         # the smallest q value evaluated to 0.02*min.  Note that negative q 
    89         # values are trimmed even for broad resolution.  Although not possible 
    90         # from the geometry, they may appear since we are using a truncated 
    91         # gaussian to represent resolution rather than a skew distribution. 
     84        # the smallest q value evaluated (in absolute) to 0.02*min 
    9285        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    93         self.q_calc = self.q_calc[self.q_calc >= cutoff] 
     86        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
    9487 
    9588        # Build weight matrix from calculated q values 
    9689        self.weight_matrix = pinhole_resolution( 
    97             self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION), 
    98             nsigma=nsigma) 
     90            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
     91        self.q_calc = abs(self.q_calc) 
    9992 
    10093    def apply(self, theory): 
    … …  
    108101    *q* points at which the data is measured. 
    109102 
    110     *qx_width* slit width in qx 
    111  
    112     *qy_width* slit height in qy 
     103    *dqx* slit width in qx 
     104 
     105    *dqy* slit height in qy 
    113106 
    114107    *q_calc* is the list of points to calculate, or None if this should 
    … …  
    161154 
    162155 
    163 def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): 
    164     r""" 
     156def pinhole_resolution(q_calc, q, q_width): 
     157    """ 
    165158    Compute the convolution matrix *W* for pinhole resolution 1-D data. 
    166159 
    … …  
    169162    *W*, the resolution smearing can be computed using *dot(W,q)*. 
    170163 
    171     Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution 
    172     function is a broadened triangle, and does not extend over the entire 
    173     range $(-\infty, +\infty)$.  It is important to impose this limitation 
    174     since some models fall so steeply that the weighted value in gaussian 
    175     tails would otherwise dominate the integral. 
    176  
    177164    *q_calc* must be increasing.  *q_width* must be greater than zero. 
    178  
    179     [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects 
    180     in Radially Symmetric Small-Angle Neutron Scattering by Numerical and 
    181     Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14. 
    182     https://doi.org/10.1107/S0021889894010095. 
    183165    """ 
    184166    # The current algorithm is a midpoint rectangle rule.  In the test case, 
    … …  
    188170    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    189171    weights = cdf[1:] - cdf[:-1] 
    190     # Limit q range to +/- 2.5 sigma 
    191     qhigh = q + nsigma*q_width 
    192     #qlow = q - nsigma*q_width  # linear limits 
    193     qlow = q*q/qhigh  # log limits 
    194     weights[q_calc[:, None] < qlow[None, :]] = 0. 
    195     weights[q_calc[:, None] > qhigh[None, :]] = 0. 
    196172    weights /= np.sum(weights, axis=0)[None, :] 
    197173    return weights 
    … …  
    518494 
    519495 
    520 def gaussian(q, q0, dq, nsigma=2.5): 
    521     """ 
    522     Return the truncated Gaussian resolution function. 
     496def gaussian(q, q0, dq): 
     497    """ 
     498    Return the Gaussian resolution function. 
    523499 
    524500    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    525501    """ 
    526     # Calculate the density of the tails; the resulting gaussian needs to be 
    527     # scaled by this amount in ordere to integrate to 1.0 
    528     two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2 
    529     return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) 
     502    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
    530503 
    531504 
    … …  
    585558 
    586559 
    587 def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5): 
     560def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
    588561    """ 
    589562    Romberg integration for pinhole resolution. 
    … …  
    705678        np.testing.assert_equal(output, self.y) 
    706679 
    707     # TODO: turn pinhole/slit demos into tests 
    708  
    709680    def test_pinhole(self): 
    710681        """ 
    … …  
    715686        theory = 12.0-1000.0*resolution.q_calc 
    716687        output = resolution.apply(theory) 
    717         # Note: answer came from output of previous run.  Non-integer 
    718         # values at ends come from the fact that q_calc does not 
    719         # extend beyond q, and so the weights don't balance. 
    720688        answer = [ 
    721             10.47037734, 9.86925860, 
    722             9., 8., 7., 6., 5., 4., 
    723             3.13074140, 2.52962266, 
     689            10.44785079, 9.84991299, 8.98101708, 
     690            7.99906585, 6.99998311, 6.00001689, 
     691            5.00093415, 4.01898292, 3.15008701, 2.55214921, 
    724692            ] 
    725693        np.testing.assert_allclose(output, answer, atol=1e-8) 
    … …  
    764732        self._compare(q, output, answer, 1e-6) 
    765733 
    766     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    767734    def test_pinhole(self): 
    768735        """ 
    … …  
    779746        self._compare(q, output, answer, 3e-4) 
    780747 
    781     @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    782748    def test_pinhole_romberg(self): 
    783749        """ 
    … …  
    795761        #                     2*np.pi/pars['radius']/200) 
    796762        #tol = 0.001 
    797         ## The default 2.5 sigma and no extra points gets 1% 
     763        ## The default 3 sigma and no extra points gets 1% 
    798764        q_calc = None  # type: np.ndarray 
    799765        tol = 0.01 
    … …  
    11141080 
    11151081    if isinstance(resolution, Slit1D): 
    1116         width, height = resolution.qx_width, resolution.qy_width 
     1082        width, height = resolution.dqx, resolution.dqy 
    11171083        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    11181084    else: 
Note: See TracChangeset for help on using the changeset viewer.