source: sasview/DataLoader/manipulations.py @ 935ae18

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 935ae18 was f265927, checked in by Jae Cho <jhjcho@…>, 15 years ago

mask feature added

  • Property mode set to 100644
File size: 34.9 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
53def flip_phi(phi):
54    """
55        Correct phi to within the 0 <= to <= 2pi range
56        @return: phi in >=0 and <=2Pi
57    """
58    Pi = math.pi
59    if phi < 0:
60        phi_out = phi  + 2*Pi
61    elif phi > 2*Pi:
62        phi_out = phi  - 2*Pi
63    else:
64        phi_out = phi
65    return phi_out
66
67def reader2D_converter(data2d=None):
68    """
69        convert old 2d format opened by IhorReader or danse_reader to new Data2D format
70        @param data2d: 2d array of Data2D object
71        @return: 1d arrays of Data2D object
72    """
73    if data2d.data==None or data2d.x_bins==None or data2d.y_bins==None:
74        raise ValueError,"Can't convert this data: data=None..."
75   
76    from DataLoader.data_info import Data2D
77
78    new_x = numpy.tile(data2d.x_bins, (len(data2d.y_bins),1))
79    new_y = numpy.tile(data2d.y_bins, (len(data2d.x_bins),1))
80    new_y = new_y.swapaxes(0,1)
81
82    new_data = data2d.data.flatten()
83    qx_data = new_x.flatten()
84    qy_data = new_y.flatten()
85    q_data = numpy.sqrt(qx_data*qx_data+qy_data*qy_data)
86    if data2d.err_data == None or numpy.any(data2d.err_data<=0): 
87        new_err_data = numpy.sqrt(numpy.abs(new_data))
88    else:
89        new_err_data = data2d.err_data.flatten()
90    mask    = numpy.ones(len(new_data), dtype = bool)
91
92    output = Data2D()
93    output = data2d
94    output.data = new_data
95    output.err_data = new_err_data
96    output.qx_data = qx_data
97    output.qy_data = qy_data
98    output.q_data = q_data
99    output.mask = mask
100
101    return output
102
103class _Slab(object):
104    """
105        Compute average I(Q) for a region of interest
106    """
107    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0, bin_width=0.001):
108        # Minimum Qx value [A-1]
109        self.x_min = x_min
110        # Maximum Qx value [A-1]
111        self.x_max = x_max
112        # Minimum Qy value [A-1]
113        self.y_min = y_min
114        # Maximum Qy value [A-1]
115        self.y_max = y_max
116        # Bin width (step size) [A-1]
117        self.bin_width = bin_width
118        # If True, I(|Q|) will be return, otherwise, negative q-values are allowed
119        self.fold = False
120       
121    def __call__(self, data2D): return NotImplemented
122       
123    def _avg(self, data2D, maj):
124        """
125             Compute average I(Q_maj) for a region of interest.
126             The major axis is defined as the axis of Q_maj.
127             The minor axis is the axis that we average over.
128             
129             @param data2D: Data2D object
130             @param maj_min: min value on the major axis
131             @return: Data1D object
132        """
133        if len(data2D.detector) != 1:
134            raise RuntimeError, "_Slab._avg: invalid number of detectors: %g" % len(data2D.detector)
135       
136        # Get data
137        data = data2D.data[numpy.isfinite(data2D.data)]
138        q_data = data2D.q_data[numpy.isfinite(data2D.data)]
139        err_data = data2D.err_data[numpy.isfinite(data2D.data)]
140        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
141        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
142             
143        # Build array of Q intervals
144        if maj=='x':
145            if self.fold: x_min = 0     
146            else:  x_min = self.x_min
147            nbins = int(math.ceil((self.x_max-x_min)/self.bin_width))
148            qbins = self.bin_width*numpy.arange(nbins)+ x_min
149        elif maj=='y':
150            if self.fold: y_min = 0     
151            else:  y_min = self.y_min
152            nbins = int(math.ceil((self.y_max-y_min)/self.bin_width))
153            qbins = self.bin_width*numpy.arange(nbins)+ y_min         
154        else:
155            raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj)
156                               
157        x  = numpy.zeros(nbins)
158        y  = numpy.zeros(nbins)
159        err_y = numpy.zeros(nbins)
160        y_counts = numpy.zeros(nbins)
161
162        # Average pixelsize in q space                                               
163        for npts in range(len(data)): 
164            # default frac                                               
165            frac_x = 0
166            frac_y = 0
167            # get ROI
168            if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]:
169                frac_x = 1
170            if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]:
171                frac_y = 1
172           
173            frac = frac_x * frac_y
174           
175            if frac == 0: continue
176
177            # binning: find axis of q
178            if maj=='x': 
179                q_value = qx_data[npts]
180                min = x_min
181            if maj=='y': 
182                q_value = qy_data[npts] 
183                min = y_min
184            if self.fold and q_value<0: q_value = -q_value   
185           
186            # bin
187            i_q = int(math.ceil((q_value-min)/self.bin_width)) - 1
188           
189            # skip outside of max bins
190            if i_q<0 or i_q>=nbins: continue
191           
192            # give it full weight
193            #frac = 1
194           
195            #TODO: find better definition of x[i_q] based on q_data
196            x[i_q]         = min +(i_q+1)*self.bin_width/2.0
197            y[i_q]         += frac * data[npts]
198           
199            if err_data == None or err_data[npts]==0.0:
200                if data[npts] <0: data[npts] = -data[npts]
201                err_y[i_q] += frac * frac * data[npts]
202            else:
203                err_y[i_q] += frac * frac * err_data[npts] * err_data[npts]
204            y_counts[i_q]  += frac
205               
206        # Average the sums       
207        for n in range(nbins):
208            err_y[n] = math.sqrt(err_y[n])
209         
210        err_y = err_y/y_counts
211        y    = y/y_counts
212       
213        idx = (numpy.isfinite(y)& numpy.isfinite(x)) 
214       
215        if not idx.any(): 
216            raise ValueError, "Average Error: No points inside ROI to average..." 
217        elif len(y[idx])!= nbins:
218            print "resulted",nbins- len(y[idx]),"empty bin(s) due to tight binning..."
219        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx])
220       
221class SlabY(_Slab):
222    """
223        Compute average I(Qy) for a region of interest
224    """
225    def __call__(self, data2D):
226        """
227             Compute average I(Qy) for a region of interest
228             
229             @param data2D: Data2D object
230             @return: Data1D object
231        """
232        return self._avg(data2D, 'y')
233       
234class SlabX(_Slab):
235    """
236        Compute average I(Qx) for a region of interest
237    """
238    def __call__(self, data2D):
239        """
240             Compute average I(Qx) for a region of interest
241             
242             @param data2D: Data2D object
243             @return: Data1D object
244        """
245        return self._avg(data2D, 'x') 
246       
247class Boxsum(object):
248    """
249        Perform the sum of counts in a 2D region of interest.
250    """
251    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
252        # Minimum Qx value [A-1]
253        self.x_min = x_min
254        # Maximum Qx value [A-1]
255        self.x_max = x_max
256        # Minimum Qy value [A-1]
257        self.y_min = y_min
258        # Maximum Qy value [A-1]
259        self.y_max = y_max
260
261    def __call__(self, data2D):
262        """
263             Perform the sum in the region of interest
264             
265             @param data2D: Data2D object
266             @return: number of counts, error on number of counts
267        """
268        y, err_y, y_counts = self._sum(data2D)
269       
270        # Average the sums
271        counts = 0 if y_counts==0 else y
272        error  = 0 if y_counts==0 else math.sqrt(err_y)
273       
274        return counts, error
275       
276    def _sum(self, data2D):
277        """
278             Perform the sum in the region of interest
279             @param data2D: Data2D object
280             @return: number of counts, error on number of counts, number of entries summed
281        """
282        if len(data2D.detector) != 1:
283            raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector)
284       
285        # Get data
286        data = data2D.data[numpy.isfinite(data2D.data)]
287        q_data = data2D.q_data[numpy.isfinite(data2D.data)]
288        err_data = data2D.err_data[numpy.isfinite(data2D.data)]
289        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
290        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
291   
292        y  = 0.0
293        err_y = 0.0
294        y_counts = 0.0
295
296        # Average pixelsize in q space                                               
297        for npts in range(len(data)): 
298            # default frac                                               
299            frac_x = 0 
300            frac_y = 0                 
301           
302            # get min and max at each points
303            qx = qx_data[npts]
304            qy = qy_data[npts]
305           
306            # get the ROI
307            if self.x_min <= qx and self.x_max > qx:
308                frac_x = 1
309            if self.y_min <= qy and self.y_max > qy:
310                frac_y = 1
311               
312            #Find the fraction along each directions           
313            frac = frac_x * frac_y
314            if frac == 0: continue
315
316            y += frac * data[npts]
317            if err_data == None or err_data[npts]==0.0:
318                if data[npts] <0: data[npts] = -data[npts]
319                err_y += frac * frac * data[npts]
320            else:
321                err_y += frac * frac * err_data[npts] * err_data[npts]
322            y_counts += frac
323               
324        return y, err_y, y_counts
325
326
327     
328class Boxavg(Boxsum):
329    """
330        Perform the average of counts in a 2D region of interest.
331    """
332    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
333        super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)
334
335    def __call__(self, data2D):
336        """
337             Perform the sum in the region of interest
338             
339             @param data2D: Data2D object
340             @return: average counts, error on average counts
341        """
342        y, err_y, y_counts = self._sum(data2D)
343       
344        # Average the sums
345        counts = 0 if y_counts==0 else y/y_counts
346        error  = 0 if y_counts==0 else math.sqrt(err_y)/y_counts
347       
348        return counts, error
349       
350def get_pixel_fraction_square(x, xmin, xmax):
351    """
352         Return the fraction of the length
353         from xmin to x.
354         
355             A            B
356         +-----------+---------+
357         xmin        x         xmax
358         
359         @param x: x-value
360         @param xmin: minimum x for the length considered
361         @param xmax: minimum x for the length considered
362         @return: (x-xmin)/(xmax-xmin) when xmin < x < xmax
363         
364    """
365    if x<=xmin:
366        return 0.0
367    if x>xmin and x<xmax:
368        return (x-xmin)/(xmax-xmin)
369    else:
370        return 1.0
371
372
373class CircularAverage(object):
374    """
375        Perform circular averaging on 2D data
376       
377        The data returned is the distribution of counts
378        as a function of Q
379    """
380    def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005):
381        # Minimum radius included in the average [A-1]
382        self.r_min = r_min
383        # Maximum radius included in the average [A-1]
384        self.r_max = r_max
385        # Bin width (step size) [A-1]
386        self.bin_width = bin_width
387
388    def __call__(self, data2D):
389        """
390            Perform circular averaging on the data
391           
392            @param data2D: Data2D object
393            @return: Data1D object
394        """
395        # Get data
396        data = data2D.data[numpy.isfinite(data2D.data)]
397        q_data = data2D.q_data[numpy.isfinite(data2D.data)]
398        err_data = data2D.err_data[numpy.isfinite(data2D.data)]
399   
400        q_data_max = numpy.max(q_data)
401
402        if len(data2D.q_data) == None:
403            raise RuntimeError, "Circular averaging: invalid q_data: %g" % data2D.q_data
404
405        # Build array of Q intervals
406        nbins = int(math.ceil((self.r_max-self.r_min)/self.bin_width))
407        qbins = self.bin_width*numpy.arange(nbins)+self.r_min
408
409        x  = numpy.zeros(nbins)
410        y  = numpy.zeros(nbins)
411        err_y = numpy.zeros(nbins)
412        y_counts = numpy.zeros(nbins)
413
414        for npt in range(len(data)):         
415            frac = 0
416           
417            # q-value at the pixel (j,i)
418            q_value = q_data[npt]
419               
420            data_n = data[npt]               
421           
422            ## No need to calculate the frac when all data are within range
423            if self.r_min >= self.r_max:
424                raise ValueError, "Limit Error: min > max ???" 
425           
426            if self.r_min <= q_value and q_value <= self.r_max: frac = 1                       
427
428            if frac == 0: continue
429                       
430            i_q = int(math.floor((q_value-self.r_min)/self.bin_width))   
431
432            # Take care of the edge case at phi = 2pi.
433            if i_q == nbins: 
434                i_q = nbins -1 
435                                   
436            y[i_q] += frac * data_n
437
438            if err_data == None or err_data[npt]==0.0:
439                if data_n <0: data_n = -data_n
440                err_y[i_q] += frac * frac * data_n
441            else:
442                err_y[i_q] += frac * frac * err_data[npt] * err_data[npt]
443            y_counts[i_q]  += frac
444       
445        ## x should be the center value of each bins       
446        x = qbins+self.bin_width/2 
447       
448        # Average the sums       
449        for n in range(nbins):
450            if err_y[n] <0: err_y[n] = -err_y[n]
451            err_y[n] = math.sqrt(err_y[n])
452                     
453        err_y = err_y/y_counts
454        y    = y/y_counts
455        idx = (numpy.isfinite(y))&(numpy.isfinite(x)) 
456       
457        if not idx.any(): 
458            raise ValueError, "Average Error: No points inside ROI to average..." 
459        elif len(y[idx])!= nbins:
460            print "resulted",nbins- len(y[idx]),"empty bin(s) due to tight binning..."
461       
462        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx])
463   
464
465class Ring(object):
466    """
467        Defines a ring on a 2D data set.
468        The ring is defined by r_min, r_max, and
469        the position of the center of the ring.
470       
471        The data returned is the distribution of counts
472        around the ring as a function of phi.
473       
474        Phi_min and phi_max should be defined between 0 and 2*pi
475        in anti-clockwise starting from the x- axis on the left-hand side
476    """
477    #Todo: remove center.
478    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0,nbins=20 ):
479        # Minimum radius
480        self.r_min = r_min
481        # Maximum radius
482        self.r_max = r_max
483        # Center of the ring in x
484        self.center_x = center_x
485        # Center of the ring in y
486        self.center_y = center_y
487        # Number of angular bins
488        self.nbins_phi = nbins
489       
490    def __call__(self, data2D):
491        """
492            Apply the ring to the data set.
493            Returns the angular distribution for a given q range
494           
495            @param data2D: Data2D object
496            @return: Data1D object
497        """
498        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
499            raise RuntimeError, "Ring averaging only take plottable_2D objects"
500       
501        Pi = math.pi
502 
503        # Get data
504        data = data2D.data[numpy.isfinite(data2D.data)]
505        q_data = data2D.q_data[numpy.isfinite(data2D.data)]
506        err_data = data2D.err_data[numpy.isfinite(data2D.data)]
507        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
508        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
509       
510        q_data_max = numpy.max(q_data)
511       
512        # Set space for 1d outputs
513        phi_bins   = numpy.zeros(self.nbins_phi)
514        phi_counts = numpy.zeros(self.nbins_phi)
515        phi_values = numpy.zeros(self.nbins_phi)
516        phi_err    = numpy.zeros(self.nbins_phi)
517       
518        for npt in range(len(data)):         
519            frac = 0
520           
521            # q-value at the point (npt)
522            q_value = q_data[npt]
523               
524            data_n = data[npt]   
525                       
526            # phi-value at the point (npt)
527            phi_value=math.atan2(qy_data[npt],qx_data[npt])+Pi
528           
529            if self.r_min <= q_value and q_value <= self.r_max: frac = 1                       
530
531            if frac == 0: continue
532           
533            # binning           
534            i_phi = int(math.floor((self.nbins_phi)*phi_value/(2*Pi)))
535           
536            # Take care of the edge case at phi = 2pi.
537            if i_phi == self.nbins_phi: 
538                i_phi =  self.nbins_phi -1 
539                                             
540            phi_bins[i_phi] += frac * data[npt]
541           
542            if err_data == None or err_data[npt] ==0.0:
543                if data_n <0: data_n = -data_n
544                phi_err[i_phi] += frac * frac * math.fabs(data_n)
545            else:
546                phi_err[i_phi] += frac * frac *err_data[npt]*err_data[npt]
547            phi_counts[i_phi] += frac
548                         
549        for i in range(self.nbins_phi):
550            phi_bins[i] = phi_bins[i] / phi_counts[i]
551            phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i]
552            phi_values[i] = 2.0*math.pi/self.nbins_phi*(1.0*i + 0.5)
553           
554        idx = (numpy.isfinite(phi_bins)) 
555
556        if not idx.any(): 
557            raise ValueError, "Average Error: No points inside ROI to average..." 
558        elif len(phi_bins[idx])!= self.nbins_phi:
559            print "resulted",self.nbins_phi- len(phi_bins[idx]),"empty bin(s) due to tight binning..."
560        return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx])
561   
562def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11):
563    """
564        Returns the fraction of the pixel defined by
565        the four corners (q_00, q_01, q_10, q_11) that
566        has q < qmax.
567       
568                q_01                q_11
569        y=1         +--------------+
570                    |              |
571                    |              |
572                    |              |
573        y=0         +--------------+
574                q_00                q_10
575       
576                    x=0            x=1
577       
578    """
579   
580    # y side for x = minx
581    x_0 = get_intercept(qmax, q_00, q_01)
582    # y side for x = maxx
583    x_1 = get_intercept(qmax, q_10, q_11)
584   
585    # x side for y = miny
586    y_0 = get_intercept(qmax, q_00, q_10)
587    # x side for y = maxy
588    y_1 = get_intercept(qmax, q_01, q_11)
589   
590    # surface fraction for a 1x1 pixel
591    frac_max = 0
592   
593    if x_0 and x_1:
594        frac_max = (x_0+x_1)/2.0
595   
596    elif y_0 and y_1:
597        frac_max = (y_0+y_1)/2.0
598   
599    elif x_0 and y_0:
600        if q_00 < q_10:
601            frac_max = x_0*y_0/2.0
602        else:
603            frac_max = 1.0-x_0*y_0/2.0
604   
605    elif x_0 and y_1:
606        if q_00 < q_10:
607            frac_max = x_0*y_1/2.0
608        else:
609            frac_max = 1.0-x_0*y_1/2.0
610   
611    elif x_1 and y_0:
612        if q_00 > q_10:
613            frac_max = x_1*y_0/2.0
614        else:
615            frac_max = 1.0-x_1*y_0/2.0
616   
617    elif x_1 and y_1:
618        if q_00 < q_10:
619            frac_max = 1.0 - (1.0-x_1)*(1.0-y_1)/2.0
620        else:
621            frac_max = (1.0-x_1)*(1.0-y_1)/2.0
622           
623    # If we make it here, there is no intercept between
624    # this pixel and the constant-q ring. We only need
625    # to know if we have to include it or exclude it.
626    elif (q_00+q_01+q_10+q_11)/4.0 < qmax:
627        frac_max = 1.0
628
629    return frac_max
630             
631def get_intercept(q, q_0, q_1):
632    """
633        Returns the fraction of the side at which the
634        q-value intercept the pixel, None otherwise.
635        The values returned is the fraction ON THE SIDE
636        OF THE LOWEST Q.
637       
638       
639       
640                A        B   
641         +-----------+--------+
642         0                    1     <--- pixel size
643         
644        Q_0 -------- Q ----- Q_1    <--- equivalent Q range
645       
646       
647        if Q_1 > Q_0, A is returned
648        if Q_1 < Q_0, B is returned
649       
650        if Q is outside the range of [Q_0, Q_1], None is returned
651         
652    """
653    if q_1 > q_0:
654        if (q > q_0 and q <= q_1):
655            return (q-q_0)/(q_1 - q_0)   
656    else:
657        if (q > q_1 and q <= q_0):
658            return (q-q_1)/(q_0 - q_1)
659    return None
660     
661class _Sector:
662    """
663        Defines a sector region on a 2D data set.
664        The sector is defined by r_min, r_max, phi_min, phi_max,
665        and the position of the center of the ring
666        where phi_min and phi_max are defined by the right and left lines wrt central line
667        and phi_max could be less than phi_min.
668       
669        Phi is defined between 0 and 2*pi in anti-clockwise starting from the x- axis on the left-hand side
670    """
671    def __init__(self, r_min, r_max, phi_min=0, phi_max=2*math.pi,nbins=20):
672        self.r_min = r_min
673        self.r_max = r_max
674        self.phi_min = phi_min
675        self.phi_max = phi_max
676        self.nbins = nbins
677       
678   
679    def _agv(self, data2D, run='phi'):
680        """
681            Perform sector averaging.
682           
683            @param data2D: Data2D object
684            @param run:  define the varying parameter ('phi' , 'q' , or 'q2')
685            @return: Data1D object
686        """
687        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
688            raise RuntimeError, "Ring averaging only take plottable_2D objects"
689        Pi = math.pi
690
691        # Get the all data & info
692        data = data2D.data[numpy.isfinite(data2D.data)]
693        q_data = data2D.q_data[numpy.isfinite(data2D.data)]
694        err_data = data2D.err_data[numpy.isfinite(data2D.data)]
695        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
696        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)]
697       
698        #set space for 1d outputs
699        x        = numpy.zeros(self.nbins)
700        y        = numpy.zeros(self.nbins)
701        y_err    = numpy.zeros(self.nbins)         
702        y_counts = numpy.zeros(self.nbins)
703                     
704        # Get the min and max into the region: 0 <= phi < 2Pi
705        phi_min = flip_phi(self.phi_min)
706        phi_max = flip_phi(self.phi_max)
707       
708        q_data_max = numpy.max(q_data)
709                     
710        for n in range(len(data)):     
711                frac = 0
712               
713                # q-value at the pixel (j,i)
714                q_value = q_data[n]
715               
716   
717                data_n = data[n]
718               
719                # Is pixel within range?
720                is_in = False
721               
722                # phi-value of the pixel (j,i)
723                phi_value=math.atan2(qy_data[n],qx_data[n])+Pi
724               
725                ## No need to calculate the frac when all data are within range
726                if self.r_min <= q_value and q_value <= self.r_max: frac = 1                       
727
728                if frac == 0: continue
729 
730                #In case of two ROIs (symmetric major and minor regions)(for 'q2')
731                if run.lower()=='q2':
732                    ## For minor sector wing
733                    # Calculate the minor wing phis
734                    phi_min_minor = flip_phi(phi_min-Pi)
735                    phi_max_minor = flip_phi(phi_max-Pi)
736                    # Check if phis of the minor ring is within 0 to 2pi
737                    if phi_min_minor > phi_max_minor:
738                        is_in = (phi_value > phi_min_minor or phi_value < phi_max_minor)
739                    else:
740                        is_in = (phi_value > phi_min_minor and phi_value < phi_max_minor)
741
742                #For all cases(i.e.,for 'q', 'q2', and 'phi')
743                #Find pixels within ROI 
744                if phi_min > phi_max: 
745                    is_in = is_in or (phi_value > phi_min or phi_value < phi_max)         
746                else:
747                    is_in = is_in or (phi_value>= phi_min  and  phi_value <phi_max)
748               
749                if not is_in: frac = 0                                                       
750                if frac == 0: continue
751               
752                # Check which type of averaging we need
753                if run.lower()=='phi': 
754                    i_bin    = int(math.floor((self.nbins)*(phi_value-self.phi_min)\
755                                              /(self.phi_max-self.phi_min)))
756                else:
757                    i_bin = int(math.floor((self.nbins)*(q_value-self.r_min)/(self.r_max-self.r_min)))
758
759                # Take care of the edge case at phi = 2pi.
760                if i_bin == self.nbins: 
761                    i_bin =  self.nbins -1
762                   
763                ## Get the total y         
764                y[i_bin] += frac * data_n
765
766                if err_data == None or err_data[n] ==0.0:
767                    if data_n<0: data_n= -data_n
768                    y_err[i_bin] += frac * frac * data_n
769                else:
770                    y_err[i_bin] += frac * frac * err_data[n]*err_data[n]
771                y_counts[i_bin] += frac
772       
773        # Organize the results
774        for i in range(self.nbins):
775            y[i] = y[i] / y_counts[i]
776            y_err[i] = math.sqrt(y_err[i]) / y_counts[i]
777           
778            # The type of averaging: phi,q2, or q
779            # Calculate x[i]should be at the center of the bin
780            if run.lower()=='phi':               
781                x[i] = (self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min
782            else:
783                x[i] = (self.r_max-self.r_min)/self.nbins*(1.0*i + 0.5)+self.r_min
784               
785        idx = (numpy.isfinite(y)& numpy.isfinite(y_err))
786       
787        if not idx.any(): 
788            raise ValueError, "Average Error: No points inside sector of ROI to average..." 
789        elif len(y[idx])!= self.nbins:
790            print "resulted",self.nbins- len(y[idx]),"empty bin(s) due to tight binning..."
791        return Data1D(x=x[idx], y=y[idx], dy=y_err[idx])
792               
793class SectorPhi(_Sector):
794    """
795        Sector average as a function of phi.
796        I(phi) is return and the data is averaged over Q.
797       
798        A sector is defined by r_min, r_max, phi_min, phi_max.
799        The number of bin in phi also has to be defined.
800    """
801    def __call__(self, data2D):
802        """
803            Perform sector average and return I(phi).
804           
805            @param data2D: Data2D object
806            @return: Data1D object
807        """
808
809        return self._agv(data2D, 'phi')
810   
811class SectorQ(_Sector):
812    """
813        Sector average as a function of Q for both symatric wings.
814        I(Q) is return and the data is averaged over phi.
815       
816        A sector is defined by r_min, r_max, phi_min, phi_max.
817        r_min, r_max, phi_min, phi_max >0. 
818        The number of bin in Q also has to be defined.
819    """
820    def __call__(self, data2D):
821        """
822            Perform sector average and return I(Q).
823           
824            @param data2D: Data2D object
825            @return: Data1D object
826        """
827        return self._agv(data2D, 'q2')
828
829class Ringcut(object):
830    """
831        Defines a ring on a 2D data set.
832        The ring is defined by r_min, r_max, and
833        the position of the center of the ring.
834       
835        The data returned is the region inside the ring
836       
837        Phi_min and phi_max should be defined between 0 and 2*pi
838        in anti-clockwise starting from the x- axis on the left-hand side
839    """
840    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0 ):
841        # Minimum radius
842        self.r_min = r_min
843        # Maximum radius
844        self.r_max = r_max
845        # Center of the ring in x
846        self.center_x = center_x
847        # Center of the ring in y
848        self.center_y = center_y
849
850       
851    def __call__(self, data2D):
852        """
853            Apply the ring to the data set.
854            Returns the angular distribution for a given q range
855           
856            @param data2D: Data2D object
857            @return: index array in the range
858        """
859        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
860            raise RuntimeError, "Ring cut only take plottable_2D objects"
861
862        # Get data
863        qx_data = data2D.qx_data
864        qy_data = data2D.qy_data
865        mask = data2D.mask
866        q_data = numpy.sqrt(qx_data*qx_data+qy_data*qy_data)
867        #q_data_max = numpy.max(q_data)
868
869        # check whether or not the data point is inside ROI
870        out = (self.r_min <= q_data) & (self.r_max >= q_data)
871
872        return (out)
873       
874
875class Boxcut(object):
876    """
877        Find a rectangular 2D region of interest.
878    """
879    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0):
880        # Minimum Qx value [A-1]
881        self.x_min = x_min
882        # Maximum Qx value [A-1]
883        self.x_max = x_max
884        # Minimum Qy value [A-1]
885        self.y_min = y_min
886        # Maximum Qy value [A-1]
887        self.y_max = y_max
888
889    def __call__(self, data2D):
890        """
891           Find a rectangular 2D region of interest.
892             
893           @param data2D: Data2D object
894           @return: mask, 1d array (len = len(data))
895               with Trues where the data points are inside ROI, otherwise False
896        """
897        mask = self._find(data2D)
898       
899        return mask
900       
901    def _find(self, data2D):
902        """
903             Find a rectangular 2D region of interest.
904             @param data2D: Data2D object
905             @return: out, 1d array (length = len(data))
906               with Trues where the data points are inside ROI, otherwise Falses
907        """
908        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
909            raise RuntimeError, "Boxcut take only plottable_2D objects"
910        # Get qx_ and qy_data
911        qx_data = data2D.qx_data
912        qy_data = data2D.qy_data
913        mask = data2D.mask
914       
915        # check whether or not the data point is inside ROI
916        outx = (self.x_min <= qx_data) & (self.x_max > qx_data)
917        outy = (self.y_min <= qy_data) & (self.y_max > qy_data)
918
919        return (outx & outy)
920
921class Sectorcut(object):
922    """
923        Defines a sector (major + minor) region on a 2D data set.
924        The sector is defined by phi_min, phi_max,
925        where phi_min and phi_max are defined by the right and left lines wrt central line.
926       
927        Phi_min and phi_max are given in units of radian
928        and (phi_max-phi_min) should not be larger than pi
929    """
930    def __init__(self,phi_min=0, phi_max=math.pi):
931        self.phi_min = phi_min
932        self.phi_max = phi_max
933             
934    def __call__(self, data2D):
935        """
936           Find a rectangular 2D region of interest.
937           
938           @param data2D: Data2D object
939           @return: mask, 1d array (len = len(data))
940               with Trues where the data points are inside ROI, otherwise False
941        """
942        mask = self._find(data2D)
943       
944        return mask
945       
946    def _find(self, data2D):
947        """
948             Find a rectangular 2D region of interest.
949             @param data2D: Data2D object
950             @return: out, 1d array (length = len(data))
951               with Trues where the data points are inside ROI, otherwise Falses
952        """
953        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:
954            raise RuntimeError, "Sectorcut take only plottable_2D objects" 
955        Pi = math.pi
956        # Get data
957        qx_data = data2D.qx_data
958        qy_data = data2D.qy_data       
959        phi_data = numpy.zeros(len(qx_data))
960
961        # get phi from data
962        phi_data = numpy.arctan2(qy_data, qx_data)
963       
964        # Get the min and max into the region: -pi <= phi < Pi
965        phi_min_major = flip_phi(self.phi_min+Pi)-Pi
966        phi_max_major = flip_phi(self.phi_max+Pi)-Pi 
967        # check for major sector
968        if phi_min_major > phi_max_major:
969            out_major = (phi_min_major <= phi_data) + (phi_max_major > phi_data)
970        else:
971            out_major = (phi_min_major <= phi_data) & (phi_max_major > phi_data)
972         
973        # minor sector
974        # Get the min and max into the region: -pi <= phi < Pi
975        phi_min_minor = flip_phi(self.phi_min)-Pi
976        phi_max_minor = flip_phi(self.phi_max)-Pi
977             
978        # check for minor sector
979        if phi_min_minor > phi_max_minor:
980            out_minor= (phi_min_minor <= phi_data) + (phi_max_minor>= phi_data) 
981        else:
982            out_minor = (phi_min_minor <= phi_data) & (phi_max_minor >= phi_data) 
983        out = out_major + out_minor
984       
985        return out
986
987if __name__ == "__main__": 
988
989    from loader import Loader
990   
991
992    d = Loader().load('test/MAR07232_rest.ASC')
993    #d = Loader().load('test/MP_New.sans')
994
995   
996    r = SectorQ(r_min=.000001, r_max=.01, phi_min=0.0, phi_max=2*math.pi)
997    o = r(d)
998   
999    s = Ring(r_min=.000001, r_max=.01) 
1000    p = s(d)
1001   
1002    for i in range(len(o.x)):
1003        print o.x[i], o.y[i], o.dy[i]
1004   
1005 
1006   
Note: See TracBrowser for help on using the repository browser.