Changeset ccc7192 in sasview


Ignore:
Timestamp:
Apr 6, 2017 11:17:21 AM (8 years ago)
Author:
Ricardo Ferraz Leal <ricleal@…>
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.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
fa4af76
Parents:
9dda8cc
Message:

Cleaning and Refactoring

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/sas/sascalc/dataloader/manipulations.py

    r9dda8cc rccc7192  
    7272    return phi_out 
    7373 
    74  
    75 def reader2D_converter(data2d=None): 
    76     """ 
    77     convert old 2d format opened by IhorReader or danse_reader 
    78     to new Data2D format 
    79  
    80     :param data2d: 2d array of Data2D object 
    81     :return: 1d arrays of Data2D object 
    82  
    83     """ 
    84     if data2d.data is None or data2d.x_bins is None or data2d.y_bins is None: 
    85         raise ValueError("Can't convert this data: data=None...") 
    86     new_x = np.tile(data2d.x_bins, (len(data2d.y_bins), 1)) 
    87     new_y = np.tile(data2d.y_bins, (len(data2d.x_bins), 1)) 
    88     new_y = new_y.swapaxes(0, 1) 
    89  
    90     new_data = data2d.data.flatten() 
    91     qx_data = new_x.flatten() 
    92     qy_data = new_y.flatten() 
    93     q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) 
    94     if data2d.err_data is None or np.any(data2d.err_data <= 0): 
    95         new_err_data = np.sqrt(np.abs(new_data)) 
    96     else: 
    97         new_err_data = data2d.err_data.flatten() 
    98     mask = np.ones(len(new_data), dtype=bool) 
    99  
    100     # TODO: make sense of the following two lines... 
    101     #from sas.sascalc.dataloader.data_info import Data2D 
    102     #output = Data2D() 
    103     output = data2d 
    104     output.data = new_data 
    105     output.err_data = new_err_data 
    106     output.qx_data = qx_data 
    107     output.qy_data = qy_data 
    108     output.q_data = q_data 
    109     output.mask = mask 
    110  
    111     return output 
    112  
    113  
    114 class _Slab(object): 
    115     """ 
    116     Compute average I(Q) for a region of interest 
    117     """ 
    118  
    119     def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, 
    120                  y_max=0.0, bin_width=0.001): 
    121         # Minimum Qx value [A-1] 
    122         self.x_min = x_min 
    123         # Maximum Qx value [A-1] 
    124         self.x_max = x_max 
    125         # Minimum Qy value [A-1] 
    126         self.y_min = y_min 
    127         # Maximum Qy value [A-1] 
    128         self.y_max = y_max 
    129         # Bin width (step size) [A-1] 
    130         self.bin_width = bin_width 
    131         # If True, I(|Q|) will be return, otherwise, 
    132         # negative q-values are allowed 
    133         self.fold = False 
    134  
    135     def __call__(self, data2D): 
    136         return NotImplemented 
    137  
    138     def _avg(self, data2D, maj): 
    139         """ 
    140         Compute average I(Q_maj) for a region of interest. 
    141         The major axis is defined as the axis of Q_maj. 
    142         The minor axis is the axis that we average over. 
    143  
    144         :param data2D: Data2D object 
    145         :param maj_min: min value on the major axis 
    146         :return: Data1D object 
    147         """ 
    148         if len(data2D.detector) > 1: 
    149             msg = "_Slab._avg: invalid number of " 
    150             msg += " detectors: %g" % len(data2D.detector) 
    151             raise RuntimeError(msg) 
    152  
    153         # Get data 
    154         data = data2D.data[np.isfinite(data2D.data)] 
    155         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    156         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    157         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    158  
    159         # Build array of Q intervals 
    160         if maj == 'x': 
    161             if self.fold: 
    162                 x_min = 0 
    163             else: 
    164                 x_min = self.x_min 
    165             nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) 
    166         elif maj == 'y': 
    167             if self.fold: 
    168                 y_min = 0 
    169             else: 
    170                 y_min = self.y_min 
    171             nbins = int(math.ceil((self.y_max - y_min) / self.bin_width)) 
    172         else: 
    173             raise RuntimeError("_Slab._avg: unrecognized axis %s" % str(maj)) 
    174  
    175         x = np.zeros(nbins) 
    176         y = np.zeros(nbins) 
    177         err_y = np.zeros(nbins) 
    178         y_counts = np.zeros(nbins) 
    179  
    180         # Average pixelsize in q space 
    181         for npts in range(len(data)): 
    182             # default frac 
    183             frac_x = 0 
    184             frac_y = 0 
    185             # get ROI 
    186             if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]: 
    187                 frac_x = 1 
    188             if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]: 
    189                 frac_y = 1 
    190             frac = frac_x * frac_y 
    191  
    192             if frac == 0: 
    193                 continue 
    194             # binning: find axis of q 
    195             if maj == 'x': 
    196                 q_value = qx_data[npts] 
    197                 min_value = x_min 
    198             if maj == 'y': 
    199                 q_value = qy_data[npts] 
    200                 min_value = y_min 
    201             if self.fold and q_value < 0: 
    202                 q_value = -q_value 
    203             # bin 
    204             i_q = int(math.ceil((q_value - min_value) / self.bin_width)) - 1 
    205  
    206             # skip outside of max bins 
    207             if i_q < 0 or i_q >= nbins: 
    208                 continue 
    209  
    210             # TODO: find better definition of x[i_q] based on q_data 
    211             # min_value + (i_q + 1) * self.bin_width / 2.0 
    212             x[i_q] += frac * q_value 
    213             y[i_q] += frac * data[npts] 
    214  
    215             if err_data is None or err_data[npts] == 0.0: 
    216                 if data[npts] < 0: 
    217                     data[npts] = -data[npts] 
    218                 err_y[i_q] += frac * frac * data[npts] 
    219             else: 
    220                 err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] 
    221             y_counts[i_q] += frac 
    222  
    223         # Average the sums 
    224         for n in range(nbins): 
    225             err_y[n] = math.sqrt(err_y[n]) 
    226  
    227         err_y = err_y / y_counts 
    228         y = y / y_counts 
    229         x = x / y_counts 
    230         idx = (np.isfinite(y) & np.isfinite(x)) 
    231  
    232         if not idx.any(): 
    233             msg = "Average Error: No points inside ROI to average..." 
    234             raise ValueError(msg) 
    235         return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) 
    236  
    237  
    238 class SlabY(_Slab): 
    239     """ 
    240     Compute average I(Qy) for a region of interest 
    241     """ 
    242  
    243     def __call__(self, data2D): 
    244         """ 
    245         Compute average I(Qy) for a region of interest 
    246  
    247         :param data2D: Data2D object 
    248         :return: Data1D object 
    249         """ 
    250         return self._avg(data2D, 'y') 
    251  
    252  
    253 class SlabX(_Slab): 
    254     """ 
    255     Compute average I(Qx) for a region of interest 
    256     """ 
    257  
    258     def __call__(self, data2D): 
    259         """ 
    260         Compute average I(Qx) for a region of interest 
    261         :param data2D: Data2D object 
    262         :return: Data1D object 
    263         """ 
    264         return self._avg(data2D, 'x') 
    265  
    266  
    267 class Boxsum(object): 
    268     """ 
    269     Perform the sum of counts in a 2D region of interest. 
    270     """ 
    271  
    272     def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
    273         # Minimum Qx value [A-1] 
    274         self.x_min = x_min 
    275         # Maximum Qx value [A-1] 
    276         self.x_max = x_max 
    277         # Minimum Qy value [A-1] 
    278         self.y_min = y_min 
    279         # Maximum Qy value [A-1] 
    280         self.y_max = y_max 
    281  
    282     def __call__(self, data2D): 
    283         """ 
    284         Perform the sum in the region of interest 
    285  
    286         :param data2D: Data2D object 
    287         :return: number of counts, error on number of counts, 
    288             number of points summed 
    289         """ 
    290         y, err_y, y_counts = self._sum(data2D) 
    291  
    292         # Average the sums 
    293         counts = 0 if y_counts == 0 else y 
    294         error = 0 if y_counts == 0 else math.sqrt(err_y) 
    295  
    296         # Added y_counts to return, SMK & PDB, 04/03/2013 
    297         return counts, error, y_counts 
    298  
    299     def _sum(self, data2D): 
    300         """ 
    301         Perform the sum in the region of interest 
    302  
    303         :param data2D: Data2D object 
    304         :return: number of counts, 
    305             error on number of counts, number of entries summed 
    306         """ 
    307         if len(data2D.detector) > 1: 
    308             msg = "Circular averaging: invalid number " 
    309             msg += "of detectors: %g" % len(data2D.detector) 
    310             raise RuntimeError(msg) 
    311         # Get data 
    312         data = data2D.data[np.isfinite(data2D.data)] 
    313         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    314         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    315         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    316  
    317         y = 0.0 
    318         err_y = 0.0 
    319         y_counts = 0.0 
    320  
    321         # Average pixelsize in q space 
    322         for npts in range(len(data)): 
    323             # default frac 
    324             frac_x = 0 
    325             frac_y = 0 
    326  
    327             # get min and max at each points 
    328             qx = qx_data[npts] 
    329             qy = qy_data[npts] 
    330  
    331             # get the ROI 
    332             if self.x_min <= qx and self.x_max > qx: 
    333                 frac_x = 1 
    334             if self.y_min <= qy and self.y_max > qy: 
    335                 frac_y = 1 
    336             # Find the fraction along each directions 
    337             frac = frac_x * frac_y 
    338             if frac == 0: 
    339                 continue 
    340             y += frac * data[npts] 
    341             if err_data is None or err_data[npts] == 0.0: 
    342                 if data[npts] < 0: 
    343                     data[npts] = -data[npts] 
    344                 err_y += frac * frac * data[npts] 
    345             else: 
    346                 err_y += frac * frac * err_data[npts] * err_data[npts] 
    347             y_counts += frac 
    348         return y, err_y, y_counts 
    349  
    350  
    351 class Boxavg(Boxsum): 
    352     """ 
    353     Perform the average of counts in a 2D region of interest. 
    354     """ 
    355  
    356     def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
    357         super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, 
    358                                      y_min=y_min, y_max=y_max) 
    359  
    360     def __call__(self, data2D): 
    361         """ 
    362         Perform the sum in the region of interest 
    363  
    364         :param data2D: Data2D object 
    365         :return: average counts, error on average counts 
    366  
    367         """ 
    368         y, err_y, y_counts = self._sum(data2D) 
    369  
    370         # Average the sums 
    371         counts = 0 if y_counts == 0 else y / y_counts 
    372         error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts 
    373  
    374         return counts, error 
    375  
    376  
    37774def get_pixel_fraction_square(x, xmin, xmax): 
    37875    """ 
    … …  
    39794        return 1.0 
    39895 
    399  
    400 class CircularAverage(object): 
    401     """ 
    402     Perform circular averaging on 2D data 
    403  
    404     The data returned is the distribution of counts 
    405     as a function of Q 
    406     """ 
    407  
    408     def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 
    409         # Minimum radius included in the average [A-1] 
    410         self.r_min = r_min 
    411         # Maximum radius included in the average [A-1] 
    412         self.r_max = r_max 
    413         # Bin width (step size) [A-1] 
    414         self.bin_width = bin_width 
    415  
    416     def __call__(self, data2D, ismask=False): 
    417         """ 
    418         Perform circular averaging on the data 
    419  
    420         :param data2D: Data2D object 
    421         :return: Data1D object 
    422         """ 
    423         # Get data W/ finite values 
    424         data = data2D.data[np.isfinite(data2D.data)] 
    425         q_data = data2D.q_data[np.isfinite(data2D.data)] 
    426         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    427         mask_data = data2D.mask[np.isfinite(data2D.data)] 
    428  
    429         dq_data = None 
    430  
    431         # Get the dq for resolution averaging 
    432         if data2D.dqx_data != None and data2D.dqy_data != None: 
    433             # The pinholes and det. pix contribution present 
    434             # in both direction of the 2D which must be subtracted when 
    435             # converting to 1D: dq_overlap should calculated ideally at 
    436             # q = 0. Note This method works on only pinhole geometry. 
    437             # Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. 
    438             z_max = max(data2D.q_data) 
    439             z_min = min(data2D.q_data) 
    440             x_max = data2D.dqx_data[data2D.q_data[z_max]] 
    441             x_min = data2D.dqx_data[data2D.q_data[z_min]] 
    442             y_max = data2D.dqy_data[data2D.q_data[z_max]] 
    443             y_min = data2D.dqy_data[data2D.q_data[z_min]] 
    444             # Find qdx at q = 0 
    445             dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
    446             # when extrapolation goes wrong 
    447             if dq_overlap_x > min(data2D.dqx_data): 
    448                 dq_overlap_x = min(data2D.dqx_data) 
    449             dq_overlap_x *= dq_overlap_x 
    450             # Find qdx at q = 0 
    451             dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
    452             # when extrapolation goes wrong 
    453             if dq_overlap_y > min(data2D.dqy_data): 
    454                 dq_overlap_y = min(data2D.dqy_data) 
    455             # get dq at q=0. 
    456             dq_overlap_y *= dq_overlap_y 
    457  
    458             dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 
    459             # Final protection of dq 
    460             if dq_overlap < 0: 
    461                 dq_overlap = y_min 
    462             dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] 
    463             dqy_data = data2D.dqy_data[np.isfinite( 
    464                 data2D.data)] - dq_overlap 
    465             # def; dqx_data = dq_r dqy_data = dq_phi 
    466             # Convert dq 2D to 1D here 
    467             dqx = dqx_data * dqx_data 
    468             dqy = dqy_data * dqy_data 
    469             dq_data = np.add(dqx, dqy) 
    470             dq_data = np.sqrt(dq_data) 
    471  
    472         #q_data_max = np.max(q_data) 
    473         if len(data2D.q_data) == None: 
    474             msg = "Circular averaging: invalid q_data: %g" % data2D.q_data 
    475             raise RuntimeError(msg) 
    476  
    477         # Build array of Q intervals 
    478         nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) 
    479  
    480         x = np.zeros(nbins) 
    481         y = np.zeros(nbins) 
    482         err_y = np.zeros(nbins) 
    483         err_x = np.zeros(nbins) 
    484         y_counts = np.zeros(nbins) 
    485  
    486         for npt in range(len(data)): 
    487  
    488             if ismask and not mask_data[npt]: 
    489                 continue 
    490  
    491             frac = 0 
    492  
    493             # q-value at the pixel (j,i) 
    494             q_value = q_data[npt] 
    495             data_n = data[npt] 
    496  
    497             # No need to calculate the frac when all data are within range 
    498             if self.r_min >= self.r_max: 
    499                 raise ValueError("Limit Error: min > max") 
    500  
    501             if self.r_min <= q_value and q_value <= self.r_max: 
    502                 frac = 1 
    503             if frac == 0: 
    504                 continue 
    505             i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 
    506  
    507             # Take care of the edge case at phi = 2pi. 
    508             if i_q == nbins: 
    509                 i_q = nbins - 1 
    510             y[i_q] += frac * data_n 
    511             # Take dqs from data to get the q_average 
    512             x[i_q] += frac * q_value 
    513             if err_data is None or err_data[npt] == 0.0: 
    514                 if data_n < 0: 
    515                     data_n = -data_n 
    516                 err_y[i_q] += frac * frac * data_n 
    517             else: 
    518                 err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] 
    519             if dq_data != None: 
    520                 # To be consistent with dq calculation in 1d reduction, 
    521                 # we need just the averages (not quadratures) because 
    522                 # it should not depend on the number of the q points 
    523                 # in the qr bins. 
    524                 err_x[i_q] += frac * dq_data[npt] 
    525             else: 
    526                 err_x = None 
    527             y_counts[i_q] += frac 
    528  
    529         # Average the sums 
    530         for n in range(nbins): 
    531             if err_y[n] < 0: 
    532                 err_y[n] = -err_y[n] 
    533             err_y[n] = math.sqrt(err_y[n]) 
    534             # if err_x != None: 
    535             #    err_x[n] = math.sqrt(err_x[n]) 
    536  
    537         err_y = err_y / y_counts 
    538         err_y[err_y == 0] = np.average(err_y) 
    539         y = y / y_counts 
    540         x = x / y_counts 
    541         idx = (np.isfinite(y)) & (np.isfinite(x)) 
    542  
    543         if err_x != None: 
    544             d_x = err_x[idx] / y_counts[idx] 
    545         else: 
    546             d_x = None 
    547  
    548         if not idx.any(): 
    549             msg = "Average Error: No points inside ROI to average..." 
    550             raise ValueError(msg) 
    551  
    552         return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) 
    553  
    554  
    555 class Ring(object): 
    556     """ 
    557     Defines a ring on a 2D data set. 
    558     The ring is defined by r_min, r_max, and 
    559     the position of the center of the ring. 
    560  
    561     The data returned is the distribution of counts 
    562     around the ring as a function of phi. 
    563  
    564     Phi_min and phi_max should be defined between 0 and 2*pi 
    565     in anti-clockwise starting from the x- axis on the left-hand side 
    566     """ 
    567     # Todo: remove center. 
    568  
    569     def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=36): 
    570         # Minimum radius 
    571         self.r_min = r_min 
    572         # Maximum radius 
    573         self.r_max = r_max 
    574         # Center of the ring in x 
    575         self.center_x = center_x 
    576         # Center of the ring in y 
    577         self.center_y = center_y 
    578         # Number of angular bins 
    579         self.nbins_phi = nbins 
    580  
    581     def __call__(self, data2D): 
    582         """ 
    583         Apply the ring to the data set. 
    584         Returns the angular distribution for a given q range 
    585  
    586         :param data2D: Data2D object 
    587  
    588         :return: Data1D object 
    589         """ 
    590         if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    591             raise RuntimeError("Ring averaging only take plottable_2D objects") 
    592  
    593         Pi = math.pi 
    594  
    595         # Get data 
    596         data = data2D.data[np.isfinite(data2D.data)] 
    597         q_data = data2D.q_data[np.isfinite(data2D.data)] 
    598         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    599         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    600         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    601  
    602         # Set space for 1d outputs 
    603         phi_bins = np.zeros(self.nbins_phi) 
    604         phi_counts = np.zeros(self.nbins_phi) 
    605         phi_values = np.zeros(self.nbins_phi) 
    606         phi_err = np.zeros(self.nbins_phi) 
    607  
    608         # Shift to apply to calculated phi values in order 
    609         # to center first bin at zero 
    610         phi_shift = Pi / self.nbins_phi 
    611  
    612         for npt in range(len(data)): 
    613             frac = 0 
    614             # q-value at the point (npt) 
    615             q_value = q_data[npt] 
    616             data_n = data[npt] 
    617  
    618             # phi-value at the point (npt) 
    619             phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi 
    620  
    621             if self.r_min <= q_value and q_value <= self.r_max: 
    622                 frac = 1 
    623             if frac == 0: 
    624                 continue 
    625             # binning 
    626             i_phi = int(math.floor((self.nbins_phi) * 
    627                                    (phi_value + phi_shift) / (2 * Pi))) 
    628  
    629             # Take care of the edge case at phi = 2pi. 
    630             if i_phi >= self.nbins_phi: 
    631                 i_phi = 0 
    632             phi_bins[i_phi] += frac * data[npt] 
    633  
    634             if err_data is None or err_data[npt] == 0.0: 
    635                 if data_n < 0: 
    636                     data_n = -data_n 
    637                 phi_err[i_phi] += frac * frac * math.fabs(data_n) 
    638             else: 
    639                 phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt] 
    640             phi_counts[i_phi] += frac 
    641  
    642         for i in range(self.nbins_phi): 
    643             phi_bins[i] = phi_bins[i] / phi_counts[i] 
    644             phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i] 
    645             phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) 
    646  
    647         idx = (np.isfinite(phi_bins)) 
    648  
    649         if not idx.any(): 
    650             msg = "Average Error: No points inside ROI to average..." 
    651             raise ValueError(msg) 
    652         # elif len(phi_bins[idx])!= self.nbins_phi: 
    653         #    print "resulted",self.nbins_phi- len(phi_bins[idx]) 
    654         #,"empty bin(s) due to tight binning..." 
    655         return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) 
    656  
     96def get_intercept(q, q_0, q_1): 
     97    """ 
     98    Returns the fraction of the side at which the 
     99    q-value intercept the pixel, None otherwise. 
     100    The values returned is the fraction ON THE SIDE 
     101    OF THE LOWEST Q. :: 
     102 
     103            A           B 
     104        +-----------+--------+    <--- pixel size 
     105        0                    1 
     106        Q_0 -------- Q ----- Q_1   <--- equivalent Q range 
     107        if Q_1 > Q_0, A is returned 
     108        if Q_1 < Q_0, B is returned 
     109        if Q is outside the range of [Q_0, Q_1], None is returned 
     110 
     111    """ 
     112    if q_1 > q_0: 
     113        if q > q_0 and q <= q_1: 
     114            return (q - q_0) / (q_1 - q_0) 
     115    else: 
     116        if q > q_1 and q <= q_0: 
     117            return (q - q_1) / (q_0 - q_1) 
     118    return None 
    657119 
    658120def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11): 
    … …  
    719181    return frac_max 
    720182 
    721  
    722 def get_intercept(q, q_0, q_1): 
    723     """ 
    724     Returns the fraction of the side at which the 
    725     q-value intercept the pixel, None otherwise. 
    726     The values returned is the fraction ON THE SIDE 
    727     OF THE LOWEST Q. :: 
    728  
    729             A           B 
    730         +-----------+--------+    <--- pixel size 
    731         0                    1 
    732         Q_0 -------- Q ----- Q_1   <--- equivalent Q range 
    733         if Q_1 > Q_0, A is returned 
    734         if Q_1 < Q_0, B is returned 
    735         if Q is outside the range of [Q_0, Q_1], None is returned 
    736  
    737     """ 
    738     if q_1 > q_0: 
    739         if q > q_0 and q <= q_1: 
    740             return (q - q_0) / (q_1 - q_0) 
     183def get_dq_data(data2D): 
     184    ''' 
     185    Get the dq for resolution averaging 
     186    The pinholes and det. pix contribution present 
     187    in both direction of the 2D which must be subtracted when 
     188    converting to 1D: dq_overlap should calculated ideally at 
     189    q = 0. Note This method works on only pinhole geometry. 
     190    Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. 
     191    ''' 
     192    z_max = max(data2D.q_data) 
     193    z_min = min(data2D.q_data) 
     194    x_max = data2D.dqx_data[data2D.q_data[z_max]] 
     195    x_min = data2D.dqx_data[data2D.q_data[z_min]] 
     196    y_max = data2D.dqy_data[data2D.q_data[z_max]] 
     197    y_min = data2D.dqy_data[data2D.q_data[z_min]] 
     198    # Find qdx at q = 0 
     199    dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
     200    # when extrapolation goes wrong 
     201    if dq_overlap_x > min(data2D.dqx_data): 
     202        dq_overlap_x = min(data2D.dqx_data) 
     203    dq_overlap_x *= dq_overlap_x 
     204    # Find qdx at q = 0 
     205    dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
     206    # when extrapolation goes wrong 
     207    if dq_overlap_y > min(data2D.dqy_data): 
     208        dq_overlap_y = min(data2D.dqy_data) 
     209    # get dq at q=0. 
     210    dq_overlap_y *= dq_overlap_y 
     211 
     212    dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 
     213    # Final protection of dq 
     214    if dq_overlap < 0: 
     215        dq_overlap = y_min 
     216    dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] 
     217    dqy_data = data2D.dqy_data[np.isfinite( 
     218        data2D.data)] - dq_overlap 
     219    # def; dqx_data = dq_r dqy_data = dq_phi 
     220    # Convert dq 2D to 1D here 
     221    dq_data = np.sqrt(dqx_data**2 + dqx_data**2)  
     222    return dq_data 
     223 
     224################################################################################ 
     225 
     226def reader2D_converter(data2d=None): 
     227    """ 
     228    convert old 2d format opened by IhorReader or danse_reader 
     229    to new Data2D format 
     230    This is mainly used by the Readers 
     231 
     232    :param data2d: 2d array of Data2D object 
     233    :return: 1d arrays of Data2D object 
     234 
     235    """ 
     236    if data2d.data is None or data2d.x_bins is None or data2d.y_bins is None: 
     237        raise ValueError("Can't convert this data: data=None...") 
     238    new_x = np.tile(data2d.x_bins, (len(data2d.y_bins), 1)) 
     239    new_y = np.tile(data2d.y_bins, (len(data2d.x_bins), 1)) 
     240    new_y = new_y.swapaxes(0, 1) 
     241 
     242    new_data = data2d.data.flatten() 
     243    qx_data = new_x.flatten() 
     244    qy_data = new_y.flatten() 
     245    q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) 
     246    if data2d.err_data is None or np.any(data2d.err_data <= 0): 
     247        new_err_data = np.sqrt(np.abs(new_data)) 
    741248    else: 
    742         if q > q_1 and q <= q_0: 
    743             return (q - q_1) / (q_0 - q_1) 
    744     return None 
    745  
     249        new_err_data = data2d.err_data.flatten() 
     250    mask = np.ones(len(new_data), dtype=bool) 
     251 
     252    # TODO: make sense of the following two lines... 
     253    #from sas.sascalc.dataloader.data_info import Data2D 
     254    #output = Data2D() 
     255    output = data2d 
     256    output.data = new_data 
     257    output.err_data = new_err_data 
     258    output.qx_data = qx_data 
     259    output.qy_data = qy_data 
     260    output.q_data = q_data 
     261    output.mask = mask 
     262 
     263    return output 
     264 
     265################################################################################ 
     266 
     267class _Slab(object): 
     268    """ 
     269    Compute average I(Q) for a region of interest 
     270    """ 
     271 
     272    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, 
     273                 y_max=0.0, bin_width=0.001): 
     274        # Minimum Qx value [A-1] 
     275        self.x_min = x_min 
     276        # Maximum Qx value [A-1] 
     277        self.x_max = x_max 
     278        # Minimum Qy value [A-1] 
     279        self.y_min = y_min 
     280        # Maximum Qy value [A-1] 
     281        self.y_max = y_max 
     282        # Bin width (step size) [A-1] 
     283        self.bin_width = bin_width 
     284        # If True, I(|Q|) will be return, otherwise, 
     285        # negative q-values are allowed 
     286        self.fold = False 
     287 
     288    def __call__(self, data2D): 
     289        return NotImplemented 
     290 
     291    def _avg(self, data2D, maj): 
     292        """ 
     293        Compute average I(Q_maj) for a region of interest. 
     294        The major axis is defined as the axis of Q_maj. 
     295        The minor axis is the axis that we average over. 
     296 
     297        :param data2D: Data2D object 
     298        :param maj_min: min value on the major axis 
     299        :return: Data1D object 
     300        """ 
     301        if len(data2D.detector) > 1: 
     302            msg = "_Slab._avg: invalid number of " 
     303            msg += " detectors: %g" % len(data2D.detector) 
     304            raise RuntimeError(msg) 
     305 
     306        # Get data 
     307        data = data2D.data[np.isfinite(data2D.data)] 
     308        err_data = data2D.err_data[np.isfinite(data2D.data)] 
     309        qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
     310        qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
     311 
     312        # Build array of Q intervals 
     313        if maj == 'x': 
     314            if self.fold: 
     315                x_min = 0 
     316            else: 
     317                x_min = self.x_min 
     318            nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) 
     319        elif maj == 'y': 
     320            if self.fold: 
     321                y_min = 0 
     322            else: 
     323                y_min = self.y_min 
     324            nbins = int(math.ceil((self.y_max - y_min) / self.bin_width)) 
     325        else: 
     326            raise RuntimeError("_Slab._avg: unrecognized axis %s" % str(maj)) 
     327 
     328        x = np.zeros(nbins) 
     329        y = np.zeros(nbins) 
     330        err_y = np.zeros(nbins) 
     331        y_counts = np.zeros(nbins) 
     332 
     333        # Average pixelsize in q space 
     334        for npts in range(len(data)): 
     335            # default frac 
     336            frac_x = 0 
     337            frac_y = 0 
     338            # get ROI 
     339            if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]: 
     340                frac_x = 1 
     341            if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]: 
     342                frac_y = 1 
     343            frac = frac_x * frac_y 
     344 
     345            if frac == 0: 
     346                continue 
     347            # binning: find axis of q 
     348            if maj == 'x': 
     349                q_value = qx_data[npts] 
     350                min_value = x_min 
     351            if maj == 'y': 
     352                q_value = qy_data[npts] 
     353                min_value = y_min 
     354            if self.fold and q_value < 0: 
     355                q_value = -q_value 
     356            # bin 
     357            i_q = int(math.ceil((q_value - min_value) / self.bin_width)) - 1 
     358 
     359            # skip outside of max bins 
     360            if i_q < 0 or i_q >= nbins: 
     361                continue 
     362 
     363            # TODO: find better definition of x[i_q] based on q_data 
     364            # min_value + (i_q + 1) * self.bin_width / 2.0 
     365            x[i_q] += frac * q_value 
     366            y[i_q] += frac * data[npts] 
     367 
     368            if err_data is None or err_data[npts] == 0.0: 
     369                if data[npts] < 0: 
     370                    data[npts] = -data[npts] 
     371                err_y[i_q] += frac * frac * data[npts] 
     372            else: 
     373                err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] 
     374            y_counts[i_q] += frac 
     375 
     376        # Average the sums 
     377        for n in range(nbins): 
     378            err_y[n] = math.sqrt(err_y[n]) 
     379 
     380        err_y = err_y / y_counts 
     381        y = y / y_counts 
     382        x = x / y_counts 
     383        idx = (np.isfinite(y) & np.isfinite(x)) 
     384 
     385        if not idx.any(): 
     386            msg = "Average Error: No points inside ROI to average..." 
     387            raise ValueError(msg) 
     388        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) 
     389 
     390 
     391class SlabY(_Slab): 
     392    """ 
     393    Compute average I(Qy) for a region of interest 
     394    """ 
     395 
     396    def __call__(self, data2D): 
     397        """ 
     398        Compute average I(Qy) for a region of interest 
     399 
     400        :param data2D: Data2D object 
     401        :return: Data1D object 
     402        """ 
     403        return self._avg(data2D, 'y') 
     404 
     405 
     406class SlabX(_Slab): 
     407    """ 
     408    Compute average I(Qx) for a region of interest 
     409    """ 
     410 
     411    def __call__(self, data2D): 
     412        """ 
     413        Compute average I(Qx) for a region of interest 
     414        :param data2D: Data2D object 
     415        :return: Data1D object 
     416        """ 
     417        return self._avg(data2D, 'x') 
     418 
     419################################################################################ 
     420 
     421class Boxsum(object): 
     422    """ 
     423    Perform the sum of counts in a 2D region of interest. 
     424    """ 
     425 
     426    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
     427        # Minimum Qx value [A-1] 
     428        self.x_min = x_min 
     429        # Maximum Qx value [A-1] 
     430        self.x_max = x_max 
     431        # Minimum Qy value [A-1] 
     432        self.y_min = y_min 
     433        # Maximum Qy value [A-1] 
     434        self.y_max = y_max 
     435 
     436    def __call__(self, data2D): 
     437        """ 
     438        Perform the sum in the region of interest 
     439 
     440        :param data2D: Data2D object 
     441        :return: number of counts, error on number of counts, 
     442            number of points summed 
     443        """ 
     444        y, err_y, y_counts = self._sum(data2D) 
     445 
     446        # Average the sums 
     447        counts = 0 if y_counts == 0 else y 
     448        error = 0 if y_counts == 0 else math.sqrt(err_y) 
     449 
     450        # Added y_counts to return, SMK & PDB, 04/03/2013 
     451        return counts, error, y_counts 
     452 
     453    def _sum(self, data2D): 
     454        """ 
     455        Perform the sum in the region of interest 
     456 
     457        :param data2D: Data2D object 
     458        :return: number of counts, 
     459            error on number of counts, number of entries summed 
     460        """ 
     461        if len(data2D.detector) > 1: 
     462            msg = "Circular averaging: invalid number " 
     463            msg += "of detectors: %g" % len(data2D.detector) 
     464            raise RuntimeError(msg) 
     465        # Get data 
     466        data = data2D.data[np.isfinite(data2D.data)] 
     467        err_data = data2D.err_data[np.isfinite(data2D.data)] 
     468        qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
     469        qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
     470 
     471        y = 0.0 
     472        err_y = 0.0 
     473        y_counts = 0.0 
     474 
     475        # Average pixelsize in q space 
     476        for npts in range(len(data)): 
     477            # default frac 
     478            frac_x = 0 
     479            frac_y = 0 
     480 
     481            # get min and max at each points 
     482            qx = qx_data[npts] 
     483            qy = qy_data[npts] 
     484 
     485            # get the ROI 
     486            if self.x_min <= qx and self.x_max > qx: 
     487                frac_x = 1 
     488            if self.y_min <= qy and self.y_max > qy: 
     489                frac_y = 1 
     490            # Find the fraction along each directions 
     491            frac = frac_x * frac_y 
     492            if frac == 0: 
     493                continue 
     494            y += frac * data[npts] 
     495            if err_data is None or err_data[npts] == 0.0: 
     496                if data[npts] < 0: 
     497                    data[npts] = -data[npts] 
     498                err_y += frac * frac * data[npts] 
     499            else: 
     500                err_y += frac * frac * err_data[npts] * err_data[npts] 
     501            y_counts += frac 
     502        return y, err_y, y_counts 
     503 
     504 
     505class Boxavg(Boxsum): 
     506    """ 
     507    Perform the average of counts in a 2D region of interest. 
     508    """ 
     509 
     510    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
     511        super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, 
     512                                     y_min=y_min, y_max=y_max) 
     513 
     514    def __call__(self, data2D): 
     515        """ 
     516        Perform the sum in the region of interest 
     517 
     518        :param data2D: Data2D object 
     519        :return: average counts, error on average counts 
     520 
     521        """ 
     522        y, err_y, y_counts = self._sum(data2D) 
     523 
     524        # Average the sums 
     525        counts = 0 if y_counts == 0 else y / y_counts 
     526        error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts 
     527 
     528        return counts, error 
     529 
     530################################################################################ 
     531 
     532class CircularAverage(object): 
     533    """ 
     534    Perform circular averaging on 2D data 
     535 
     536    The data returned is the distribution of counts 
     537    as a function of Q 
     538    """ 
     539 
     540    def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 
     541        # Minimum radius included in the average [A-1] 
     542        self.r_min = r_min 
     543        # Maximum radius included in the average [A-1] 
     544        self.r_max = r_max 
     545        # Bin width (step size) [A-1] 
     546        self.bin_width = bin_width 
     547 
     548    def __call__(self, data2D, ismask=False): 
     549        """ 
     550        Perform circular averaging on the data 
     551 
     552        :param data2D: Data2D object 
     553        :return: Data1D object 
     554        """ 
     555        # Get data W/ finite values 
     556        data = data2D.data[np.isfinite(data2D.data)] 
     557        q_data = data2D.q_data[np.isfinite(data2D.data)] 
     558        err_data = data2D.err_data[np.isfinite(data2D.data)] 
     559        mask_data = data2D.mask[np.isfinite(data2D.data)] 
     560 
     561        dq_data = None 
     562        if data2D.dqx_data != None and data2D.dqy_data != None: 
     563            dq_data = get_dq_data(data2D) 
     564 
     565        #q_data_max = np.max(q_data) 
     566        if len(data2D.q_data) == None: 
     567            msg = "Circular averaging: invalid q_data: %g" % data2D.q_data 
     568            raise RuntimeError(msg) 
     569 
     570        # Build array of Q intervals 
     571        nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) 
     572 
     573        x = np.zeros(nbins) 
     574        y = np.zeros(nbins) 
     575        err_y = np.zeros(nbins) 
     576        err_x = np.zeros(nbins) 
     577        y_counts = np.zeros(nbins) 
     578 
     579        for npt in range(len(data)): 
     580 
     581            if ismask and not mask_data[npt]: 
     582                continue 
     583 
     584            frac = 0 
     585 
     586            # q-value at the pixel (j,i) 
     587            q_value = q_data[npt] 
     588            data_n = data[npt] 
     589 
     590            # No need to calculate the frac when all data are within range 
     591            if self.r_min >= self.r_max: 
     592                raise ValueError("Limit Error: min > max") 
     593 
     594            if self.r_min <= q_value and q_value <= self.r_max: 
     595                frac = 1 
     596            if frac == 0: 
     597                continue 
     598            i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 
     599 
     600            # Take care of the edge case at phi = 2pi. 
     601            if i_q == nbins: 
     602                i_q = nbins - 1 
     603            y[i_q] += frac * data_n 
     604            # Take dqs from data to get the q_average 
     605            x[i_q] += frac * q_value 
     606            if err_data is None or err_data[npt] == 0.0: 
     607                if data_n < 0: 
     608                    data_n = -data_n 
     609                err_y[i_q] += frac * frac * data_n 
     610            else: 
     611                err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] 
     612            if dq_data != None: 
     613                # To be consistent with dq calculation in 1d reduction, 
     614                # we need just the averages (not quadratures) because 
     615                # it should not depend on the number of the q points 
     616                # in the qr bins. 
     617                err_x[i_q] += frac * dq_data[npt] 
     618            else: 
     619                err_x = None 
     620            y_counts[i_q] += frac 
     621 
     622        # Average the sums 
     623        for n in range(nbins): 
     624            if err_y[n] < 0: 
     625                err_y[n] = -err_y[n] 
     626            err_y[n] = math.sqrt(err_y[n]) 
     627            # if err_x != None: 
     628            #    err_x[n] = math.sqrt(err_x[n]) 
     629 
     630        err_y = err_y / y_counts 
     631        err_y[err_y == 0] = np.average(err_y) 
     632        y = y / y_counts 
     633        x = x / y_counts 
     634        idx = (np.isfinite(y)) & (np.isfinite(x)) 
     635 
     636        if err_x != None: 
     637            d_x = err_x[idx] / y_counts[idx] 
     638        else: 
     639            d_x = None 
     640 
     641        if not idx.any(): 
     642            msg = "Average Error: No points inside ROI to average..." 
     643            raise ValueError(msg) 
     644 
     645        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) 
     646 
     647################################################################################ 
     648 
     649class Ring(object): 
     650    """ 
     651    Defines a ring on a 2D data set. 
     652    The ring is defined by r_min, r_max, and 
     653    the position of the center of the ring. 
     654 
     655    The data returned is the distribution of counts 
     656    around the ring as a function of phi. 
     657 
     658    Phi_min and phi_max should be defined between 0 and 2*pi 
     659    in anti-clockwise starting from the x- axis on the left-hand side 
     660    """ 
     661    # Todo: remove center. 
     662 
     663    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=36): 
     664        # Minimum radius 
     665        self.r_min = r_min 
     666        # Maximum radius 
     667        self.r_max = r_max 
     668        # Center of the ring in x 
     669        self.center_x = center_x 
     670        # Center of the ring in y 
     671        self.center_y = center_y 
     672        # Number of angular bins 
     673        self.nbins_phi = nbins 
     674 
     675    def __call__(self, data2D): 
     676        """ 
     677        Apply the ring to the data set. 
     678        Returns the angular distribution for a given q range 
     679 
     680        :param data2D: Data2D object 
     681 
     682        :return: Data1D object 
     683        """ 
     684        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
     685            raise RuntimeError("Ring averaging only take plottable_2D objects") 
     686 
     687        Pi = math.pi 
     688 
     689        # Get data 
     690        data = data2D.data[np.isfinite(data2D.data)] 
     691        q_data = data2D.q_data[np.isfinite(data2D.data)] 
     692        err_data = data2D.err_data[np.isfinite(data2D.data)] 
     693        qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
     694        qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
     695 
     696        # Set space for 1d outputs 
     697        phi_bins = np.zeros(self.nbins_phi) 
     698        phi_counts = np.zeros(self.nbins_phi) 
     699        phi_values = np.zeros(self.nbins_phi) 
     700        phi_err = np.zeros(self.nbins_phi) 
     701 
     702        # Shift to apply to calculated phi values in order 
     703        # to center first bin at zero 
     704        phi_shift = Pi / self.nbins_phi 
     705 
     706        for npt in range(len(data)): 
     707            frac = 0 
     708            # q-value at the point (npt) 
     709            q_value = q_data[npt] 
     710            data_n = data[npt] 
     711 
     712            # phi-value at the point (npt) 
     713            phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi 
     714 
     715            if self.r_min <= q_value and q_value <= self.r_max: 
     716                frac = 1 
     717            if frac == 0: 
     718                continue 
     719            # binning 
     720            i_phi = int(math.floor((self.nbins_phi) * 
     721                                   (phi_value + phi_shift) / (2 * Pi))) 
     722 
     723            # Take care of the edge case at phi = 2pi. 
     724            if i_phi >= self.nbins_phi: 
     725                i_phi = 0 
     726            phi_bins[i_phi] += frac * data[npt] 
     727 
     728            if err_data is None or err_data[npt] == 0.0: 
     729                if data_n < 0: 
     730                    data_n = -data_n 
     731                phi_err[i_phi] += frac * frac * math.fabs(data_n) 
     732            else: 
     733                phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt] 
     734            phi_counts[i_phi] += frac 
     735 
     736        for i in range(self.nbins_phi): 
     737            phi_bins[i] = phi_bins[i] / phi_counts[i] 
     738            phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i] 
     739            phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) 
     740 
     741        idx = (np.isfinite(phi_bins)) 
     742 
     743        if not idx.any(): 
     744            msg = "Average Error: No points inside ROI to average..." 
     745            raise ValueError(msg) 
     746        # elif len(phi_bins[idx])!= self.nbins_phi: 
     747        #    print "resulted",self.nbins_phi- len(phi_bins[idx]) 
     748        #,"empty bin(s) due to tight binning..." 
     749        return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) 
     750 
     751################################################################################ 
    746752 
    747753class _Sector(object): 
    … …  
    776782        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    777783            raise RuntimeError("Ring averaging only take plottable_2D objects") 
    778         Pi = math.pi 
    779784 
    780785        # Get the all data & info 
    … …  
    784789        qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    785790        qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
     791 
    786792        dq_data = None 
    787  
    788         # Get the dq for resolution averaging 
    789793        if data2D.dqx_data != None and data2D.dqy_data != None: 
    790             # The pinholes and det. pix contribution present 
    791             # in both direction of the 2D which must be subtracted when 
    792             # converting to 1D: dq_overlap should calculated ideally at 
    793             # q = 0. 
    794             # Extrapolate dqy(perp) at q = 0 
    795             z_max = max(data2D.q_data) 
    796             z_min = min(data2D.q_data) 
    797             x_max = data2D.dqx_data[data2D.q_data[z_max]] 
    798             x_min = data2D.dqx_data[data2D.q_data[z_min]] 
    799             y_max = data2D.dqy_data[data2D.q_data[z_max]] 
    800             y_min = data2D.dqy_data[data2D.q_data[z_min]] 
    801             # Find qdx at q = 0 
    802             dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
    803             # when extrapolation goes wrong 
    804             if dq_overlap_x > min(data2D.dqx_data): 
    805                 dq_overlap_x = min(data2D.dqx_data) 
    806             dq_overlap_x *= dq_overlap_x 
    807             # Find qdx at q = 0 
    808             dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
    809             # when extrapolation goes wrong 
    810             if dq_overlap_y > min(data2D.dqy_data): 
    811                 dq_overlap_y = min(data2D.dqy_data) 
    812             # get dq at q=0. 
    813             dq_overlap_y *= dq_overlap_y 
    814  
    815             dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 
    816             if dq_overlap < 0: 
    817                 dq_overlap = y_min 
    818             dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] 
    819             dqy_data = data2D.dqy_data[np.isfinite( 
    820                 data2D.data)] - dq_overlap 
    821             # def; dqx_data = dq_r dqy_data = dq_phi 
    822             # Convert dq 2D to 1D here 
    823             dqx = dqx_data * dqx_data 
    824             dqy = dqy_data * dqy_data 
    825             dq_data = np.add(dqx, dqy) 
    826             dq_data = np.sqrt(dq_data) 
     794            dq_data = get_dq_data(data2D) 
    827795 
    828796        # set space for 1d outputs 
    … …  
    838806 
    839807        for n in range(len(data)): 
    840             frac = 0 
    841808 
    842809            # q-value at the pixel (j,i) 
    … …  
    848815 
    849816            # phi-value of the pixel (j,i) 
    850             phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 
     817            phi_value = math.atan2(qy_data[n], qx_data[n]) + math.pi 
    851818 
    852819            # No need to calculate the frac when all data are within range 
    853             if self.r_min <= q_value and q_value <= self.r_max: 
    854                 frac = 1 
    855             if frac == 0: 
     820            if self.r_min > q_value or q_value > self.r_max: 
    856821                continue 
     822 
    857823            # In case of two ROIs (symmetric major and minor regions)(for 'q2') 
    858824            if run.lower() == 'q2': 
    859825                # For minor sector wing 
    860826                # Calculate the minor wing phis 
    861                 phi_min_minor = flip_phi(phi_min - Pi) 
    862                 phi_max_minor = flip_phi(phi_max - Pi) 
     827                phi_min_minor = flip_phi(phi_min - math.pi) 
     828                phi_max_minor = flip_phi(phi_max - math.pi) 
    863829                # Check if phis of the minor ring is within 0 to 2pi 
    864830                if phi_min_minor > phi_max_minor: 
    … …  
    879845 
    880846            if not is_in: 
    881                 frac = 0 
    882             if frac == 0: 
    883847                continue 
    884848            # Check which type of averaging we need 
    … …  
    897861 
    898862            # Get the total y 
    899             y[i_bin] += frac * data_n 
    900             x[i_bin] += frac * q_value 
     863            y[i_bin] += data_n 
     864            x[i_bin] += q_value 
    901865            if err_data[n] == None or err_data[n] == 0.0: 
    902866                if data_n < 0: 
    903867                    data_n = -data_n 
    904                 y_err[i_bin] += frac * frac * data_n 
     868                y_err[i_bin] += data_n 
    905869            else: 
    906                 y_err[i_bin] += frac * frac * err_data[n] * err_data[n] 
     870                y_err[i_bin] += err_data[n] * err_data[n] 
    907871 
    908872            if dq_data != None: 
    … …  
    911875                # it should not depend on the number of the q points 
    912876                # in the qr bins. 
    913                 x_err[i_bin] += frac * dq_data[n] 
     877                x_err[i_bin] += dq_data[n] 
    914878            else: 
    915879                x_err = None 
    916             y_counts[i_bin] += frac 
     880            y_counts[i_bin] += 1 
    917881 
    918882        # Organize the results 
    … …  
    988952        return self._agv(data2D, 'q2') 
    989953 
     954################################################################################ 
    990955 
    991956class Ringcut(object): 
    … …  
    1032997        return out 
    1033998 
     999################################################################################ 
    10341000 
    10351001class Boxcut(object): 
    … …  
    10811047        return outx & outy 
    10821048 
     1049################################################################################ 
    10831050 
    10841051class Sectorcut(object): 
Note: See TracChangeset for help on using the changeset viewer.