source: sasview/DataLoader/manipulations.py @ 1fbb343

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 1fbb343 was ee3b2c0c, checked in by Jae Cho <jhjcho@…>, 16 years ago

reversed the last change causing another problem

  • Property mode set to 100644
File size: 43.2 KB
Line 
1"""
2    Data manipulations for 2D data sets.
3    Using the meta data information, various types of averaging
4    are performed in Q-space
5"""
6
7"""
8This software was developed by the University of Tennessee as part of the
9Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
10project funded by the US National Science Foundation.
11
12See the license text in license.txt
13
14copyright 2008, University of Tennessee
15"""
16#TODO: copy the meta data from the 2D object to the resulting 1D object
17
18from data_info import plottable_2D, Data1D
19import math
20import numpy
21
22def get_q(dx, dy, det_dist, wavelength):
23    """
24        @param dx: x-distance from beam center [mm]
25        @param dy: y-distance from beam center [mm]
26        @return: q-value at the given position
27    """
28    # Distance from beam center in the plane of detector
29    plane_dist = math.sqrt(dx*dx + dy*dy)
30    # Half of the scattering angle
31    theta      = 0.5*math.atan(plane_dist/det_dist)
32    return (4.0*math.pi/wavelength)*math.sin(theta)
33
34def get_q_compo(dx, dy, det_dist, wavelength,compo=None):
35    #This reduces tiny error at very large q.
36    #Implementation of this func is not started yet.<--ToDo
37    if dy==0:
38        if dx>=0:
39            angle_xy=0
40        else:
41            angle_xy=math.pi
42    else:
43        angle_xy=math.atan(dx/dy)
44       
45    if compo=="x":
46        out=get_q(dx, dy, det_dist, wavelength)*cos(angle_xy)
47    elif compo=="y":
48        out=get_q(dx, dy, det_dist, wavelength)*sin(angle_xy)
49    else:
50        out=get_q(dx, dy, det_dist, wavelength)
51    return out
52   
53class _Slab(object):
54    """
55        Compute average I(Q) for a region of interest
56    """
57    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0, bin_width=0.001):
58        # Minimum Qx value [A-1]
59        self.x_min = x_min
60        # Maximum Qx value [A-1]
61        self.x_max = x_max
62        # Minimum Qy value [A-1]
63        self.y_min = y_min
64        # Maximum Qy value [A-1]
65        self.y_max = y_max
66        # Bin width (step size) [A-1]
67        self.bin_width = bin_width
68        # If True, I(|Q|) will be return, otherwise, negative q-values are allowed
69        self.fold = False
70       
71    def __call__(self, data2D): return NotImplemented
72       
73    def _avg(self, data2D, maj):
74        """
75             Compute average I(Q_maj) for a region of interest.
76             The major axis is defined as the axis of Q_maj.
77             The minor axis is the axis that we average over.
78             
79             @param data2D: Data2D object
80             @param maj_min: min value on the major axis
81             @return: Data1D object
82        """
83        if len(data2D.detector) != 1:
84            raise RuntimeError, "_Slab._avg: invalid number of detectors: %g" % len(data2D.detector)
85       
86        pixel_width_x = data2D.detector[0].pixel_size.x
87        pixel_width_y = data2D.detector[0].pixel_size.y
88        det_dist    = data2D.detector[0].distance
89        wavelength  = data2D.source.wavelength
90        center_x    = data2D.detector[0].beam_center.x/pixel_width_x
91        center_y    = data2D.detector[0].beam_center.y/pixel_width_y
92               
93        # Build array of Q intervals
94        if maj=='x':
95            nbins = int(math.ceil((self.x_max-self.x_min)/self.bin_width))
96            qbins = self.bin_width*numpy.arange(nbins)+self.x_min
97        elif maj=='y':
98            nbins = int(math.ceil((self.y_max-self.y_min)/self.bin_width))
99            qbins = self.bin_width*numpy.arange(nbins)+self.y_min           
100        else:
101            raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj)
102                               
103        x  = numpy.zeros(nbins)
104        y  = numpy.zeros(nbins)
105        err_y = numpy.zeros(nbins)
106        y_counts = numpy.zeros(nbins)
107                                               
108        for i in range(numpy.size(data2D.data,1)):
109            # Min and max x-value for the pixel
110            minx = pixel_width_x*(i - center_x)
111            maxx = pixel_width_x*(i+1.0 - center_x)
112           
113            qxmin = get_q(minx, 0.0, det_dist, wavelength)
114            qxmax = get_q(maxx, 0.0, det_dist, wavelength)
115           
116            # Get the count fraction in x for that pixel
117            frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
118            frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
119            frac_x = frac_max - frac_min
120           
121            if frac_x == 0: 
122                continue
123           
124            if maj=='x':
125                dx = pixel_width_x*(i+0.5 - center_x)
126                q_value = get_q(dx, 0.0, det_dist, wavelength)
127                if self.fold==False and dx<0:
128                    q_value = -q_value
129                i_q = int(math.ceil((q_value-self.x_min)/self.bin_width)) - 1
130               
131                if i_q<0 or i_q>=nbins:
132                    continue
133                       
134            for j in range(numpy.size(data2D.data,0)):
135                # Min and max y-value for the pixel
136                miny = pixel_width_y*(j - center_y)
137                maxy = pixel_width_y*(j+1.0 - center_y)
138
139                qymin = get_q(0.0, miny, det_dist, wavelength)
140                qymax = get_q(0.0, maxy, det_dist, wavelength)
141               
142                # Get the count fraction in x for that pixel
143                frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
144                frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
145                frac_y = frac_max - frac_min
146               
147                frac = frac_x * frac_y
148               
149                if frac == 0:
150                    continue
151
152                if maj=='y':
153                    dy = pixel_width_y*(j+0.5 - center_y)
154                    q_value = get_q(0.0, dy, det_dist, wavelength)
155                    if self.fold==False and dy<0:
156                        q_value = -q_value
157                    i_q = int(math.ceil((q_value-self.y_min)/self.bin_width)) - 1
158                   
159                    if i_q<0 or i_q>=nbins:
160                        continue
161           
162                x[i_q]          = q_value
163                y[i_q]         += frac * data2D.data[j][i]
164                if data2D.err_data == None or data2D.err_data[j][i]==0.0:
165                    err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
166                else:
167                    err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
168                y_counts[i_q]  += frac
169
170        # Average the sums
171        for i in range(nbins):
172            if y_counts[i]>0:
173                err_y[i] = math.sqrt(err_y[i])/y_counts[i]
174                y[i]     = y[i]/y_counts[i]
175       
176        return Data1D(x=x, y=y, dy=err_y)
177       
178class SlabY(_Slab):
179    """
180        Compute average I(Qy) for a region of interest
181    """
182    def __call__(self, data2D):
183        """
184             Compute average I(Qy) for a region of interest
185             
186             @param data2D: Data2D object
187             @return: Data1D object
188        """
189        return self._avg(data2D, 'y')
190       
191class SlabX(_Slab):
192    """
193        Compute average I(Qx) for a region of interest
194    """
195    def __call__(self, data2D):
196        """
197             Compute average I(Qx) for a region of interest
198             
199             @param data2D: Data2D object
200             @return: Data1D object
201        """
202        return self._avg(data2D, 'x') 
203       
204class Boxsum(object):
205    """
206        Perform the sum of counts in a 2D region of interest.
207    """
208    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
209        # Minimum Qx value [A-1]
210        self.x_min = x_min
211        # Maximum Qx value [A-1]
212        self.x_max = x_max
213        # Minimum Qy value [A-1]
214        self.y_min = y_min
215        # Maximum Qy value [A-1]
216        self.y_max = y_max
217
218    def __call__(self, data2D):
219        """
220             Perform the sum in the region of interest
221             
222             @param data2D: Data2D object
223             @return: number of counts, error on number of counts
224        """
225        y, err_y, y_counts = self._sum(data2D)
226       
227        # Average the sums
228        counts = 0 if y_counts==0 else y
229        error  = 0 if y_counts==0 else math.sqrt(err_y)
230       
231        return counts, error
232       
233    def _sum(self, data2D):
234        """
235             Perform the sum in the region of interest
236             @param data2D: Data2D object
237             @return: number of counts, error on number of counts, number of entries summed
238        """
239        if len(data2D.detector) != 1:
240            raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
241       
242        pixel_width_x = data2D.detector[0].pixel_size.x
243        pixel_width_y = data2D.detector[0].pixel_size.y
244        det_dist    = data2D.detector[0].distance
245        wavelength  = data2D.source.wavelength
246        center_x    = data2D.detector[0].beam_center.x/pixel_width_x
247        center_y    = data2D.detector[0].beam_center.y/pixel_width_y
248               
249        y  = 0.0
250        err_y = 0.0
251        y_counts = 0.0
252        sign=1
253        for i in range(numpy.size(data2D.data,1)):
254            # Min and max x-value for the pixel
255            minx = pixel_width_x*(i - center_x)
256            maxx = pixel_width_x*(i+1.0 - center_x)
257            if minx>=0:
258                sign=1
259            else:
260                sign=-1           
261            qxmin = sign*get_q(minx, 0.0, det_dist, wavelength)
262            if maxx>=0:
263                sign=1
264            else:
265                sign=-1
266            qxmax = sign*get_q(maxx, 0.0, det_dist, wavelength)
267           
268            # Get the count fraction in x for that pixel
269            frac_min = get_pixel_fraction_square(self.x_min, qxmin, qxmax)
270            frac_max = get_pixel_fraction_square(self.x_max, qxmin, qxmax)
271            frac_x = frac_max - frac_min
272           
273            for j in range(numpy.size(data2D.data,0)):
274                # Min and max y-value for the pixel
275                miny = pixel_width_y*(j - center_y)
276                maxy = pixel_width_y*(j+1.0 - center_y)
277                if miny>=0:
278                    sign=1
279                else:
280                    sign=-1           
281
282                qymin = sign*get_q(0.0, miny, det_dist, wavelength)
283                if maxy>=0:
284                    sign=1
285                else:
286                    sign=-1           
287               
288                qymax = sign*get_q(0.0, maxy, det_dist, wavelength)
289               
290                # Get the count fraction in x for that pixel
291                frac_min = get_pixel_fraction_square(self.y_min, qymin, qymax)
292                frac_max = get_pixel_fraction_square(self.y_max, qymin, qymax)
293                frac_y = frac_max - frac_min
294               
295                frac = frac_x * frac_y
296
297                y += frac * data2D.data[j][i]
298                if data2D.err_data == None or data2D.err_data[j][i]==0.0:
299                    err_y += frac * frac * math.fabs(data2D.data[j][i])
300                else:
301                    err_y += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
302                y_counts += frac
303       
304        return y, err_y, y_counts
305     
306class Boxavg(Boxsum):
307    """
308        Perform the average of counts in a 2D region of interest.
309    """
310    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
311        super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
312
313    def __call__(self, data2D):
314        """
315             Perform the sum in the region of interest
316             
317             @param data2D: Data2D object
318             @return: average counts, error on average counts
319        """
320        y, err_y, y_counts = self._sum(data2D)
321       
322        # Average the sums
323        counts = 0 if y_counts==0 else y/y_counts
324        error  = 0 if y_counts==0 else math.sqrt(err_y)/y_counts
325       
326        return counts, error
327       
328def get_pixel_fraction_square(x, xmin, xmax):
329    """
330         Return the fraction of the length
331         from xmin to x.
332         
333             A            B
334         +-----------+---------+
335         xmin        x         xmax
336         
337         @param x: x-value
338         @param xmin: minimum x for the length considered
339         @param xmax: minimum x for the length considered
340         @return: (x-xmin)/(xmax-xmin) when xmin < x < xmax
341         
342    """
343    if x<=xmin:
344        return 0.0
345    if x>xmin and x<xmax:
346        return (x-xmin)/(xmax-xmin)
347    else:
348        return 1.0
349
350
351class CircularAverage(object):
352    """
353        Perform circular averaging on 2D data
354       
355        The data returned is the distribution of counts
356        as a function of Q
357    """
358    def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.001):
359        # Minimum radius included in the average [A-1]
360        self.r_min = r_min
361        # Maximum radius included in the average [A-1]
362        self.r_max = r_max
363        # Bin width (step size) [A-1]
364        self.bin_width = bin_width
365
366    def __call__(self, data2D):
367        """
368            Perform circular averaging on the data
369           
370            @param data2D: Data2D object
371            @return: Data1D object
372        """
373        if len(data2D.detector) != 1:
374            raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
375       
376        pixel_width_x = data2D.detector[0].pixel_size.x
377        pixel_width_y = data2D.detector[0].pixel_size.y
378        det_dist    = data2D.detector[0].distance
379        wavelength  = data2D.source.wavelength
380        center_x    = data2D.detector[0].beam_center.x/pixel_width_x
381        center_y    = data2D.detector[0].beam_center.y/pixel_width_y
382       
383        # Find out the maximum Q range
384        xwidth = numpy.size(data2D.data,1)*pixel_width_x
385        dx_max = xwidth - data2D.detector[0].beam_center.x
386        if xwidth-dx_max>dx_max:
387            dx_max = xwidth-dx_max
388           
389        ywidth = numpy.size(data2D.data,0)*pixel_width_y
390        dy_max = ywidth - data2D.detector[0].beam_center.y
391        if ywidth-dy_max>dy_max:
392            dy_max = ywidth-dy_max
393       
394        qmax = get_q(dx_max, dy_max, det_dist, wavelength)
395       
396        # Build array of Q intervals
397        nbins = int(math.ceil((qmax-self.r_min)/self.bin_width))
398        qbins = self.bin_width*numpy.arange(nbins)+self.r_min
399       
400        x  = numpy.zeros(nbins)
401        y  = numpy.zeros(nbins)
402        err_y = numpy.zeros(nbins)
403        y_counts = numpy.zeros(nbins)
404       
405        for i in range(numpy.size(data2D.data,1)):
406            dx = pixel_width_x*(i+0.5 - center_x)
407           
408            # Min and max x-value for the pixel
409            minx = pixel_width_x*(i - center_x)
410            maxx = pixel_width_x*(i+1.0 - center_x)
411           
412            for j in range(numpy.size(data2D.data,0)):
413                dy = pixel_width_y*(j+0.5 - center_y)
414           
415                q_value = get_q(dx, dy, det_dist, wavelength)
416           
417                # Min and max y-value for the pixel
418                miny = pixel_width_y*(j - center_y)
419                maxy = pixel_width_y*(j+1.0 - center_y)
420               
421                # Calculate the q-value for each corner
422                # q_[x min or max][y min or max]
423                q_00 = get_q(minx, miny, det_dist, wavelength)
424                q_01 = get_q(minx, maxy, det_dist, wavelength)
425                q_10 = get_q(maxx, miny, det_dist, wavelength)
426                q_11 = get_q(maxx, maxy, det_dist, wavelength)
427               
428                # Look for intercept between each side of the pixel
429                # and the constant-q ring for qmax
430                frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11)
431               
432                # Look for intercept between each side of the pixel
433                # and the constant-q ring for qmin
434                frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11)
435               
436                # We are interested in the region between qmin and qmax
437                # therefore the fraction of the surface of the pixel
438                # that we will use to calculate the number of counts to
439                # include is given by:
440                frac = frac_max - frac_min
441
442                i_q = int(math.ceil((q_value-self.r_min)/self.bin_width)) - 1
443           
444                x[i_q]          = q_value
445                y[i_q]         += frac * data2D.data[j][i]
446                if data2D.err_data == None or data2D.err_data[j][i]==0.0:
447                    err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i])
448                else:
449                    err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
450                y_counts[i_q]  += frac
451       
452        # Average the sums
453        for i in range(nbins):
454            if y_counts[i]>0:
455                err_y[i] = math.sqrt(err_y[i])/y_counts[i]
456                y[i]     = y[i]/y_counts[i]
457       
458        return Data1D(x=x, y=y, dy=err_y)
459   
460
461class Ring(object):
462    """
463        Defines a ring on a 2D data set.
464        The ring is defined by r_min, r_max, and
465        the position of the center of the ring.
466       
467        The data returned is the distribution of counts
468        around the ring as a function of phi.
469       
470    """
471   
472    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0,nbins=20 ):
473        # Minimum radius
474        self.r_min = r_min
475        # Maximum radius
476        self.r_max = r_max
477        # Center of the ring in x
478        self.center_x = center_x
479        # Center of the ring in y
480        self.center_y = center_y
481        # Number of angular bins
482        self.nbins_phi = nbins_phi
483       
484    def __call__(self, data2D):
485        """
486            Apply the ring to the data set.
487            Returns the angular distribution for a given q range
488           
489            @param data2D: Data2D object
490            @return: Data1D object
491        """
492        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
493            raise RuntimeError, "Ring averaging only take plottable_2D objects"
494       
495        data = data2D.data
496        qmin = self.r_min
497        qmax = self.r_max
498       
499        if len(data2D.detector) != 1:
500            raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
501        pixel_width_x = data2D.detector[0].pixel_size.x
502        pixel_width_y = data2D.detector[0].pixel_size.y
503        det_dist = data2D.detector[0].distance
504        wavelength = data2D.source.wavelength
505        #center_x = self.center_x/pixel_width_x
506        #center_y = self.center_y/pixel_width_y
507        center_x    = data2D.detector[0].beam_center.x/pixel_width_x
508        center_y    = data2D.detector[0].beam_center.y/pixel_width_y
509       
510       
511        phi_bins   = numpy.zeros(self.nbins_phi)
512        phi_counts = numpy.zeros(self.nbins_phi)
513        phi_values = numpy.zeros(self.nbins_phi)
514        phi_err    = numpy.zeros(self.nbins_phi)
515       
516        for i in range(numpy.size(data,1)):
517            dx = pixel_width_x*(i+0.5 - center_x)
518           
519            # Min and max x-value for the pixel
520            minx = pixel_width_x*(i - center_x)
521            maxx = pixel_width_x*(i+1.0 - center_x)
522           
523            for j in range(numpy.size(data,0)):
524                dy = pixel_width_y*(j+0.5 - center_y)
525           
526                q_value = get_q(dx, dy, det_dist, wavelength)
527           
528                # Min and max y-value for the pixel
529                miny = pixel_width_y*(j - center_y)
530                maxy = pixel_width_y*(j+1.0 - center_y)
531               
532                # Calculate the q-value for each corner
533                # q_[x min or max][y min or max]
534                q_00 = get_q(minx, miny, det_dist, wavelength)
535                q_01 = get_q(minx, maxy, det_dist, wavelength)
536                q_10 = get_q(maxx, miny, det_dist, wavelength)
537                q_11 = get_q(maxx, maxy, det_dist, wavelength)
538               
539                # Look for intercept between each side of the pixel
540                # and the constant-q ring for qmax
541                frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
542               
543                # Look for intercept between each side of the pixel
544                # and the constant-q ring for qmin
545                frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
546               
547                # We are interested in the region between qmin and qmax
548                # therefore the fraction of the surface of the pixel
549                # that we will use to calculate the number of counts to
550                # include is given by:
551               
552                frac = frac_max - frac_min
553
554                i_phi = int(math.ceil(self.nbins_phi*(math.atan2(dy, dx)+math.pi)/(2.0*math.pi))) - 1
555           
556                phi_bins[i_phi] += frac * data[j][i]
557               
558                if data2D.err_data == None or data2D.err_data[j][i]==0.0:
559                    phi_err[i_phi] += frac * frac * math.fabs(data2D.data[j][i])
560                else:
561                    phi_err[i_phi] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
562                phi_counts[i_phi] += frac
563       
564        for i in range(self.nbins_phi):
565            phi_bins[i] = phi_bins[i] / phi_counts[i]
566            phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
567            phi_values[i] = 2.0*math.pi/self.nbins_phi*(1.0*i + 0.5)-math.pi # move the pi back to -pi <-->+pi
568           
569        return Data1D(x=phi_values, y=phi_bins, dy=phi_err)
570   
571def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
572    """
573        Returns the fraction of the pixel defined by
574        the four corners (q_00, q_01, q_10, q_11) that
575        has q < qmax.
576       
577                q_01                q_11
578        y=1         +--------------+
579                    |              |
580                    |              |
581                    |              |
582        y=0         +--------------+
583                q_00                q_01
584       
585                    x=0            x=1
586       
587    """
588   
589    # y side for x = minx
590    x_0 = get_intercept(qmax, q_00, q_01)
591    # y side for x = maxx
592    x_1 = get_intercept(qmax, q_10, q_11)
593   
594    # x side for y = miny
595    y_0 = get_intercept(qmax, q_00, q_10)
596    # x side for y = maxy
597    y_1 = get_intercept(qmax, q_01, q_11)
598   
599    # surface fraction for a 1x1 pixel
600    frac_max = 0
601   
602    if x_0 and x_1:
603        frac_max = (x_0+x_1)/2.0
604   
605    elif y_0 and y_1:
606        frac_max = (y_0+y_1)/2.0
607   
608    elif x_0 and y_0:
609        if q_00 < q_10:
610            frac_max = x_0*y_0/2.0
611        else:
612            frac_max = 1.0-x_0*y_0/2.0
613   
614    elif x_0 and y_1:
615        if q_00 < q_10:
616            frac_max = x_0*y_1/2.0
617        else:
618            frac_max = 1.0-x_0*y_1/2.0
619   
620    elif x_1 and y_0:
621        if q_00 > q_10:
622            frac_max = x_1*y_0/2.0
623        else:
624            frac_max = 1.0-x_1*y_0/2.0
625   
626    elif x_1 and y_1:
627        if q_00 < q_10:
628            frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
629        else:
630            frac_max = (1.0-x_1)*(1.0-y_1)/2.0
631           
632    # If we make it here, there is no intercept between
633    # this pixel and the constant-q ring. We only need
634    # to know if we have to include it or exclude it.
635    elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
636        frac_max = 1.0
637   
638    return frac_max
639             
640def get_intercept(q, q_0, q_1):
641    """
642        Returns the fraction of the side at which the
643        q-value intercept the pixel, None otherwise.
644        The values returned is the fraction ON THE SIDE
645        OF THE LOWEST Q.
646       
647       
648       
649                A        B   
650         +-----------+--------+
651         0                    1     <--- pixel size
652         
653        Q_0 -------- Q ----- Q_1    <--- equivalent Q range
654       
655       
656        if Q_1 > Q_0, A is returned
657        if Q_1 < Q_0, B is returned
658       
659        if Q is outside the range of [Q_0, Q_1], None is returned
660         
661    """
662    if q_1 > q_0:
663        if (q > q_0 and q <= q_1):
664            return (q-q_0)/(q_1 - q_0)   
665    else:
666        if (q > q_1 and q <= q_0):
667            return (q-q_1)/(q_0 - q_1)
668    return None
669   
670#This class can be removed.
671class _Sectorold:
672    """
673        Defines a sector region on a 2D data set.
674        The sector is defined by r_min, r_max, phi_min, phi_max,
675        and the position of the center of the ring.         
676        Phi is defined between 0 and 2pi
677    """
678    def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
679        self.r_min = r_min
680        self.r_max = r_max
681        self.phi_min = phi_min
682        self.phi_max = phi_max
683        self.nbins = nbins
684       
685    def _agv(self, data2D, run='phi'):
686        """
687            Perform sector averaging.
688           
689            @param data2D: Data2D object
690            @param run:  define the varying parameter ('phi' or 'q')
691            @return: Data1D object
692        """
693        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
694            raise RuntimeError, "Ring averaging only take plottable_2D objects"
695                   
696        data = data2D.data     
697        qmax = self.r_max
698        qmin = self.r_min
699       
700        if len(data2D.detector) != 1:
701            raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
702        pixel_width_x = data2D.detector[0].pixel_size.x
703        pixel_width_y = data2D.detector[0].pixel_size.y
704        det_dist      = data2D.detector[0].distance
705        wavelength    = data2D.source.wavelength
706        center_x      = data2D.detector[0].beam_center.x/pixel_width_x
707        center_y      = data2D.detector[0].beam_center.y/pixel_width_y
708       
709        y        = numpy.zeros(self.nbins)
710        y_counts = numpy.zeros(self.nbins)
711        x        = numpy.zeros(self.nbins)
712        y_err    = numpy.zeros(self.nbins)
713       
714        for i in range(numpy.size(data,1)):
715            dx = pixel_width_x*(i+0.5 - center_x)
716           
717            # Min and max x-value for the pixel
718            minx = pixel_width_x*(i - center_x)
719            maxx = pixel_width_x*(i+1.0 - center_x)
720           
721            for j in range(numpy.size(data,0)):
722                dy = pixel_width_y*(j+0.5 - center_y)
723           
724                q_value = get_q(dx, dy, det_dist, wavelength)
725
726                # Min and max y-value for the pixel
727                miny = pixel_width_y*(j - center_y)
728                maxy = pixel_width_y*(j+1.0 - center_y)
729               
730                # Calculate the q-value for each corner
731                # q_[x min or max][y min or max]
732                q_00 = get_q(minx, miny, det_dist, wavelength)
733                q_01 = get_q(minx, maxy, det_dist, wavelength)
734                q_10 = get_q(maxx, miny, det_dist, wavelength)
735                q_11 = get_q(maxx, maxy, det_dist, wavelength)
736               
737                # Look for intercept between each side of the pixel
738                # and the constant-q ring for qmax
739                frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
740               
741                # Look for intercept between each side of the pixel
742                # and the constant-q ring for qmin
743                frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
744               
745                # We are interested in the region between qmin and qmax
746                # therefore the fraction of the surface of the pixel
747                # that we will use to calculate the number of counts to
748                # include is given by:
749               
750                frac = frac_max - frac_min
751
752                # Compute phi and check whether it's within the limits
753                phi_value=math.atan2(dy,dx)+math.pi
754 #               if phi_value<self.phi_min or phi_value>self.phi_max:               
755                if phi_value<self.phi_min or phi_value>self.phi_max:
756                    continue
757                                                   
758                # Check which type of averaging we need
759                if run.lower()=='phi': 
760                    i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
761                else:
762                    # If we don't need this pixel, skip the rest of the work
763                    #TODO: an improvement here would be to compute the average
764                    # Q for the pixel from the part that is covered by
765                    # the ring defined by q_min/q_max rather than the complete
766                    # pixel
767                    if q_value<self.r_min or q_value>self.r_max:
768                        continue
769                    i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1
770           
771                try:
772                    y[i_bin] += frac * data[j][i]
773                except:
774                    import sys
775                    print sys.exc_value
776                    print i_bin, frac
777               
778                if data2D.err_data == None or data2D.err_data[j][i]==0.0:
779                    y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
780                else:
781                    y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
782                y_counts[i_bin] += frac
783       
784        for i in range(self.nbins):
785            y[i] = y[i] / y_counts[i]
786            y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
787            # Check which type of averaging we need
788            if run.lower()=='phi':
789                x[i] = (self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min
790            else:
791                x[i] = (self.r_max-self.r_min)/self.nbins*(1.0*i + 0.5)+self.r_min
792           
793        return Data1D(x=x, y=y, dy=y_err)
794       
795class _Sector:
796    """
797        Defines a sector region on a 2D data set.
798        The sector is defined by r_min, r_max, phi_min, phi_max,
799        and the position of the center of the ring
800        where phi_min and phi_max are defined by the right and left lines wrt central line
801        and phi_max could be less than phi_min.
802       
803        Phi is defined between 0 and 2pi
804    """
805    def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20):
806        self.r_min = r_min
807        self.r_max = r_max
808        self.phi_min = phi_min
809        self.phi_max = phi_max
810        self.nbins = nbins
811       
812    def _agv(self, data2D, run='phi'):
813        """
814            Perform sector averaging.
815           
816            @param data2D: Data2D object
817            @param run:  define the varying parameter ('phi' or 'q')
818            @return: Data1D object
819        """
820        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
821            raise RuntimeError, "Ring averaging only take plottable_2D objects"
822                   
823        data = data2D.data     
824        qmax = self.r_max
825        qmin = self.r_min
826       
827        if len(data2D.detector) != 1:
828            raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector)
829        pixel_width_x = data2D.detector[0].pixel_size.x
830        pixel_width_y = data2D.detector[0].pixel_size.y
831        det_dist      = data2D.detector[0].distance
832        wavelength    = data2D.source.wavelength
833        center_x      = data2D.detector[0].beam_center.x/pixel_width_x
834        center_y      = data2D.detector[0].beam_center.y/pixel_width_y
835       
836        y        = numpy.zeros(self.nbins)
837        y_counts = numpy.zeros(self.nbins)
838        x        = numpy.zeros(self.nbins)
839        y_err    = numpy.zeros(self.nbins)
840       
841        # This If finds qmax within ROI defined by sector lines
842        temp=0 #to find qmax within ROI or phimax and phimin
843        temp0=1000000
844        temp1=0
845        for i in range(numpy.size(data,1)): 
846            dx = pixel_width_x*(i+0.5 - center_x)                 
847            for j in range(numpy.size(data,0)):
848                   
849                dy = pixel_width_y*(j+0.5 - center_y)
850                q_value = get_q(dx, dy, det_dist, wavelength)
851                # Compute phi and check whether it's within the limits
852                phi_value=math.atan2(dy,dx)+math.pi
853                if self.phi_max>2*math.pi:
854                    self.phi_max=self.phi_max-2*math.pi
855                if self.phi_min<0:
856                    self.phi_max=self.phi_max+2*math.pi
857               
858                #In case of two ROI (symmetric major and minor regions)(for 'q2')
859                if run.lower()=='q2':
860                    if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
861                        temp_max=self.phi_max+math.pi
862                        temp_min=self.phi_min+math.pi
863                    else:
864                        temp_max=self.phi_max
865                        temp_min=self.phi_min
866                       
867                    if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)):
868                        if (phi_value<temp_min  or phi_value>temp_max):
869                             if (phi_value<temp_min-math.pi  or phi_value>temp_max-math.pi):
870                                 continue
871                    if (self.phi_max<self.phi_min):
872                        tmp_max=self.phi_max+math.pi
873                        tmp_min=self.phi_min-math.pi
874                    else:
875                        tmp_max=self.phi_max
876                        tmp_min=self.phi_min
877                    if (tmp_min<math.pi and tmp_max>math.pi):
878                        if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
879                            continue
880                #In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'.
881                elif run.lower()=='q': 
882                    if (self.phi_max>=self.phi_min):
883                        if (phi_value<self.phi_min  or phi_value>self.phi_max):
884                            continue
885                    else:
886                        if (phi_value<self.phi_min and phi_value>self.phi_max):
887                            continue   
888                if q_value<qmin or q_value>qmax:
889                        continue                   
890                       
891                if run.lower()=='phi':
892                    if temp1<phi_value:
893                        temp1=phi_value
894                    if temp0>phi_value:
895                        temp0=phi_value   
896                                                                                         
897                elif temp<q_value:
898                    temp=q_value
899                   
900        if run.lower()=='phi':
901            self.phi_max=temp1
902            self.phi_min=temp0
903        else:
904            qmax=temp
905        #Beam center is already corrected, but the calculation below assumed it was not.
906        # Thus Beam center shifted back to uncorrected value. ToDo: cleanup this codes.
907        center_x=center_x+0.5 
908        center_y=center_y+0.5         
909        for i in range(numpy.size(data,1)):
910            dx = pixel_width_x*(i+0.5 - center_x)
911           
912            # Min and max x-value for the pixel
913            minx = pixel_width_x*(i - center_x)
914            maxx = pixel_width_x*(i+1.0 - center_x)
915           
916            for j in range(numpy.size(data,0)):
917                dy = pixel_width_y*(j+0.5 - center_y)
918           
919                q_value = get_q(dx, dy, det_dist, wavelength)
920
921                # Min and max y-value for the pixel
922                miny = pixel_width_y*(j - center_y)
923                maxy = pixel_width_y*(j+1.0 - center_y)
924               
925                # Calculate the q-value for each corner
926                # q_[x min or max][y min or max]
927                q_00 = get_q(minx, miny, det_dist, wavelength)
928                q_01 = get_q(minx, maxy, det_dist, wavelength)
929                q_10 = get_q(maxx, miny, det_dist, wavelength)
930                q_11 = get_q(maxx, maxy, det_dist, wavelength)
931               
932                # Compute phi and check whether it's within the limits
933                phi_value=math.atan2(dy,dx)+math.pi
934                if self.phi_max>2*math.pi:
935                    self.phi_max=self.phi_max-2*math.pi
936                if self.phi_min<0:
937                    self.phi_max=self.phi_max+2*math.pi
938                   
939                # Look for intercept between each side of the pixel
940                # and the constant-q ring for qmax
941                frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11)
942               
943                # Look for intercept between each side of the pixel
944                # and the constant-q ring for qmin
945                frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11)
946               
947                # We are interested in the region between qmin and qmax
948                # therefore the fraction of the surface of the pixel
949                # that we will use to calculate the number of counts to
950                # include is given by:
951               
952                frac = frac_max - frac_min
953
954                #In case of two ROI (symmetric major and minor regions)(for 'q2')
955                if run.lower()=='q2':
956                    if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)):
957                        temp_max=self.phi_max+math.pi
958                        temp_min=self.phi_min+math.pi
959                    else:
960                        temp_max=self.phi_max
961                        temp_min=self.phi_min
962                       
963                    if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)):
964                        if (phi_value<temp_min  or phi_value>temp_max):
965                            if (phi_value<temp_min-math.pi  or phi_value>temp_max-math.pi):
966                                continue
967                    if (self.phi_max<self.phi_min):
968                        tmp_max=self.phi_max+math.pi
969                        tmp_min=self.phi_min-math.pi
970                    else:
971                        tmp_max=self.phi_max
972                        tmp_min=self.phi_min
973                    if (tmp_min<math.pi and tmp_max>math.pi):
974                        if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)):
975                            continue
976                #In case of one ROI (major only)(i.e.,for 'q' and 'phi')
977                else: 
978                    if (self.phi_max>=self.phi_min):
979                        if (phi_value<self.phi_min  or phi_value>self.phi_max):
980                            continue
981                           
982                    else:
983                        if (phi_value<self.phi_min and phi_value>self.phi_max):
984                            continue
985                #print "qmax=",qmax,qmin       
986
987                if q_value<qmin or q_value>qmax:
988                    continue
989                                                   
990                # Check which type of averaging we need
991                if run.lower()=='phi': 
992                    i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1
993                else:
994                    # If we don't need this pixel, skip the rest of the work
995                    #TODO: an improvement here would be to compute the average
996                    # Q for the pixel from the part that is covered by
997                    # the ring defined by q_min/q_max rather than the complete
998                    # pixel
999                    i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1
1000                           
1001                try:
1002                    y[i_bin] += frac * data[j][i]
1003                except:
1004                    import sys
1005                    print sys.exc_value
1006                    print i_bin, frac
1007               
1008                if data2D.err_data == None or data2D.err_data[j][i]==0.0:
1009                    y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i])
1010                else:
1011                    y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i]
1012                y_counts[i_bin] += frac
1013       
1014        for i in range(self.nbins):
1015            y[i] = y[i] / y_counts[i]
1016            y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
1017            # Check which type of averaging we need
1018            if run.lower()=='phi':
1019                #Calculate x[i] and put back the origin of angle back to the right hand side (from the left).
1020                x[i] = ((self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min-2*math.pi)*180/math.pi
1021                if x[i]<0:
1022                    x[i]=180+x[i]
1023            else:
1024                x[i] = (qmax-qmin)/self.nbins*(1.0*i + 0.5)+qmin
1025           
1026        return Data1D(x=x, y=y, dy=y_err)
1027               
1028class SectorPhi(_Sector):
1029    """
1030        Sector average as a function of phi.
1031        I(phi) is return and the data is averaged over Q.
1032       
1033        A sector is defined by r_min, r_max, phi_min, phi_max.
1034        The number of bin in phi also has to be defined.
1035    """
1036    def __call__(self, data2D):
1037        """
1038            Perform sector average and return I(phi).
1039           
1040            @param data2D: Data2D object
1041            @return: Data1D object
1042        """
1043        return self._agv(data2D, 'phi')
1044
1045class SectorQold(_Sector):
1046    """
1047        Sector average as a function of Q.
1048        I(Q) is return and the data is averaged over phi.
1049       
1050        A sector is defined by r_min, r_max, phi_min, phi_max.
1051        The number of bin in Q also has to be defined.
1052    """
1053    def __call__(self, data2D):
1054        """
1055            Perform sector average and return I(Q).
1056           
1057            @param data2D: Data2D object
1058            @return: Data1D object
1059        """
1060        return self._agv(data2D, 'q')
1061   
1062class SectorQ(_Sector):
1063    """
1064        Sector average as a function of Q for both symatric wings.
1065        I(Q) is return and the data is averaged over phi.
1066       
1067        A sector is defined by r_min, r_max, phi_min, phi_max.
1068        r_min, r_max, phi_min, phi_max >0. 
1069        The number of bin in Q also has to be defined.
1070    """
1071    def __call__(self, data2D):
1072        """
1073            Perform sector average and return I(Q).
1074           
1075            @param data2D: Data2D object
1076            @return: Data1D object
1077        """
1078        return self._agv(data2D, 'q2')
1079if __name__ == "__main__": 
1080
1081    from loader import Loader
1082   
1083
1084    d = Loader().load('test/MAR07232_rest.ASC')
1085    #d = Loader().load('test/MP_New.sans')
1086
1087   
1088    r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
1089    o = r(d)
1090   
1091    s = Ring(r_min=.000001, r_max=.01) 
1092    p = s(d)
1093   
1094    for i in range(len(o.x)):
1095        print o.x[i], o.y[i], o.dy[i]
1096   
1097 
1098   
Note: See TracBrowser for help on using the repository browser.