source: sasview/DataLoader/qsmearing.py @ 7d0c1a8

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 7d0c1a8 was d00f8ff, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

First version of slit smearing and resolution smearing

  • Property mode set to 100644
File size: 9.1 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"""
10import numpy
11import math
12import scipy.special
13
14def smear_selection(data1D):
15    """
16        Creates the right type of smearer according
17        to the data
18    """
19    pass
20
21class _BaseSmearer(object):
22   
23    def __init__(self):
24        self.nbins = 0
25        self._weights = None
26       
27    def _compute_matrix(self): return NotImplemented
28
29    def __call__(self, iq):
30        """
31            Return the smeared I(q) value at the given q.
32            The smeared I(q) is computed using a predetermined
33            smearing matrix for a particular binning.
34       
35            @param q: I(q) array
36            @return: smeared I(q)
37        """
38        # Sanity check
39        if len(iq) != self.nbins:
40            raise RuntimeError, "Invalid I(q) vector: inconsistent array length %s <> %s" % (len(iq), self.nbins)
41       
42        if self._weights == None:
43            self._compute_matrix()
44           
45        iq_smeared = numpy.zeros(self.nbins)
46        # Loop over q-values
47        for q_i in range(self.nbins):
48            sum = 0.0
49            counts = 0.0   
50           
51            for i in range(self.nbins):
52                sum += iq[i] * self._weights[q_i][i]
53                counts += self._weights[q_i][i]
54               
55            iq_smeared[q_i] = sum/counts
56               
57        return iq_smeared   
58   
59class _SlitSmearer(_BaseSmearer):
60    """
61        Slit smearing for I(q) array
62    """
63   
64    def __init__(self, nbins=None, width=None, height=None, min=None, max=None):
65        """
66            Initialization
67           
68            @param iq: I(q) array [cm-1]
69            @param width: slit width [A-1]
70            @param height: slit height [A-1]
71            @param min: Q_min [A-1]
72            @param max: Q_max [A-1]
73        """
74        ## Slit width in Q units
75        self.width  = width
76        ## Slit height in Q units
77        self.height = height
78        ## Q_min (Min Q-value for I(q))
79        self.min    = min
80        ## Q_max (Max Q_value for I(q))
81        self.max    = max
82        ## Number of Q bins
83        self.nbins  = nbins
84        ## Number of points used in the smearing computation
85        self.npts   = 1000
86        ## Smearing matrix
87        self._weights = None
88       
89    def _compute_matrix(self):
90        """
91            Compute the smearing matrix for the current I(q) array
92        """
93        weights = numpy.zeros([self.nbins, self.nbins])
94       
95        # Loop over all q-values
96        for i in range(self.nbins):
97            q = self.min + i*(self.max-self.min)/float(self.nbins-1.0)
98           
99            # For each q-value, compute the weight of each other q-bin
100            # in the I(q) array
101            npts_h = self.npts if self.height>0 else 1
102            npts_w = self.npts if self.width>0 else 1
103           
104            # If both height and width are great than zero,
105            # modify the number of points in each direction so
106            # that the total number of points is still what
107            # the user would expect (downgrade resolution)
108            if npts_h>1 and npts_w>1:
109                npts_h = int(math.ceil(math.sqrt(self.npts)))
110                npts_w = npts_h
111               
112            for k in range(npts_h):
113                if npts_h==1:
114                    shift_h = 0
115                else:
116                    shift_h = self.height/(float(npts_h-1.0)) * k
117                for j in range(npts_w):
118                    if npts_w==1:
119                        shift_w = 0
120                    else:
121                        shift_w = self.width/(float(npts_w-1.0)) * j
122                    q_shifted = math.sqrt( ((q-shift_w)*(q-shift_w) + shift_h*shift_h) )
123                    q_i = int(math.floor( (q_shifted-self.min)/((self.max-self.min)/(self.nbins -1.0)) ))
124                   
125                    # Skip the entries outside our I(q) range
126                    #TODO: be careful with edge effect
127                    if q_i<self.nbins:
128                        weights[i][q_i] = weights[i][q_i]+1.0
129                               
130        self._weights = weights
131        return self._weights
132
133class SlitSmearer(_SlitSmearer):
134    """
135        Adaptor for slit smearing class and SANS data
136    """
137    def __init__(self, data1D):
138        """
139            Assumption: equally spaced bins of increasing q-values.
140           
141            @param data1D: data used to set the smearing parameters
142        """
143        # Initialization from parent class
144        super(SlitSmearer, self).__init__()
145       
146        ## Slit width
147        self.width = 0
148        if data1D.dxw is not None and len(data1D.dxw)==len(data1D.x):
149            self.width = data1D.dxw[0]
150            # Sanity check
151            for value in data1D.dxw:
152                if value != self.width:
153                    raise RuntimeError, "Slit smearing parameters must be the same for all data"
154               
155        ## Slit height
156        self.height = 0
157        if data1D.dxl is not None and len(data1D.dxl)==len(data1D.x):
158            self.height = data1D.dxl[0]
159            # Sanity check
160            for value in data1D.dxl:
161                if value != self.height:
162                    raise RuntimeError, "Slit smearing parameters must be the same for all data"
163       
164        ## Number of Q bins
165        self.nbins = len(data1D.x)
166        ## Minimum Q
167        self.min = data1D.x[0]
168        ## Maximum
169        self.max = data1D.x[len(data1D.x)-1]       
170
171
172class _QSmearer(_BaseSmearer):
173    """
174        Perform Gaussian Q smearing
175    """
176       
177    def __init__(self, nbins=None, width=None, min=None, max=None):
178        """
179            Initialization
180           
181            @param nbins: number of Q bins
182            @param width: standard deviation in Q [A-1]
183            @param min: Q_min [A-1]
184            @param max: Q_max [A-1]
185        """
186        ## Standard deviation in Q [A-1]
187        self.width  = width
188        ## Q_min (Min Q-value for I(q))
189        self.min    = min
190        ## Q_max (Max Q_value for I(q))
191        self.max    = max
192        ## Number of Q bins
193        self.nbins  = nbins
194        ## Smearing matrix
195        self._weights = None
196       
197    def _compute_matrix(self):
198        """
199            Compute the smearing matrix for the current I(q) array
200        """
201        weights = numpy.zeros([self.nbins, self.nbins])
202       
203        # Loop over all q-values
204        step = (self.max-self.min)/float(self.nbins-1.0)
205        for i in range(self.nbins):
206            q = self.min + i*step
207            q_min = q - 0.5*step
208            q_max = q + 0.5*step
209           
210            for j in range(self.nbins):
211                q_j = self.min + j*step
212               
213                # Compute the fraction of the Gaussian contributing
214                # to the q bin between q_min and q_max
215                value =  scipy.special.erf( (q_max-q_j)/(math.sqrt(2.0)*self.width) ) 
216                value -=scipy.special.erf( (q_min-q_j)/(math.sqrt(2.0)*self.width) ) 
217
218                weights[i][j] += value
219                               
220        self._weights = weights
221        return self._weights
222       
223class QSmearer(_QSmearer):
224    """
225        Adaptor for Gaussian Q smearing class and SANS data
226    """
227    def __init__(self, data1D):
228        """
229            Assumption: equally spaced bins of increasing q-values.
230           
231            @param data1D: data used to set the smearing parameters
232        """
233        # Initialization from parent class
234        super(QSmearer, self).__init__()
235       
236        ## Slit width
237        self.width = 0
238        if data1D.dx is not None and len(data1D.dx)==len(data1D.x):
239            self.width = data1D.dx[0]
240            # Sanity check
241            for value in data1D.dx:
242                if value != self.width:
243                    raise RuntimeError, "dQ must be the same for all data"
244       
245        ## Number of Q bins
246        self.nbins = len(data1D.x)
247        ## Minimum Q
248        self.min = data1D.x[0]
249        ## Maximum
250        self.max = data1D.x[len(data1D.x)-1]       
251
252
253if __name__ == '__main__':
254    x = 0.001*numpy.arange(1,11)
255    y = 12.0-numpy.arange(1,11)
256    print x
257    #for i in range(10): print i, 0.001 + i*0.008/9.0
258    #for i in range(100): print i, int(math.floor( (i/ (100/9.0)) ))
259
260   
261    s = _SlitSmearer(nbins=10, width=0.0, height=0.005, min=0.001, max=0.010)
262    #s = _QSmearer(nbins=10, width=0.001, min=0.001, max=0.010)
263    s._compute_matrix()
264
265    sy = s(y)
266    print sy
267   
268    if True:
269        for i in range(10):
270            print x[i], sy[i]
271            #print q, ' : ', s.weight(q), s._compute_iq(q)
272            #print q, ' : ', s(q), s._compute_iq(q)
273            #s._compute_iq(q)
274
275
276
277
Note: See TracBrowser for help on using the repository browser.