Changeset 3c67340 in sasview


Ignore:
Timestamp:
Jan 25, 2010 6:22:03 PM (14 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
393f0f3
Parents:
2f569b3
Message:

dataloader: revert changes to manipulation until the issues are properly investigated

File:
1 edited

Legend:

Unmodified
Added
Removed
  • DataLoader/manipulations.py

    r2f569b3 r3c67340  
    104104        y  = numpy.zeros(nbins) 
    105105        err_y = numpy.zeros(nbins) 
    106         x_counts = numpy.zeros(nbins) 
    107106        y_counts = numpy.zeros(nbins) 
    108107                                                 
     
    161160                        continue 
    162161             
    163                 x[i_q]          += q_value 
     162                x[i_q]          = q_value 
    164163                y[i_q]         += frac * data2D.data[j][i] 
    165164                if data2D.err_data == None or data2D.err_data[j][i]==0.0: 
     
    167166                else: 
    168167                    err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 
    169                 x_counts[i_q] += 1 
    170168                y_counts[i_q]  += frac 
    171169                 
    172170        # Average the sums 
    173171        for i in range(nbins): 
    174             if x_counts[i]>0 and y_counts[i]>0: 
     172            if y_counts[i]>0: 
    175173                err_y[i] = math.sqrt(err_y[i])/y_counts[i] 
    176                 x[i]     = x[i]/x_counts[i] 
    177174                y[i]     = y[i]/y_counts[i] 
    178175         
     
    359356        as a function of Q 
    360357    """ 
    361     def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 
     358    def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.001): 
    362359        # Minimum radius included in the average [A-1] 
    363360        self.r_min = r_min 
     
    381378        det_dist    = data2D.detector[0].distance 
    382379        wavelength  = data2D.source.wavelength 
    383         center_x    = data2D.detector[0].beam_center.x/pixel_width_x 
    384         center_y    = data2D.detector[0].beam_center.y/pixel_width_y 
     380        center_x    = data2D.detector[0].beam_center.x/pixel_width_x+0.5 
     381        center_y    = data2D.detector[0].beam_center.y/pixel_width_y+0.5 
    385382         
    386383        # Find out the maximum Q range 
     
    396393         
    397394        qmax = get_q(dx_max, dy_max, det_dist, wavelength) 
    398  
     395         
    399396        # Build array of Q intervals 
    400397        nbins = int(math.ceil((qmax-self.r_min)/self.bin_width)) 
     
    402399         
    403400        x  = numpy.zeros(nbins) 
    404         x_counts  = numpy.zeros(nbins) 
    405401        y  = numpy.zeros(nbins) 
    406402        err_y = numpy.zeros(nbins) 
    407403        y_counts = numpy.zeros(nbins) 
    408  
     404         
    409405        for i in range(numpy.size(data2D.data,1)): 
    410              
     406            dx = pixel_width_x*(i+0.5 - center_x) 
    411407             
    412408            # Min and max x-value for the pixel 
     
    415411             
    416412            for j in range(numpy.size(data2D.data,0)): 
    417                  
     413                dy = pixel_width_y*(j+0.5 - center_y) 
     414             
     415                q_value = get_q(dx, dy, det_dist, wavelength) 
    418416             
    419417                # Min and max y-value for the pixel 
     
    428426                q_11 = get_q(maxx, maxy, det_dist, wavelength) 
    429427                 
    430                 #number of sub-bin: default = 1 (ie., 1 bin: no sub-bin) 
    431                 nsubbins = 2 
    432  
    433                 #if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \ 
    434                 #    or is_intercept(self.r_min, q_00, q_01, q_10, q_11): 
    435                 #    # Number of sub-bins near the boundary of ROI in x-direction 
    436                 #    nsubbins = 3 
    437                  
    438                 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 
    439                 for x_subbins in range(nsubbins): 
    440                     for y_subbins in range(nsubbins): 
    441                          
    442                         dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 
    443                         dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 
    444                          
    445                         q_value = get_q(dx, dy, det_dist, wavelength) 
    446                         # Min and max x-value for the subpixel 
    447                         minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x) 
    448                         maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x) 
    449                         # Min and max y-value for the subpixel 
    450                         miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y) 
    451                         maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y) 
    452                         # Look for intercept between each side of the pixel 
    453                         # and the constant-q ring for qmax 
    454                         # Calculate the q-value for each corner of sub pixel 
    455                         # q_[x min or max][y min or max] 
    456                         q_00 = get_q(minx, miny, det_dist, wavelength) 
    457                         q_01 = get_q(minx, maxy, det_dist, wavelength) 
    458                         q_10 = get_q(maxx, miny, det_dist, wavelength) 
    459                         q_11 = get_q(maxx, maxy, det_dist, wavelength) 
    460                         frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11)/(nsubbins*nsubbins) 
    461                  
    462                         # Look for intercept between each side of the pixel 
    463                         # and the constant-q ring for qmin 
    464                         frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11)/(nsubbins*nsubbins) 
    465                  
    466                         # We are interested in the region between qmin and qmax 
    467                         # therefore the fraction of the surface of the pixel 
    468                         # that we will use to calculate the number of counts to  
    469                         # include is given by: 
    470                         frac = (frac_max - frac_min) #??? Todo: What if frac_max - fracmin <0? 
    471                          
    472                         
    473                         i_q = int(math.floor((q_value-self.r_min)/self.bin_width)) #- 1     
    474  
    475                         if q_value > qmax or q_value < self.r_min: 
    476                             continue   
    477                                                  
    478                         x[i_q]         += q_value 
    479                         x_counts[i_q]  += 1 
    480          
    481                         y[i_q]         += frac * data2D.data[j][i] 
    482          
    483                         if data2D.err_data == None or data2D.err_data[j][i]==0.0: 
    484                             err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i]) 
    485                         else: 
    486                             err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 
    487                         y_counts[i_q]  += frac 
    488                                                       
     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                if q_value > qmax or q_value < self.r_min: 
     444                    continue 
     445                     
     446                x[i_q]          = q_value 
     447                y[i_q]         += frac * data2D.data[j][i] 
     448                if data2D.err_data == None or data2D.err_data[j][i]==0.0: 
     449                    err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i]) 
     450                else: 
     451                    err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 
     452                y_counts[i_q]  += frac 
     453                 
    489454        # Average the sums 
    490455        nzero = 0 
    491          
    492456        for i in range(nbins): 
    493             #print x_counts[i],x[i],y[i] 
    494             if y_counts[i]>0 and x_counts[i]>0: 
     457            if y_counts[i]>0: 
    495458                err_y[i] = math.sqrt(err_y[i])/y_counts[i] 
    496459                y[i]     = y[i]/y_counts[i] 
    497                 x[i] = (qmax-self.r_min)/nbins*(1.0*i + 0.5)+ self.r_min#x[i]/x_counts[i] 
    498  
    499460            else: 
    500461                nzero += 1 
     
    505466        j=0 
    506467        for i in range(nbins): 
    507             if err_y[i] != 0 and y_counts[i]>0 and x_counts[i]>0: 
    508                 #Make sure wintin range 
    509                 if x[i] <= self.r_max and x[i] >= self.r_min: 
    510                     tx[j]  = x[i] 
    511                     ty[j]  = y[i] 
    512                     terr_y[j] = err_y[i] 
    513                     j+=1        
     468            if err_y[i] != 0 : 
     469                tx[j]  = x[i] 
     470                ty[j]  = y[i] 
     471                terr_y[j] = err_y[i] 
     472                j+=1 
     473                 
    514474        return Data1D(x=tx, y=ty, dy=terr_y) 
    515475     
     
    691651    elif (q_00+q_01+q_10+q_11)/4.0 < qmax: 
    692652        frac_max = 1.0 
    693  
     653    
    694654    return frac_max 
    695655              
     
    723683            return (q-q_1)/(q_0 - q_1) 
    724684    return None 
    725       
    726 class _Sector: 
     685     
     686#This class can be removed. 
     687class _Sectorold: 
    727688    """ 
    728689        Defines a sector region on a 2D data set. 
    729690        The sector is defined by r_min, r_max, phi_min, phi_max, 
    730         and the position of the center of the ring  
    731         where phi_min and phi_max are defined by the right and left lines wrt central line 
    732         and phi_max could be less than phi_min.  
    733         
     691        and the position of the center of the ring.          
    734692        Phi is defined between 0 and 2pi 
    735693    """ 
     
    768726        y_counts = numpy.zeros(self.nbins) 
    769727        x        = numpy.zeros(self.nbins) 
    770         x_counts = numpy.zeros(self.nbins) 
    771         y_err    = numpy.zeros(self.nbins)        
     728        y_err    = numpy.zeros(self.nbins) 
     729         
     730        for i in range(numpy.size(data,1)): 
     731            dx = pixel_width_x*(i+0.5 - center_x) 
     732             
     733            # Min and max x-value for the pixel 
     734            minx = pixel_width_x*(i - center_x) 
     735            maxx = pixel_width_x*(i+1.0 - center_x) 
     736             
     737            for j in range(numpy.size(data,0)): 
     738                dy = pixel_width_y*(j+0.5 - center_y) 
     739             
     740                q_value = get_q(dx, dy, det_dist, wavelength) 
     741 
     742                # Min and max y-value for the pixel 
     743                miny = pixel_width_y*(j - center_y) 
     744                maxy = pixel_width_y*(j+1.0 - center_y) 
     745                 
     746                # Calculate the q-value for each corner 
     747                # q_[x min or max][y min or max] 
     748                q_00 = get_q(minx, miny, det_dist, wavelength) 
     749                q_01 = get_q(minx, maxy, det_dist, wavelength) 
     750                q_10 = get_q(maxx, miny, det_dist, wavelength) 
     751                q_11 = get_q(maxx, maxy, det_dist, wavelength) 
     752                 
     753                # Look for intercept between each side of the pixel 
     754                # and the constant-q ring for qmax 
     755                frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 
     756                 
     757                # Look for intercept between each side of the pixel 
     758                # and the constant-q ring for qmin 
     759                frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 
     760                 
     761                # We are interested in the region between qmin and qmax 
     762                # therefore the fraction of the surface of the pixel 
     763                # that we will use to calculate the number of counts to  
     764                # include is given by: 
     765                 
     766                frac = frac_max - frac_min 
     767 
     768                # Compute phi and check whether it's within the limits 
     769                phi_value=math.atan2(dy,dx)+math.pi 
     770 #               if phi_value<self.phi_min or phi_value>self.phi_max:                 
     771                if phi_value<self.phi_min or phi_value>self.phi_max: 
     772                    continue 
     773                                                     
     774                # Check which type of averaging we need 
     775                if run.lower()=='phi':  
     776                    i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 
     777                else: 
     778                    # If we don't need this pixel, skip the rest of the work 
     779                    #TODO: an improvement here would be to compute the average 
     780                    # Q for the pixel from the part that is covered by 
     781                    # the ring defined by q_min/q_max rather than the complete 
     782                    # pixel  
     783                    if q_value<self.r_min or q_value>self.r_max: 
     784                        continue 
     785                    i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1 
     786             
     787                try: 
     788                    y[i_bin] += frac * data[j][i] 
     789                except: 
     790                    import sys 
     791                    print sys.exc_value 
     792                    print i_bin, frac 
     793                 
     794                if data2D.err_data == None or data2D.err_data[j][i]==0.0: 
     795                    y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 
     796                else: 
     797                    y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 
     798                y_counts[i_bin] += frac 
     799         
     800        for i in range(self.nbins): 
     801            y[i] = y[i] / y_counts[i] 
     802            y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 
     803            # Check which type of averaging we need 
     804            if run.lower()=='phi': 
     805                x[i] = (self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min 
     806            else: 
     807                x[i] = (self.r_max-self.r_min)/self.nbins*(1.0*i + 0.5)+self.r_min 
     808             
     809        return Data1D(x=x, y=y, dy=y_err) 
     810         
     811class _Sector: 
     812    """ 
     813        Defines a sector region on a 2D data set. 
     814        The sector is defined by r_min, r_max, phi_min, phi_max, 
     815        and the position of the center of the ring  
     816        where phi_min and phi_max are defined by the right and left lines wrt central line 
     817        and phi_max could be less than phi_min.  
     818        
     819        Phi is defined between 0 and 2pi 
     820    """ 
     821    def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20): 
     822        self.r_min = r_min 
     823        self.r_max = r_max 
     824        self.phi_min = phi_min 
     825        self.phi_max = phi_max 
     826        self.nbins = nbins 
     827         
     828    def _agv(self, data2D, run='phi'): 
     829        """ 
     830            Perform sector averaging. 
     831             
     832            @param data2D: Data2D object 
     833            @param run:  define the varying parameter ('phi' or 'q') 
     834            @return: Data1D object 
     835        """ 
     836        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
     837            raise RuntimeError, "Ring averaging only take plottable_2D objects" 
     838                    
     839        data = data2D.data       
     840        qmax = self.r_max 
     841        qmin = self.r_min 
     842         
     843        if len(data2D.detector) != 1: 
     844            raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector) 
     845        pixel_width_x = data2D.detector[0].pixel_size.x 
     846        pixel_width_y = data2D.detector[0].pixel_size.y 
     847        det_dist      = data2D.detector[0].distance 
     848        wavelength    = data2D.source.wavelength 
     849        center_x      = data2D.detector[0].beam_center.x/pixel_width_x 
     850        center_y      = data2D.detector[0].beam_center.y/pixel_width_y 
     851         
     852        y        = numpy.zeros(self.nbins) 
     853        y_counts = numpy.zeros(self.nbins) 
     854        x        = numpy.zeros(self.nbins) 
     855        y_err    = numpy.zeros(self.nbins) 
    772856         
    773857        # This If finds qmax within ROI defined by sector lines 
     
    775859        temp0=1000000 
    776860        temp1=0 
    777         for i in range(numpy.size(data,1)):                  
     861        for i in range(numpy.size(data,1)):   
     862            dx = pixel_width_x*(i+0.5 - center_x)                   
    778863            for j in range(numpy.size(data,0)): 
    779                 #number of sub-bin: default = 2 (ie., 1 bin: no sub-bin) 
    780                 nsubbins = 1 
    781  
    782                 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 
    783                 for x_subbins in range(nsubbins): 
    784                     for y_subbins in range(nsubbins): 
     864                     
     865                dy = pixel_width_y*(j+0.5 - center_y) 
     866                q_value = get_q(dx, dy, det_dist, wavelength) 
     867                # Compute phi and check whether it's within the limits 
     868                phi_value=math.atan2(dy,dx)+math.pi 
     869                if self.phi_max>2*math.pi: 
     870                    self.phi_max=self.phi_max-2*math.pi 
     871                if self.phi_min<0: 
     872                    self.phi_max=self.phi_max+2*math.pi 
     873                 
     874                #In case of two ROI (symmetric major and minor wings)(for 'q2') 
     875                if run.lower()=='q2': 
     876                    if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 
     877                        temp_max=self.phi_max+math.pi 
     878                        temp_min=self.phi_min+math.pi 
     879                    else: 
     880                        temp_max=self.phi_max 
     881                        temp_min=self.phi_min 
     882                        
     883                    if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 
     884                        if (phi_value<temp_min  or phi_value>temp_max): 
     885                             if (phi_value<temp_min-math.pi  or phi_value>temp_max-math.pi): 
     886                                 continue 
     887                    if (self.phi_max<self.phi_min): 
     888                        tmp_max=self.phi_max+math.pi 
     889                        tmp_min=self.phi_min-math.pi 
     890                    else: 
     891                        tmp_max=self.phi_max 
     892                        tmp_min=self.phi_min 
     893                    if (tmp_min<math.pi and tmp_max>math.pi): 
     894                        if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 
     895                            continue 
     896                #In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'. 
     897                elif run.lower()=='q':  
     898                    if (self.phi_max>=self.phi_min): 
     899                        if (phi_value<self.phi_min  or phi_value>self.phi_max): 
     900                            continue 
     901                    else: 
     902                        if (phi_value<self.phi_min and phi_value>self.phi_max): 
     903                            continue    
     904                if q_value<qmin or q_value>qmax: 
     905                        continue                     
    785906                         
    786                         dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 
    787                         dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 
    788                          
    789                         q_value = get_q(dx, dy, det_dist, wavelength) 
    790  
    791                         # Compute phi and check whether it's within the limits 
    792                         phi_value=math.atan2(dy,dx)+math.pi 
    793                         if self.phi_max>2*math.pi: 
    794                             self.phi_max=self.phi_max-2*math.pi 
    795                         if self.phi_min<0: 
    796                             self.phi_max=self.phi_max+2*math.pi 
    797                          
    798                         #In case of two ROI (symmetric major and minor wings)(for 'q2') 
    799                         if run.lower()=='q2': 
    800                             if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 
    801                                 temp_max=self.phi_max+math.pi 
    802                                 temp_min=self.phi_min+math.pi 
    803                             else: 
    804                                 temp_max=self.phi_max 
    805                                 temp_min=self.phi_min 
    806                                 
    807                             if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 
    808                                 if (phi_value<temp_min  or phi_value>temp_max): 
    809                                      if (phi_value<temp_min-math.pi  or phi_value>temp_max-math.pi): 
    810                                          continue 
    811                             if (self.phi_max<self.phi_min): 
    812                                 tmp_max=self.phi_max+math.pi 
    813                                 tmp_min=self.phi_min-math.pi 
    814                             else: 
    815                                 tmp_max=self.phi_max 
    816                                 tmp_min=self.phi_min 
    817                             if (tmp_min<math.pi and tmp_max>math.pi): 
    818                                 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 
    819                                     continue 
    820                         #In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'. 
    821                         elif run.lower()=='q':  
    822                             if (self.phi_max>=self.phi_min): 
    823                                 if (phi_value<self.phi_min  or phi_value>self.phi_max): 
    824                                     continue 
    825                             else: 
    826                                 if (phi_value<self.phi_min and phi_value>self.phi_max): 
    827                                     continue    
    828                         if q_value<qmin or q_value>qmax: 
    829                                 continue                     
    830                                  
    831                         if run.lower()=='phi': 
    832                             if temp1<phi_value: 
    833                                 temp1=phi_value 
    834                             if temp0>phi_value: 
    835                                 temp0=phi_value    
    836                                                                                                   
    837                         elif temp<q_value: 
    838                             temp=q_value 
     907                if run.lower()=='phi': 
     908                    if temp1<phi_value: 
     909                        temp1=phi_value 
     910                    if temp0>phi_value: 
     911                        temp0=phi_value    
     912                                                                                          
     913                elif temp<q_value: 
     914                    temp=q_value 
    839915                     
    840916        if run.lower()=='phi': 
     
    845921        #Beam center is already corrected, but the calculation below assumed it was not. 
    846922        # Thus Beam center shifted back to uncorrected value. ToDo: cleanup the mess. 
    847         #center_x=center_x+0.5  
    848         #center_y=center_y+0.5          
    849         for i in range(numpy.size(data,1)):           
     923        center_x=center_x+0.5  
     924        center_y=center_y+0.5          
     925        for i in range(numpy.size(data,1)): 
     926            dx = pixel_width_x*(i+0.5 - center_x) 
     927             
     928            # Min and max x-value for the pixel 
     929            minx = pixel_width_x*(i - center_x) 
     930            maxx = pixel_width_x*(i+1.0 - center_x) 
     931             
    850932            for j in range(numpy.size(data,0)): 
    851                 #number of sub-bin: default = 1 (ie., 1 bin: no sub-bin) 
    852                 nsubbins = 1 
    853  
    854                 #if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \ 
    855                 #    or is_intercept(self.r_min, q_00, q_01, q_10, q_11): 
    856                 #    # Number of sub-bins near the boundary of ROI in x-direction 
    857                 #    nsubbins = 3 
    858                  
    859                 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 
    860                 for x_subbins in range(nsubbins): 
    861                     for y_subbins in range(nsubbins): 
    862                          
    863                         dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 
    864                         dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 
    865                          
    866                         q_value = get_q(dx, dy, det_dist, wavelength) 
    867                         # Min and max x-value for the subpixel 
    868                         minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x) 
    869                         maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x) 
    870                         # Min and max y-value for the subpixel 
    871                         miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y) 
    872                         maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y) 
    873              
    874                  
    875                         # Calculate the q-value for each corner 
    876                         # q_[x min or max][y min or max] 
    877                         q_00 = get_q(minx, miny, det_dist, wavelength) 
    878                         q_01 = get_q(minx, maxy, det_dist, wavelength) 
    879                         q_10 = get_q(maxx, miny, det_dist, wavelength) 
    880                         q_11 = get_q(maxx, maxy, det_dist, wavelength) 
    881                          
    882                         # Compute phi and check whether it's within the limits 
    883                         phi_value=math.atan2(dy,dx)+math.pi 
    884                         if self.phi_max>2*math.pi: 
    885                             self.phi_max=self.phi_max-2*math.pi 
    886                         if self.phi_min<0: 
    887                             self.phi_max=self.phi_max+2*math.pi 
     933                dy = pixel_width_y*(j+0.5 - center_y) 
     934             
     935                q_value = get_q(dx, dy, det_dist, wavelength) 
     936 
     937                # Min and max y-value for the pixel 
     938                miny = pixel_width_y*(j - center_y) 
     939                maxy = pixel_width_y*(j+1.0 - center_y) 
     940                 
     941                # Calculate the q-value for each corner 
     942                # q_[x min or max][y min or max] 
     943                q_00 = get_q(minx, miny, det_dist, wavelength) 
     944                q_01 = get_q(minx, maxy, det_dist, wavelength) 
     945                q_10 = get_q(maxx, miny, det_dist, wavelength) 
     946                q_11 = get_q(maxx, maxy, det_dist, wavelength) 
     947                 
     948                # Compute phi and check whether it's within the limits 
     949                phi_value=math.atan2(dy,dx)+math.pi 
     950                if self.phi_max>2*math.pi: 
     951                    self.phi_max=self.phi_max-2*math.pi 
     952                if self.phi_min<0: 
     953                    self.phi_max=self.phi_max+2*math.pi 
     954                     
     955                # Look for intercept between each side of the pixel 
     956                # and the constant-q ring for qmax 
     957                frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 
     958                 
     959                # Look for intercept between each side of the pixel 
     960                # and the constant-q ring for qmin 
     961                frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 
     962                 
     963                # We are interested in the region between qmin and qmax 
     964                # therefore the fraction of the surface of the pixel 
     965                # that we will use to calculate the number of counts to  
     966                # include is given by: 
     967                 
     968                frac = frac_max - frac_min 
     969 
     970                #In case of two ROI (symmetric major and minor regions)(for 'q2') 
     971                if run.lower()=='q2': 
     972                    if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 
     973                        temp_max=self.phi_max+math.pi 
     974                        temp_min=self.phi_min+math.pi 
     975                    else: 
     976                        temp_max=self.phi_max 
     977                        temp_min=self.phi_min 
     978                        
     979                    if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 
     980                        if (phi_value<temp_min  or phi_value>temp_max): 
     981                            if (phi_value<temp_min-math.pi  or phi_value>temp_max-math.pi): 
     982                                continue 
     983                    if (self.phi_max<self.phi_min): 
     984                        tmp_max=self.phi_max+math.pi 
     985                        tmp_min=self.phi_min-math.pi 
     986                    else: 
     987                        tmp_max=self.phi_max 
     988                        tmp_min=self.phi_min 
     989                    if (tmp_min<math.pi and tmp_max>math.pi): 
     990                        if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 
     991                            continue 
     992                #In case of one ROI (major only)(i.e.,for 'q' and 'phi') 
     993                else:  
     994                    if (self.phi_max>=self.phi_min): 
     995                        if (phi_value<self.phi_min  or phi_value>self.phi_max): 
     996                            continue 
    888997                             
    889                         # Look for intercept between each side of the pixel 
    890                         # and the constant-q ring for qmax 
    891                         frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 
    892                          
    893                         # Look for intercept between each side of the pixel 
    894                         # and the constant-q ring for qmin 
    895                         frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 
    896                          
    897                         # We are interested in the region between qmin and qmax 
    898                         # therefore the fraction of the surface of the pixel 
    899                         # that we will use to calculate the number of counts to  
    900                         # include is given by: 
    901                          
    902                         frac = frac_max - frac_min 
    903          
    904                         #In case of two ROI (symmetric major and minor regions)(for 'q2') 
    905                         if run.lower()=='q2': 
    906                             if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 
    907                                 temp_max=self.phi_max+math.pi 
    908                                 temp_min=self.phi_min+math.pi 
    909                             else: 
    910                                 temp_max=self.phi_max 
    911                                 temp_min=self.phi_min 
    912                                 
    913                             if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 
    914                                 if (phi_value<temp_min  or phi_value>temp_max): 
    915                                     if (phi_value<temp_min-math.pi  or phi_value>temp_max-math.pi): 
    916                                         continue 
    917                             if (self.phi_max<self.phi_min): 
    918                                 tmp_max=self.phi_max+math.pi 
    919                                 tmp_min=self.phi_min-math.pi 
    920                             else: 
    921                                 tmp_max=self.phi_max 
    922                                 tmp_min=self.phi_min 
    923                             if (tmp_min<math.pi and tmp_max>math.pi): 
    924                                 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 
    925                                     continue 
    926                         #In case of one ROI (major only)(i.e.,for 'q' and 'phi') 
    927                         else:  
    928                             if (self.phi_max>=self.phi_min): 
    929                                 if (phi_value<self.phi_min  or phi_value>self.phi_max): 
    930                                     continue 
    931                                                                
    932                                 # Check which type of averaging we need 
    933                                 if run.lower()=='phi':  
    934                                     i_bin = int(math.floor(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min)))  
    935                                 else: 
    936                                     # If we don't need this pixel, skip the rest of the work 
    937                                     #TODO: an improvement here would be to compute the average 
    938                                     # Q for the pixel from the part that is covered by 
    939                                     # the ring defined by q_min/q_max rather than the complete 
    940                                     # pixel  
    941                                     i_bin = int(math.floor(self.nbins*(q_value-qmin)/(qmax-qmin))) 
    942                                      
    943                             else: 
    944                                 if (phi_value<self.phi_min and phi_value>self.phi_max): 
    945                                     continue 
    946                         #print "qmax=",qmax,qmin        
    947          
    948                         if q_value<qmin or q_value>qmax: 
     998                    else: 
     999                        if (phi_value<self.phi_min and phi_value>self.phi_max): 
    9491000                            continue 
    950                                                              
    951                         # Check which type of averaging we need 
    952                         if run.lower()=='phi':  
    953                             i_bin = int(math.floor(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) 
    954                         else: 
    955                             # If we don't need this pixel, skip the rest of the work 
    956                             #TODO: an improvement here would be to compute the average 
    957                             # Q for the pixel from the part that is covered by 
    958                             # the ring defined by q_min/q_max rather than the complete 
    959                             # pixel  
    960                             i_bin = int(math.floor(self.nbins*(q_value-qmin)/(qmax-qmin))) 
    961                              
    962                         #Make sure that no out of bound happens 
    963                         #This is needed because of a numerical error in math.ceil(). 
    964                         if i_bin >= self.nbins: 
    965                             i_bin = self.nbins - 1 
    966                         elif i_bin < 0: 
    967                             i_bin = 0 
    968                                    
    969                         try: 
    970                             y[i_bin] += frac * data[j][i] 
    971                         except: 
    972                             import sys 
    973                             print sys.exc_value 
    974                             print i_bin, frac 
    975                          
    976                         if data2D.err_data == None or data2D.err_data[j][i]==0.0: 
    977                             y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 
    978                         else: 
    979                             y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 
    980                         y_counts[i_bin] += frac 
    981          
    982  
     1001                #print "qmax=",qmax,qmin        
     1002 
     1003                if q_value<qmin or q_value>qmax: 
     1004                    continue 
     1005                                                     
     1006                # Check which type of averaging we need 
     1007                if run.lower()=='phi':  
     1008                    i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 
     1009                else: 
     1010                    # If we don't need this pixel, skip the rest of the work 
     1011                    #TODO: an improvement here would be to compute the average 
     1012                    # Q for the pixel from the part that is covered by 
     1013                    # the ring defined by q_min/q_max rather than the complete 
     1014                    # pixel  
     1015                    i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1 
     1016                            
     1017                try: 
     1018                    y[i_bin] += frac * data[j][i] 
     1019                except: 
     1020                    import sys 
     1021                    print sys.exc_value 
     1022                    print i_bin, frac 
     1023                 
     1024                if data2D.err_data == None or data2D.err_data[j][i]==0.0: 
     1025                    y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 
     1026                else: 
     1027                    y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 
     1028                y_counts[i_bin] += frac 
     1029         
    9831030        for i in range(self.nbins): 
    9841031            y[i] = y[i] / y_counts[i] 
     
    10111058        """ 
    10121059        return self._agv(data2D, 'phi') 
     1060 
     1061class SectorQold(_Sector): 
     1062    """ 
     1063        Sector average as a function of Q. 
     1064        I(Q) is return and the data is averaged over phi. 
     1065         
     1066        A sector is defined by r_min, r_max, phi_min, phi_max. 
     1067        The number of bin in Q also has to be defined. 
     1068    """ 
     1069    def __call__(self, data2D): 
     1070        """ 
     1071            Perform sector average and return I(Q). 
     1072             
     1073            @param data2D: Data2D object 
     1074            @return: Data1D object 
     1075        """ 
     1076        return self._agv(data2D, 'q') 
    10131077     
    10141078class SectorQ(_Sector): 
Note: See TracChangeset for help on using the changeset viewer.