source: sasview/DataLoader/qsmearing.py @ ca10d8e

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since ca10d8e was c0d9981, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

preparing for release

  • Property mode set to 100644
File size: 10.8 KB
Line 
1"""
2This software was developed by the University of Tennessee as part of the
3Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4project funded by the US National Science Foundation.
5
6See the license text in license.txt
7
8copyright 2008, University of Tennessee
9"""
10
11
12import numpy
13import math
14import scipy.special
15
16def smear_selection(data1D):
17    """
18        Creates the right type of smearer according
19        to the data.
20   
21        The canSAS format has a rule that either
22        slit smearing data OR resolution smearing data
23        is available.
24       
25        For the present purpose, we choose the one that
26        has none-zero data. If both slit and resolution
27        smearing arrays are filled with good data
28        (which should not happen), then we choose the
29        resolution smearing data.
30       
31        @param data1D: Data1D object
32    """
33    # Sanity check. If we are not dealing with a SANS Data1D
34    # object, just return None
35    if data1D.__class__.__name__ != 'Data1D' \
36        or not hasattr(data1D, "dxl") or not hasattr(data1D, "dxw"):
37        return None
38   
39    # Look for resolution smearing data
40    _found_resolution = False
41    if data1D.dx is not None and len(data1D.dx)==len(data1D.x):
42       
43        # Check that we have non-zero data
44        if data1D.dx[0]>0.0:
45            _found_resolution = True
46       
47    # If we found resolution smearing data, return a QSmearer
48    if _found_resolution == True:
49        return QSmearer(data1D)
50
51    # Look for slit smearing data
52    _found_slit = False
53    if data1D.dxl is not None and len(data1D.dxl)==len(data1D.x) \
54        and data1D.dxw is not None and len(data1D.dxw)==len(data1D.x):
55       
56        # Check that we have non-zero data
57        if data1D.dxl[0]>0.0 or data1D.dxw[0]>0.0:
58            _found_slit = True
59       
60        # Sanity check: all data should be the same as a function of Q
61        for item in data1D.dxl:
62            if data1D.dxl[0] != item:
63                _found_resolution = False
64                break
65           
66        for item in data1D.dxw:
67            if data1D.dxw[0] != item:
68                _found_resolution = False
69                break
70           
71    # If we found slit smearing data, return a slit smearer
72    if _found_slit == True:
73        return SlitSmearer(data1D)
74   
75    return None
76           
77
78class _BaseSmearer(object):
79   
80    def __init__(self):
81        self.nbins = 0
82        self._weights = None
83       
84    def _compute_matrix(self): return NotImplemented
85
86    def __call__(self, iq):
87        """
88            Return the smeared I(q) value at the given q.
89            The smeared I(q) is computed using a predetermined
90            smearing matrix for a particular binning.
91       
92            @param q: I(q) array
93            @return: smeared I(q)
94        """
95        # Sanity check
96        if len(iq) != self.nbins:
97            raise RuntimeError, "Invalid I(q) vector: inconsistent array length %s <> %s" % (len(iq), self.nbins)
98       
99        if self._weights == None:
100            self._compute_matrix()
101           
102        iq_smeared = numpy.zeros(self.nbins)
103        # Loop over q-values
104        for q_i in range(self.nbins):
105            sum = 0.0
106            counts = 0.0   
107           
108            for i in range(self.nbins):
109                sum += iq[i] * self._weights[q_i][i]
110                counts += self._weights[q_i][i]
111               
112            iq_smeared[q_i] = sum/counts
113               
114        return iq_smeared   
115   
116class _SlitSmearer(_BaseSmearer):
117    """
118        Slit smearing for I(q) array
119    """
120   
121    def __init__(self, nbins=None, width=None, height=None, min=None, max=None):
122        """
123            Initialization
124           
125            @param iq: I(q) array [cm-1]
126            @param width: slit width [A-1]
127            @param height: slit height [A-1]
128            @param min: Q_min [A-1]
129            @param max: Q_max [A-1]
130        """
131        ## Slit width in Q units
132        self.width  = width
133        ## Slit height in Q units
134        self.height = height
135        ## Q_min (Min Q-value for I(q))
136        self.min    = min
137        ## Q_max (Max Q_value for I(q))
138        self.max    = max
139        ## Number of Q bins
140        self.nbins  = nbins
141        ## Number of points used in the smearing computation
142        self.npts   = 1000
143        ## Smearing matrix
144        self._weights = None
145       
146    def _compute_matrix(self):
147        """
148            Compute the smearing matrix for the current I(q) array
149        """
150        weights = numpy.zeros([self.nbins, self.nbins])
151       
152        # Loop over all q-values
153        for i in range(self.nbins):
154            q = self.min + i*(self.max-self.min)/float(self.nbins-1.0)
155           
156            # For each q-value, compute the weight of each other q-bin
157            # in the I(q) array
158            npts_h = self.npts if self.height>0 else 1
159            npts_w = self.npts if self.width>0 else 1
160           
161            # If both height and width are great than zero,
162            # modify the number of points in each direction so
163            # that the total number of points is still what
164            # the user would expect (downgrade resolution)
165            if npts_h>1 and npts_w>1:
166                npts_h = int(math.ceil(math.sqrt(self.npts)))
167                npts_w = npts_h
168               
169            for k in range(npts_h):
170                if npts_h==1:
171                    shift_h = 0
172                else:
173                    shift_h = self.height/(float(npts_h-1.0)) * k
174                for j in range(npts_w):
175                    if npts_w==1:
176                        shift_w = 0
177                    else:
178                        shift_w = self.width/(float(npts_w-1.0)) * j
179                    q_shifted = math.sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) )
180                    q_i = int(math.floor( (q_shifted-self.min)/((self.max-self.min)/(self.nbins -1.0)) ))
181                   
182                    # Skip the entries outside our I(q) range
183                    #TODO: be careful with edge effect
184                    if q_i<self.nbins:
185                        weights[i][q_i] = weights[i][q_i]+1.0
186                               
187        self._weights = weights
188        return self._weights
189
190class SlitSmearer(_SlitSmearer):
191    """
192        Adaptor for slit smearing class and SANS data
193    """
194    def __init__(self, data1D):
195        """
196            Assumption: equally spaced bins of increasing q-values.
197           
198            @param data1D: data used to set the smearing parameters
199        """
200        # Initialization from parent class
201        super(SlitSmearer, self).__init__()
202       
203        ## Slit width
204        self.width = 0
205        if data1D.dxw is not None and len(data1D.dxw)==len(data1D.x):
206            self.width = data1D.dxw[0]
207            # Sanity check
208            for value in data1D.dxw:
209                if value != self.width:
210                    raise RuntimeError, "Slit smearing parameters must be the same for all data"
211               
212        ## Slit height
213        self.height = 0
214        if data1D.dxl is not None and len(data1D.dxl)==len(data1D.x):
215            self.height = data1D.dxl[0]
216            # Sanity check
217            for value in data1D.dxl:
218                if value != self.height:
219                    raise RuntimeError, "Slit smearing parameters must be the same for all data"
220       
221        ## Number of Q bins
222        self.nbins = len(data1D.x)
223        ## Minimum Q
224        self.min = data1D.x[0]
225        ## Maximum
226        self.max = data1D.x[len(data1D.x)-1]       
227
228
229class _QSmearer(_BaseSmearer):
230    """
231        Perform Gaussian Q smearing
232    """
233       
234    def __init__(self, nbins=None, width=None, min=None, max=None):
235        """
236            Initialization
237           
238            @param nbins: number of Q bins
239            @param width: array standard deviation in Q [A-1]
240            @param min: Q_min [A-1]
241            @param max: Q_max [A-1]
242        """
243        ## Standard deviation in Q [A-1]
244        self.width  = width
245        ## Q_min (Min Q-value for I(q))
246        self.min    = min
247        ## Q_max (Max Q_value for I(q))
248        self.max    = max
249        ## Number of Q bins
250        self.nbins  = nbins
251        ## Smearing matrix
252        self._weights = None
253       
254    def _compute_matrix(self):
255        """
256            Compute the smearing matrix for the current I(q) array
257        """
258        weights = numpy.zeros([self.nbins, self.nbins])
259       
260        # Loop over all q-values
261        step = (self.max-self.min)/float(self.nbins-1.0)
262        for i in range(self.nbins):
263            q = self.min + i*step
264            q_min = q - 0.5*step
265            q_max = q + 0.5*step
266           
267            for j in range(self.nbins):
268                q_j = self.min + j*step
269               
270                # Compute the fraction of the Gaussian contributing
271                # to the q bin between q_min and q_max
272                value =  scipy.special.erf( (q_max-q_j)/(math.sqrt(2.0)*self.width[j]) ) 
273                value -=scipy.special.erf( (q_min-q_j)/(math.sqrt(2.0)*self.width[j]) ) 
274
275                weights[i][j] += value
276                               
277        self._weights = weights
278        return self._weights
279       
280class QSmearer(_QSmearer):
281    """
282        Adaptor for Gaussian Q smearing class and SANS data
283    """
284    def __init__(self, data1D):
285        """
286            Assumption: equally spaced bins of increasing q-values.
287           
288            @param data1D: data used to set the smearing parameters
289        """
290        # Initialization from parent class
291        super(QSmearer, self).__init__()
292       
293        ## Resolution
294        self.width = numpy.zeros(len(data1D.x))
295        if data1D.dx is not None and len(data1D.dx)==len(data1D.x):
296            self.width = data1D.dx
297       
298        ## Number of Q bins
299        self.nbins = len(data1D.x)
300        ## Minimum Q
301        self.min = data1D.x[0]
302        ## Maximum
303        self.max = data1D.x[len(data1D.x)-1]       
304
305
306if __name__ == '__main__':
307    x = 0.001*numpy.arange(1,11)
308    y = 12.0-numpy.arange(1,11)
309    print x
310    #for i in range(10): print i, 0.001 + i*0.008/9.0
311    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
312
313   
314    s = _SlitSmearer(nbins=10, width=0.0, height=0.005, min=0.001, max=0.010)
315    #s = _QSmearer(nbins=10, width=0.001, min=0.001, max=0.010)
316    s._compute_matrix()
317
318    sy = s(y)
319    print sy
320   
321    if True:
322        for i in range(10):
323            print x[i], sy[i]
324            #print q, ' : ', s.weight(q), s._compute_iq(q)
325            #print q, ' : ', s(q), s._compute_iq(q)
326            #s._compute_iq(q)
327
328
329
330
Note: See TracBrowser for help on using the repository browser.