[f60a8c2] | 1 | """ |
---|
| 2 | Data manipulations for 2D data sets. |
---|
| 3 | Using the meta data information, various types of averaging |
---|
| 4 | are performed in Q-space |
---|
| 5 | """ |
---|
[0997158f] | 6 | ##################################################################### |
---|
| 7 | #This software was developed by the University of Tennessee as part of the |
---|
| 8 | #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
[f60a8c2] | 9 | #project funded by the US National Science Foundation. |
---|
[0997158f] | 10 | #See the license text in license.txt |
---|
| 11 | #copyright 2008, University of Tennessee |
---|
| 12 | ###################################################################### |
---|
| 13 | |
---|
[76e2369] | 14 | #TODO: copy the meta data from the 2D object to the resulting 1D object |
---|
| 15 | import math |
---|
[9a5097c] | 16 | import numpy as np |
---|
[76e2369] | 17 | |
---|
[a7a5886] | 18 | #from data_info import plottable_2D |
---|
| 19 | from data_info import Data1D |
---|
| 20 | |
---|
| 21 | |
---|
[76e2369] | 22 | def get_q(dx, dy, det_dist, wavelength): |
---|
| 23 | """ |
---|
[0997158f] | 24 | :param dx: x-distance from beam center [mm] |
---|
| 25 | :param dy: y-distance from beam center [mm] |
---|
| 26 | :return: q-value at the given position |
---|
[76e2369] | 27 | """ |
---|
| 28 | # Distance from beam center in the plane of detector |
---|
[c8a6c3d7] | 29 | plane_dist = math.sqrt(dx * dx + dy * dy) |
---|
[76e2369] | 30 | # Half of the scattering angle |
---|
[c8a6c3d7] | 31 | theta = 0.5 * math.atan(plane_dist / det_dist) |
---|
| 32 | return (4.0 * math.pi / wavelength) * math.sin(theta) |
---|
[acb37d9] | 33 | |
---|
[f60a8c2] | 34 | |
---|
[a7a5886] | 35 | def get_q_compo(dx, dy, det_dist, wavelength, compo=None): |
---|
[0997158f] | 36 | """ |
---|
| 37 | This reduces tiny error at very large q. |
---|
| 38 | Implementation of this func is not started yet.<--ToDo |
---|
| 39 | """ |
---|
[a7a5886] | 40 | if dy == 0: |
---|
| 41 | if dx >= 0: |
---|
| 42 | angle_xy = 0 |
---|
[acb37d9] | 43 | else: |
---|
[a7a5886] | 44 | angle_xy = math.pi |
---|
[acb37d9] | 45 | else: |
---|
[c8a6c3d7] | 46 | angle_xy = math.atan(dx / dy) |
---|
| 47 | |
---|
[a7a5886] | 48 | if compo == "x": |
---|
| 49 | out = get_q(dx, dy, det_dist, wavelength) * math.cos(angle_xy) |
---|
| 50 | elif compo == "y": |
---|
| 51 | out = get_q(dx, dy, det_dist, wavelength) * math.sin(angle_xy) |
---|
[acb37d9] | 52 | else: |
---|
[a7a5886] | 53 | out = get_q(dx, dy, det_dist, wavelength) |
---|
[acb37d9] | 54 | return out |
---|
[095ab1b] | 55 | |
---|
[f60a8c2] | 56 | |
---|
[095ab1b] | 57 | def flip_phi(phi): |
---|
| 58 | """ |
---|
[0997158f] | 59 | Correct phi to within the 0 <= to <= 2pi range |
---|
[c8a6c3d7] | 60 | |
---|
[0997158f] | 61 | :return: phi in >=0 and <=2Pi |
---|
[095ab1b] | 62 | """ |
---|
| 63 | Pi = math.pi |
---|
| 64 | if phi < 0: |
---|
[f60a8c2] | 65 | phi_out = phi + (2 * Pi) |
---|
[a7a5886] | 66 | elif phi > (2 * Pi): |
---|
[f60a8c2] | 67 | phi_out = phi - (2 * Pi) |
---|
[095ab1b] | 68 | else: |
---|
[f60a8c2] | 69 | phi_out = phi |
---|
[095ab1b] | 70 | return phi_out |
---|
| 71 | |
---|
[f60a8c2] | 72 | |
---|
[095ab1b] | 73 | def reader2D_converter(data2d=None): |
---|
| 74 | """ |
---|
[a7a5886] | 75 | convert old 2d format opened by IhorReader or danse_reader |
---|
| 76 | to new Data2D format |
---|
[c8a6c3d7] | 77 | |
---|
[0997158f] | 78 | :param data2d: 2d array of Data2D object |
---|
| 79 | :return: 1d arrays of Data2D object |
---|
[c8a6c3d7] | 80 | |
---|
[095ab1b] | 81 | """ |
---|
[36d69e1] | 82 | if data2d.data is None or data2d.x_bins is None or data2d.y_bins is None: |
---|
[a7a5886] | 83 | raise ValueError, "Can't convert this data: data=None..." |
---|
[9a5097c] | 84 | new_x = np.tile(data2d.x_bins, (len(data2d.y_bins), 1)) |
---|
| 85 | new_y = np.tile(data2d.y_bins, (len(data2d.x_bins), 1)) |
---|
[a7a5886] | 86 | new_y = new_y.swapaxes(0, 1) |
---|
[095ab1b] | 87 | |
---|
| 88 | new_data = data2d.data.flatten() |
---|
| 89 | qx_data = new_x.flatten() |
---|
| 90 | qy_data = new_y.flatten() |
---|
[9a5097c] | 91 | q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) |
---|
| 92 | if data2d.err_data == None or np.any(data2d.err_data <= 0): |
---|
| 93 | new_err_data = np.sqrt(np.abs(new_data)) |
---|
[ed2276f] | 94 | |
---|
[dde2d44] | 95 | else: |
---|
| 96 | new_err_data = data2d.err_data.flatten() |
---|
[9a5097c] | 97 | mask = np.ones(len(new_data), dtype=bool) |
---|
[095ab1b] | 98 | |
---|
[f60a8c2] | 99 | #TODO: make sense of the following two lines... |
---|
[b699768] | 100 | #from sas.sascalc.dataloader.data_info import Data2D |
---|
[c8a6c3d7] | 101 | #output = Data2D() |
---|
[095ab1b] | 102 | output = data2d |
---|
| 103 | output.data = new_data |
---|
| 104 | output.err_data = new_err_data |
---|
| 105 | output.qx_data = qx_data |
---|
| 106 | output.qy_data = qy_data |
---|
| 107 | output.q_data = q_data |
---|
| 108 | output.mask = mask |
---|
| 109 | |
---|
| 110 | return output |
---|
| 111 | |
---|
[f60a8c2] | 112 | |
---|
[70975f3] | 113 | class _Slab(object): |
---|
| 114 | """ |
---|
[0997158f] | 115 | Compute average I(Q) for a region of interest |
---|
[70975f3] | 116 | """ |
---|
[a7a5886] | 117 | def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, |
---|
| 118 | y_max=0.0, bin_width=0.001): |
---|
[70975f3] | 119 | # Minimum Qx value [A-1] |
---|
| 120 | self.x_min = x_min |
---|
| 121 | # Maximum Qx value [A-1] |
---|
| 122 | self.x_max = x_max |
---|
| 123 | # Minimum Qy value [A-1] |
---|
| 124 | self.y_min = y_min |
---|
| 125 | # Maximum Qy value [A-1] |
---|
| 126 | self.y_max = y_max |
---|
| 127 | # Bin width (step size) [A-1] |
---|
| 128 | self.bin_width = bin_width |
---|
[a7a5886] | 129 | # If True, I(|Q|) will be return, otherwise, |
---|
| 130 | # negative q-values are allowed |
---|
[70975f3] | 131 | self.fold = False |
---|
[c8a6c3d7] | 132 | |
---|
[a7a5886] | 133 | def __call__(self, data2D): |
---|
| 134 | return NotImplemented |
---|
[c8a6c3d7] | 135 | |
---|
[70975f3] | 136 | def _avg(self, data2D, maj): |
---|
| 137 | """ |
---|
[0997158f] | 138 | Compute average I(Q_maj) for a region of interest. |
---|
| 139 | The major axis is defined as the axis of Q_maj. |
---|
| 140 | The minor axis is the axis that we average over. |
---|
[c8a6c3d7] | 141 | |
---|
[0997158f] | 142 | :param data2D: Data2D object |
---|
| 143 | :param maj_min: min value on the major axis |
---|
| 144 | :return: Data1D object |
---|
[70975f3] | 145 | """ |
---|
[b2b36932] | 146 | if len(data2D.detector) > 1: |
---|
[a7a5886] | 147 | msg = "_Slab._avg: invalid number of " |
---|
| 148 | msg += " detectors: %g" % len(data2D.detector) |
---|
| 149 | raise RuntimeError, msg |
---|
[c8a6c3d7] | 150 | |
---|
[f60a8c2] | 151 | # Get data |
---|
[9a5097c] | 152 | data = data2D.data[np.isfinite(data2D.data)] |
---|
| 153 | err_data = data2D.err_data[np.isfinite(data2D.data)] |
---|
| 154 | qx_data = data2D.qx_data[np.isfinite(data2D.data)] |
---|
| 155 | qy_data = data2D.qy_data[np.isfinite(data2D.data)] |
---|
[c8a6c3d7] | 156 | |
---|
[70975f3] | 157 | # Build array of Q intervals |
---|
[a7a5886] | 158 | if maj == 'x': |
---|
| 159 | if self.fold: |
---|
[f60a8c2] | 160 | x_min = 0 |
---|
| 161 | else: |
---|
| 162 | x_min = self.x_min |
---|
| 163 | nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) |
---|
[a7a5886] | 164 | elif maj == 'y': |
---|
[f60a8c2] | 165 | if self.fold: |
---|
| 166 | y_min = 0 |
---|
| 167 | else: |
---|
| 168 | y_min = self.y_min |
---|
[c8a6c3d7] | 169 | nbins = int(math.ceil((self.y_max - y_min) / self.bin_width)) |
---|
[70975f3] | 170 | else: |
---|
| 171 | raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj) |
---|
[c8a6c3d7] | 172 | |
---|
[9a5097c] | 173 | x = np.zeros(nbins) |
---|
| 174 | y = np.zeros(nbins) |
---|
| 175 | err_y = np.zeros(nbins) |
---|
| 176 | y_counts = np.zeros(nbins) |
---|
[70975f3] | 177 | |
---|
[f60a8c2] | 178 | # Average pixelsize in q space |
---|
| 179 | for npts in range(len(data)): |
---|
| 180 | # default frac |
---|
[095ab1b] | 181 | frac_x = 0 |
---|
| 182 | frac_y = 0 |
---|
| 183 | # get ROI |
---|
| 184 | if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]: |
---|
| 185 | frac_x = 1 |
---|
| 186 | if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]: |
---|
| 187 | frac_y = 1 |
---|
| 188 | frac = frac_x * frac_y |
---|
[c8a6c3d7] | 189 | |
---|
[a7a5886] | 190 | if frac == 0: |
---|
| 191 | continue |
---|
[095ab1b] | 192 | # binning: find axis of q |
---|
[f60a8c2] | 193 | if maj == 'x': |
---|
[095ab1b] | 194 | q_value = qx_data[npts] |
---|
[c8a6c3d7] | 195 | min_value = x_min |
---|
[f60a8c2] | 196 | if maj == 'y': |
---|
| 197 | q_value = qy_data[npts] |
---|
[c8a6c3d7] | 198 | min_value = y_min |
---|
[a7a5886] | 199 | if self.fold and q_value < 0: |
---|
[f60a8c2] | 200 | q_value = -q_value |
---|
[095ab1b] | 201 | # bin |
---|
[c8a6c3d7] | 202 | i_q = int(math.ceil((q_value - min_value) / self.bin_width)) - 1 |
---|
| 203 | |
---|
[095ab1b] | 204 | # skip outside of max bins |
---|
[a7a5886] | 205 | if i_q < 0 or i_q >= nbins: |
---|
| 206 | continue |
---|
[c8a6c3d7] | 207 | |
---|
[095ab1b] | 208 | #TODO: find better definition of x[i_q] based on q_data |
---|
[3c3a440] | 209 | # min_value + (i_q + 1) * self.bin_width / 2.0 |
---|
| 210 | x[i_q] += frac * q_value |
---|
[a7a5886] | 211 | y[i_q] += frac * data[npts] |
---|
[c8a6c3d7] | 212 | |
---|
[a7a5886] | 213 | if err_data == None or err_data[npts] == 0.0: |
---|
[f60a8c2] | 214 | if data[npts] < 0: |
---|
| 215 | data[npts] = -data[npts] |
---|
[c6f95bb] | 216 | err_y[i_q] += frac * frac * data[npts] |
---|
[095ab1b] | 217 | else: |
---|
| 218 | err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] |
---|
[f60a8c2] | 219 | y_counts[i_q] += frac |
---|
[c8a6c3d7] | 220 | |
---|
[f60a8c2] | 221 | # Average the sums |
---|
[095ab1b] | 222 | for n in range(nbins): |
---|
| 223 | err_y[n] = math.sqrt(err_y[n]) |
---|
[c8a6c3d7] | 224 | |
---|
[a7a5886] | 225 | err_y = err_y / y_counts |
---|
[f60a8c2] | 226 | y = y / y_counts |
---|
| 227 | x = x / y_counts |
---|
[9a5097c] | 228 | idx = (np.isfinite(y) & np.isfinite(x)) |
---|
[c8a6c3d7] | 229 | |
---|
| 230 | if not idx.any(): |
---|
[f60a8c2] | 231 | msg = "Average Error: No points inside ROI to average..." |
---|
[a7a5886] | 232 | raise ValueError, msg |
---|
[095ab1b] | 233 | return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) |
---|
[c8a6c3d7] | 234 | |
---|
| 235 | |
---|
[70975f3] | 236 | class SlabY(_Slab): |
---|
| 237 | """ |
---|
[0997158f] | 238 | Compute average I(Qy) for a region of interest |
---|
[70975f3] | 239 | """ |
---|
| 240 | def __call__(self, data2D): |
---|
| 241 | """ |
---|
[0997158f] | 242 | Compute average I(Qy) for a region of interest |
---|
[c8a6c3d7] | 243 | |
---|
[0997158f] | 244 | :param data2D: Data2D object |
---|
| 245 | :return: Data1D object |
---|
[70975f3] | 246 | """ |
---|
| 247 | return self._avg(data2D, 'y') |
---|
[c8a6c3d7] | 248 | |
---|
| 249 | |
---|
[70975f3] | 250 | class SlabX(_Slab): |
---|
| 251 | """ |
---|
[0997158f] | 252 | Compute average I(Qx) for a region of interest |
---|
[70975f3] | 253 | """ |
---|
| 254 | def __call__(self, data2D): |
---|
| 255 | """ |
---|
[0997158f] | 256 | Compute average I(Qx) for a region of interest |
---|
| 257 | :param data2D: Data2D object |
---|
| 258 | :return: Data1D object |
---|
[70975f3] | 259 | """ |
---|
[f60a8c2] | 260 | return self._avg(data2D, 'x') |
---|
| 261 | |
---|
| 262 | |
---|
[f8d0ee7] | 263 | class Boxsum(object): |
---|
| 264 | """ |
---|
[0997158f] | 265 | Perform the sum of counts in a 2D region of interest. |
---|
[f8d0ee7] | 266 | """ |
---|
| 267 | def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): |
---|
| 268 | # Minimum Qx value [A-1] |
---|
| 269 | self.x_min = x_min |
---|
| 270 | # Maximum Qx value [A-1] |
---|
| 271 | self.x_max = x_max |
---|
| 272 | # Minimum Qy value [A-1] |
---|
| 273 | self.y_min = y_min |
---|
| 274 | # Maximum Qy value [A-1] |
---|
| 275 | self.y_max = y_max |
---|
| 276 | |
---|
| 277 | def __call__(self, data2D): |
---|
| 278 | """ |
---|
[f60a8c2] | 279 | Perform the sum in the region of interest |
---|
[c8a6c3d7] | 280 | |
---|
[0997158f] | 281 | :param data2D: Data2D object |
---|
[d555416] | 282 | :return: number of counts, error on number of counts, |
---|
| 283 | number of points summed |
---|
[f8d0ee7] | 284 | """ |
---|
| 285 | y, err_y, y_counts = self._sum(data2D) |
---|
[c8a6c3d7] | 286 | |
---|
[f8d0ee7] | 287 | # Average the sums |
---|
[a7a5886] | 288 | counts = 0 if y_counts == 0 else y |
---|
[f60a8c2] | 289 | error = 0 if y_counts == 0 else math.sqrt(err_y) |
---|
[c8a6c3d7] | 290 | |
---|
[d555416] | 291 | # Added y_counts to return, SMK & PDB, 04/03/2013 |
---|
| 292 | return counts, error, y_counts |
---|
[c8a6c3d7] | 293 | |
---|
[f8d0ee7] | 294 | def _sum(self, data2D): |
---|
| 295 | """ |
---|
[f60a8c2] | 296 | Perform the sum in the region of interest |
---|
[c8a6c3d7] | 297 | |
---|
[0997158f] | 298 | :param data2D: Data2D object |
---|
[f60a8c2] | 299 | :return: number of counts, |
---|
[a7a5886] | 300 | error on number of counts, number of entries summed |
---|
[f8d0ee7] | 301 | """ |
---|
[b2b36932] | 302 | if len(data2D.detector) > 1: |
---|
[a7a5886] | 303 | msg = "Circular averaging: invalid number " |
---|
| 304 | msg += "of detectors: %g" % len(data2D.detector) |
---|
| 305 | raise RuntimeError, msg |
---|
[f60a8c2] | 306 | # Get data |
---|
[9a5097c] | 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)] |
---|
[c8a6c3d7] | 311 | |
---|
[f60a8c2] | 312 | y = 0.0 |
---|
[f8d0ee7] | 313 | err_y = 0.0 |
---|
| 314 | y_counts = 0.0 |
---|
| 315 | |
---|
[f60a8c2] | 316 | # Average pixelsize in q space |
---|
| 317 | for npts in range(len(data)): |
---|
| 318 | # default frac |
---|
| 319 | frac_x = 0 |
---|
| 320 | frac_y = 0 |
---|
[c8a6c3d7] | 321 | |
---|
[095ab1b] | 322 | # get min and max at each points |
---|
| 323 | qx = qx_data[npts] |
---|
| 324 | qy = qy_data[npts] |
---|
[c8a6c3d7] | 325 | |
---|
[095ab1b] | 326 | # get the ROI |
---|
| 327 | if self.x_min <= qx and self.x_max > qx: |
---|
| 328 | frac_x = 1 |
---|
| 329 | if self.y_min <= qy and self.y_max > qy: |
---|
| 330 | frac_y = 1 |
---|
[f60a8c2] | 331 | #Find the fraction along each directions |
---|
[095ab1b] | 332 | frac = frac_x * frac_y |
---|
[a7a5886] | 333 | if frac == 0: |
---|
| 334 | continue |
---|
[095ab1b] | 335 | y += frac * data[npts] |
---|
[a7a5886] | 336 | if err_data == None or err_data[npts] == 0.0: |
---|
| 337 | if data[npts] < 0: |
---|
| 338 | data[npts] = -data[npts] |
---|
[c6f95bb] | 339 | err_y += frac * frac * data[npts] |
---|
[095ab1b] | 340 | else: |
---|
| 341 | err_y += frac * frac * err_data[npts] * err_data[npts] |
---|
[f60a8c2] | 342 | y_counts += frac |
---|
[f8d0ee7] | 343 | return y, err_y, y_counts |
---|
[095ab1b] | 344 | |
---|
| 345 | |
---|
[f8d0ee7] | 346 | class Boxavg(Boxsum): |
---|
| 347 | """ |
---|
[0997158f] | 348 | Perform the average of counts in a 2D region of interest. |
---|
[f8d0ee7] | 349 | """ |
---|
| 350 | def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): |
---|
[a7a5886] | 351 | super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, |
---|
[c8a6c3d7] | 352 | y_min=y_min, y_max=y_max) |
---|
[f8d0ee7] | 353 | |
---|
| 354 | def __call__(self, data2D): |
---|
| 355 | """ |
---|
[f60a8c2] | 356 | Perform the sum in the region of interest |
---|
[c8a6c3d7] | 357 | |
---|
[0997158f] | 358 | :param data2D: Data2D object |
---|
| 359 | :return: average counts, error on average counts |
---|
[c8a6c3d7] | 360 | |
---|
[f8d0ee7] | 361 | """ |
---|
| 362 | y, err_y, y_counts = self._sum(data2D) |
---|
[c8a6c3d7] | 363 | |
---|
[f8d0ee7] | 364 | # Average the sums |
---|
[f60a8c2] | 365 | counts = 0 if y_counts == 0 else y / y_counts |
---|
| 366 | error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts |
---|
[c8a6c3d7] | 367 | |
---|
[f8d0ee7] | 368 | return counts, error |
---|
[c8a6c3d7] | 369 | |
---|
| 370 | |
---|
[f8d0ee7] | 371 | def get_pixel_fraction_square(x, xmin, xmax): |
---|
| 372 | """ |
---|
[f60a8c2] | 373 | Return the fraction of the length |
---|
[0997158f] | 374 | from xmin to x.:: |
---|
[c8a6c3d7] | 375 | |
---|
[0997158f] | 376 | A B |
---|
| 377 | +-----------+---------+ |
---|
| 378 | xmin x xmax |
---|
[c8a6c3d7] | 379 | |
---|
[0997158f] | 380 | :param x: x-value |
---|
| 381 | :param xmin: minimum x for the length considered |
---|
| 382 | :param xmax: minimum x for the length considered |
---|
| 383 | :return: (x-xmin)/(xmax-xmin) when xmin < x < xmax |
---|
[c8a6c3d7] | 384 | |
---|
[f8d0ee7] | 385 | """ |
---|
[a7a5886] | 386 | if x <= xmin: |
---|
[f8d0ee7] | 387 | return 0.0 |
---|
[a7a5886] | 388 | if x > xmin and x < xmax: |
---|
| 389 | return (x - xmin) / (xmax - xmin) |
---|
[f8d0ee7] | 390 | else: |
---|
| 391 | return 1.0 |
---|
| 392 | |
---|
[76e2369] | 393 | |
---|
| 394 | class CircularAverage(object): |
---|
| 395 | """ |
---|
[0997158f] | 396 | Perform circular averaging on 2D data |
---|
[c8a6c3d7] | 397 | |
---|
[0997158f] | 398 | The data returned is the distribution of counts |
---|
| 399 | as a function of Q |
---|
[76e2369] | 400 | """ |
---|
[095ab1b] | 401 | def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): |
---|
[76e2369] | 402 | # Minimum radius included in the average [A-1] |
---|
| 403 | self.r_min = r_min |
---|
| 404 | # Maximum radius included in the average [A-1] |
---|
| 405 | self.r_max = r_max |
---|
| 406 | # Bin width (step size) [A-1] |
---|
| 407 | self.bin_width = bin_width |
---|
| 408 | |
---|
[8f12385] | 409 | def __call__(self, data2D, ismask=False): |
---|
[76e2369] | 410 | """ |
---|
[0997158f] | 411 | Perform circular averaging on the data |
---|
[c8a6c3d7] | 412 | |
---|
[0997158f] | 413 | :param data2D: Data2D object |
---|
| 414 | :return: Data1D object |
---|
[76e2369] | 415 | """ |
---|
[729bcf6] | 416 | # Get data W/ finite values |
---|
[9a5097c] | 417 | data = data2D.data[np.isfinite(data2D.data)] |
---|
| 418 | q_data = data2D.q_data[np.isfinite(data2D.data)] |
---|
| 419 | err_data = data2D.err_data[np.isfinite(data2D.data)] |
---|
| 420 | mask_data = data2D.mask[np.isfinite(data2D.data)] |
---|
[c8a6c3d7] | 421 | |
---|
[342a506] | 422 | dq_data = None |
---|
[c8a6c3d7] | 423 | |
---|
[729bcf6] | 424 | # Get the dq for resolution averaging |
---|
[342a506] | 425 | if data2D.dqx_data != None and data2D.dqy_data != None: |
---|
[f60a8c2] | 426 | # The pinholes and det. pix contribution present |
---|
[729bcf6] | 427 | # in both direction of the 2D which must be subtracted when |
---|
| 428 | # converting to 1D: dq_overlap should calculated ideally at |
---|
[f60a8c2] | 429 | # q = 0. Note This method works on only pinhole geometry. |
---|
[729bcf6] | 430 | # Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. |
---|
| 431 | z_max = max(data2D.q_data) |
---|
| 432 | z_min = min(data2D.q_data) |
---|
| 433 | x_max = data2D.dqx_data[data2D.q_data[z_max]] |
---|
[f60a8c2] | 434 | x_min = data2D.dqx_data[data2D.q_data[z_min]] |
---|
[729bcf6] | 435 | y_max = data2D.dqy_data[data2D.q_data[z_max]] |
---|
[f60a8c2] | 436 | y_min = data2D.dqy_data[data2D.q_data[z_min]] |
---|
[729bcf6] | 437 | # Find qdx at q = 0 |
---|
| 438 | dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) |
---|
| 439 | # when extrapolation goes wrong |
---|
| 440 | if dq_overlap_x > min(data2D.dqx_data): |
---|
| 441 | dq_overlap_x = min(data2D.dqx_data) |
---|
[f60a8c2] | 442 | dq_overlap_x *= dq_overlap_x |
---|
[729bcf6] | 443 | # Find qdx at q = 0 |
---|
| 444 | dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) |
---|
| 445 | # when extrapolation goes wrong |
---|
| 446 | if dq_overlap_y > min(data2D.dqy_data): |
---|
| 447 | dq_overlap_y = min(data2D.dqy_data) |
---|
| 448 | # get dq at q=0. |
---|
| 449 | dq_overlap_y *= dq_overlap_y |
---|
| 450 | |
---|
[9a5097c] | 451 | dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) |
---|
[729bcf6] | 452 | # Final protection of dq |
---|
| 453 | if dq_overlap < 0: |
---|
| 454 | dq_overlap = y_min |
---|
[9a5097c] | 455 | dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] |
---|
| 456 | dqy_data = data2D.dqy_data[np.isfinite(data2D.data)] - dq_overlap |
---|
[729bcf6] | 457 | # def; dqx_data = dq_r dqy_data = dq_phi |
---|
| 458 | # Convert dq 2D to 1D here |
---|
[f60a8c2] | 459 | dqx = dqx_data * dqx_data |
---|
[729bcf6] | 460 | dqy = dqy_data * dqy_data |
---|
[9a5097c] | 461 | dq_data = np.add(dqx, dqy) |
---|
| 462 | dq_data = np.sqrt(dq_data) |
---|
[c8a6c3d7] | 463 | |
---|
[9a5097c] | 464 | #q_data_max = np.max(q_data) |
---|
[095ab1b] | 465 | if len(data2D.q_data) == None: |
---|
[a7a5886] | 466 | msg = "Circular averaging: invalid q_data: %g" % data2D.q_data |
---|
| 467 | raise RuntimeError, msg |
---|
[095ab1b] | 468 | |
---|
[76e2369] | 469 | # Build array of Q intervals |
---|
[a7a5886] | 470 | nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) |
---|
[095ab1b] | 471 | |
---|
[9a5097c] | 472 | x = np.zeros(nbins) |
---|
| 473 | y = np.zeros(nbins) |
---|
| 474 | err_y = np.zeros(nbins) |
---|
| 475 | err_x = np.zeros(nbins) |
---|
| 476 | y_counts = np.zeros(nbins) |
---|
[095ab1b] | 477 | |
---|
[f60a8c2] | 478 | for npt in range(len(data)): |
---|
[c8a6c3d7] | 479 | |
---|
[8f12385] | 480 | if ismask and not mask_data[npt]: |
---|
[f60a8c2] | 481 | continue |
---|
[c8a6c3d7] | 482 | |
---|
[095ab1b] | 483 | frac = 0 |
---|
[c8a6c3d7] | 484 | |
---|
[095ab1b] | 485 | # q-value at the pixel (j,i) |
---|
[f60a8c2] | 486 | q_value = q_data[npt] |
---|
| 487 | data_n = data[npt] |
---|
[c8a6c3d7] | 488 | |
---|
[095ab1b] | 489 | ## No need to calculate the frac when all data are within range |
---|
| 490 | if self.r_min >= self.r_max: |
---|
[f60a8c2] | 491 | raise ValueError, "Limit Error: min > max" |
---|
[c8a6c3d7] | 492 | |
---|
[a7a5886] | 493 | if self.r_min <= q_value and q_value <= self.r_max: |
---|
[f60a8c2] | 494 | frac = 1 |
---|
[a7a5886] | 495 | if frac == 0: |
---|
[c8a6c3d7] | 496 | continue |
---|
[f60a8c2] | 497 | i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) |
---|
[095ab1b] | 498 | |
---|
[f60a8c2] | 499 | # Take care of the edge case at phi = 2pi. |
---|
| 500 | if i_q == nbins: |
---|
| 501 | i_q = nbins - 1 |
---|
[095ab1b] | 502 | y[i_q] += frac * data_n |
---|
[729bcf6] | 503 | # Take dqs from data to get the q_average |
---|
| 504 | x[i_q] += frac * q_value |
---|
[a7a5886] | 505 | if err_data == None or err_data[npt] == 0.0: |
---|
| 506 | if data_n < 0: |
---|
| 507 | data_n = -data_n |
---|
[c6f95bb] | 508 | err_y[i_q] += frac * frac * data_n |
---|
[8ba103f] | 509 | else: |
---|
[095ab1b] | 510 | err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] |
---|
[342a506] | 511 | if dq_data != None: |
---|
[f60a8c2] | 512 | # To be consistent with dq calculation in 1d reduction, |
---|
| 513 | # we need just the averages (not quadratures) because |
---|
| 514 | # it should not depend on the number of the q points |
---|
[729bcf6] | 515 | # in the qr bins. |
---|
| 516 | err_x[i_q] += frac * dq_data[npt] |
---|
[342a506] | 517 | else: |
---|
| 518 | err_x = None |
---|
[f60a8c2] | 519 | y_counts[i_q] += frac |
---|
[c8a6c3d7] | 520 | |
---|
[f60a8c2] | 521 | # Average the sums |
---|
[095ab1b] | 522 | for n in range(nbins): |
---|
[f60a8c2] | 523 | if err_y[n] < 0: |
---|
| 524 | err_y[n] = -err_y[n] |
---|
[095ab1b] | 525 | err_y[n] = math.sqrt(err_y[n]) |
---|
[729bcf6] | 526 | #if err_x != None: |
---|
| 527 | # err_x[n] = math.sqrt(err_x[n]) |
---|
[c8a6c3d7] | 528 | |
---|
[a7a5886] | 529 | err_y = err_y / y_counts |
---|
[9a5097c] | 530 | err_y[err_y == 0] = np.average(err_y) |
---|
[f60a8c2] | 531 | y = y / y_counts |
---|
| 532 | x = x / y_counts |
---|
[9a5097c] | 533 | idx = (np.isfinite(y)) & (np.isfinite(x)) |
---|
[c8a6c3d7] | 534 | |
---|
[342a506] | 535 | if err_x != None: |
---|
| 536 | d_x = err_x[idx] / y_counts[idx] |
---|
| 537 | else: |
---|
| 538 | d_x = None |
---|
| 539 | |
---|
[f60a8c2] | 540 | if not idx.any(): |
---|
| 541 | msg = "Average Error: No points inside ROI to average..." |
---|
[a7a5886] | 542 | raise ValueError, msg |
---|
[c8a6c3d7] | 543 | |
---|
[342a506] | 544 | return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) |
---|
[c8a6c3d7] | 545 | |
---|
[76e2369] | 546 | |
---|
| 547 | class Ring(object): |
---|
| 548 | """ |
---|
[0997158f] | 549 | Defines a ring on a 2D data set. |
---|
| 550 | The ring is defined by r_min, r_max, and |
---|
| 551 | the position of the center of the ring. |
---|
[c8a6c3d7] | 552 | |
---|
[0997158f] | 553 | The data returned is the distribution of counts |
---|
| 554 | around the ring as a function of phi. |
---|
[c8a6c3d7] | 555 | |
---|
[f60a8c2] | 556 | Phi_min and phi_max should be defined between 0 and 2*pi |
---|
[0997158f] | 557 | in anti-clockwise starting from the x- axis on the left-hand side |
---|
[76e2369] | 558 | """ |
---|
[095ab1b] | 559 | #Todo: remove center. |
---|
[400155b] | 560 | def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=36): |
---|
[76e2369] | 561 | # Minimum radius |
---|
| 562 | self.r_min = r_min |
---|
| 563 | # Maximum radius |
---|
| 564 | self.r_max = r_max |
---|
| 565 | # Center of the ring in x |
---|
| 566 | self.center_x = center_x |
---|
| 567 | # Center of the ring in y |
---|
| 568 | self.center_y = center_y |
---|
| 569 | # Number of angular bins |
---|
[8ba103f] | 570 | self.nbins_phi = nbins |
---|
[400155b] | 571 | |
---|
[c8a6c3d7] | 572 | |
---|
[76e2369] | 573 | def __call__(self, data2D): |
---|
| 574 | """ |
---|
[0997158f] | 575 | Apply the ring to the data set. |
---|
| 576 | Returns the angular distribution for a given q range |
---|
[3c3a440] | 577 | |
---|
[0997158f] | 578 | :param data2D: Data2D object |
---|
[3c3a440] | 579 | |
---|
[0997158f] | 580 | :return: Data1D object |
---|
[76e2369] | 581 | """ |
---|
| 582 | if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: |
---|
| 583 | raise RuntimeError, "Ring averaging only take plottable_2D objects" |
---|
[c8a6c3d7] | 584 | |
---|
[095ab1b] | 585 | Pi = math.pi |
---|
[c8a6c3d7] | 586 | |
---|
[095ab1b] | 587 | # Get data |
---|
[9a5097c] | 588 | data = data2D.data[np.isfinite(data2D.data)] |
---|
| 589 | q_data = data2D.q_data[np.isfinite(data2D.data)] |
---|
| 590 | err_data = data2D.err_data[np.isfinite(data2D.data)] |
---|
| 591 | qx_data = data2D.qx_data[np.isfinite(data2D.data)] |
---|
| 592 | qy_data = data2D.qy_data[np.isfinite(data2D.data)] |
---|
[c8a6c3d7] | 593 | |
---|
[095ab1b] | 594 | # Set space for 1d outputs |
---|
[9a5097c] | 595 | phi_bins = np.zeros(self.nbins_phi) |
---|
| 596 | phi_counts = np.zeros(self.nbins_phi) |
---|
| 597 | phi_values = np.zeros(self.nbins_phi) |
---|
| 598 | phi_err = np.zeros(self.nbins_phi) |
---|
[c8a6c3d7] | 599 | |
---|
[3c3a440] | 600 | # Shift to apply to calculated phi values in order |
---|
| 601 | # to center first bin at zero |
---|
[ddc192a] | 602 | phi_shift = Pi / self.nbins_phi |
---|
[400155b] | 603 | |
---|
[f60a8c2] | 604 | for npt in range(len(data)): |
---|
[095ab1b] | 605 | frac = 0 |
---|
| 606 | # q-value at the point (npt) |
---|
| 607 | q_value = q_data[npt] |
---|
[f60a8c2] | 608 | data_n = data[npt] |
---|
[c8a6c3d7] | 609 | |
---|
[095ab1b] | 610 | # phi-value at the point (npt) |
---|
[a7a5886] | 611 | phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi |
---|
[c8a6c3d7] | 612 | |
---|
[a7a5886] | 613 | if self.r_min <= q_value and q_value <= self.r_max: |
---|
[f60a8c2] | 614 | frac = 1 |
---|
[a7a5886] | 615 | if frac == 0: |
---|
| 616 | continue |
---|
[3c3a440] | 617 | # binning |
---|
| 618 | i_phi = int(math.floor((self.nbins_phi) * \ |
---|
| 619 | (phi_value + phi_shift) / (2 * Pi))) |
---|
[c8a6c3d7] | 620 | |
---|
[f60a8c2] | 621 | # Take care of the edge case at phi = 2pi. |
---|
[400155b] | 622 | if i_phi >= self.nbins_phi: |
---|
[c8a6c3d7] | 623 | i_phi = 0 |
---|
[095ab1b] | 624 | phi_bins[i_phi] += frac * data[npt] |
---|
[c8a6c3d7] | 625 | |
---|
[a7a5886] | 626 | if err_data == None or err_data[npt] == 0.0: |
---|
| 627 | if data_n < 0: |
---|
| 628 | data_n = -data_n |
---|
[095ab1b] | 629 | phi_err[i_phi] += frac * frac * math.fabs(data_n) |
---|
| 630 | else: |
---|
[a7a5886] | 631 | phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt] |
---|
[095ab1b] | 632 | phi_counts[i_phi] += frac |
---|
[c8a6c3d7] | 633 | |
---|
[76e2369] | 634 | for i in range(self.nbins_phi): |
---|
| 635 | phi_bins[i] = phi_bins[i] / phi_counts[i] |
---|
| 636 | phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i] |
---|
[400155b] | 637 | phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) |
---|
[c8a6c3d7] | 638 | |
---|
[9a5097c] | 639 | idx = (np.isfinite(phi_bins)) |
---|
[095ab1b] | 640 | |
---|
[a7a5886] | 641 | if not idx.any(): |
---|
[f60a8c2] | 642 | msg = "Average Error: No points inside ROI to average..." |
---|
[a7a5886] | 643 | raise ValueError, msg |
---|
| 644 | #elif len(phi_bins[idx])!= self.nbins_phi: |
---|
| 645 | # print "resulted",self.nbins_phi- len(phi_bins[idx]) |
---|
| 646 | #,"empty bin(s) due to tight binning..." |
---|
[095ab1b] | 647 | return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) |
---|
[c8a6c3d7] | 648 | |
---|
| 649 | |
---|
[76e2369] | 650 | def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11): |
---|
| 651 | """ |
---|
[0997158f] | 652 | Returns the fraction of the pixel defined by |
---|
[f60a8c2] | 653 | the four corners (q_00, q_01, q_10, q_11) that |
---|
[0997158f] | 654 | has q < qmax.:: |
---|
[3c3a440] | 655 | |
---|
[76e2369] | 656 | q_01 q_11 |
---|
| 657 | y=1 +--------------+ |
---|
| 658 | | | |
---|
| 659 | | | |
---|
| 660 | | | |
---|
| 661 | y=0 +--------------+ |
---|
[bb0b12c] | 662 | q_00 q_10 |
---|
[3c3a440] | 663 | |
---|
[76e2369] | 664 | x=0 x=1 |
---|
[3c3a440] | 665 | |
---|
[76e2369] | 666 | """ |
---|
| 667 | # y side for x = minx |
---|
| 668 | x_0 = get_intercept(qmax, q_00, q_01) |
---|
| 669 | # y side for x = maxx |
---|
| 670 | x_1 = get_intercept(qmax, q_10, q_11) |
---|
[c8a6c3d7] | 671 | |
---|
[76e2369] | 672 | # x side for y = miny |
---|
| 673 | y_0 = get_intercept(qmax, q_00, q_10) |
---|
| 674 | # x side for y = maxy |
---|
| 675 | y_1 = get_intercept(qmax, q_01, q_11) |
---|
[c8a6c3d7] | 676 | |
---|
[76e2369] | 677 | # surface fraction for a 1x1 pixel |
---|
| 678 | frac_max = 0 |
---|
[c8a6c3d7] | 679 | |
---|
[76e2369] | 680 | if x_0 and x_1: |
---|
[a7a5886] | 681 | frac_max = (x_0 + x_1) / 2.0 |
---|
[76e2369] | 682 | elif y_0 and y_1: |
---|
[a7a5886] | 683 | frac_max = (y_0 + y_1) / 2.0 |
---|
[76e2369] | 684 | elif x_0 and y_0: |
---|
| 685 | if q_00 < q_10: |
---|
[a7a5886] | 686 | frac_max = x_0 * y_0 / 2.0 |
---|
[76e2369] | 687 | else: |
---|
[a7a5886] | 688 | frac_max = 1.0 - x_0 * y_0 / 2.0 |
---|
[76e2369] | 689 | elif x_0 and y_1: |
---|
| 690 | if q_00 < q_10: |
---|
[a7a5886] | 691 | frac_max = x_0 * y_1 / 2.0 |
---|
[76e2369] | 692 | else: |
---|
[a7a5886] | 693 | frac_max = 1.0 - x_0 * y_1 / 2.0 |
---|
[76e2369] | 694 | elif x_1 and y_0: |
---|
| 695 | if q_00 > q_10: |
---|
[a7a5886] | 696 | frac_max = x_1 * y_0 / 2.0 |
---|
[76e2369] | 697 | else: |
---|
[a7a5886] | 698 | frac_max = 1.0 - x_1 * y_0 / 2.0 |
---|
[76e2369] | 699 | elif x_1 and y_1: |
---|
| 700 | if q_00 < q_10: |
---|
[a7a5886] | 701 | frac_max = 1.0 - (1.0 - x_1) * (1.0 - y_1) / 2.0 |
---|
[76e2369] | 702 | else: |
---|
[a7a5886] | 703 | frac_max = (1.0 - x_1) * (1.0 - y_1) / 2.0 |
---|
[c8a6c3d7] | 704 | |
---|
[76e2369] | 705 | # If we make it here, there is no intercept between |
---|
| 706 | # this pixel and the constant-q ring. We only need |
---|
| 707 | # to know if we have to include it or exclude it. |
---|
[c8a6c3d7] | 708 | elif (q_00 + q_01 + q_10 + q_11) / 4.0 < qmax: |
---|
[76e2369] | 709 | frac_max = 1.0 |
---|
[095ab1b] | 710 | |
---|
[76e2369] | 711 | return frac_max |
---|
[c8a6c3d7] | 712 | |
---|
| 713 | |
---|
[76e2369] | 714 | def get_intercept(q, q_0, q_1): |
---|
| 715 | """ |
---|
[0997158f] | 716 | Returns the fraction of the side at which the |
---|
| 717 | q-value intercept the pixel, None otherwise. |
---|
| 718 | The values returned is the fraction ON THE SIDE |
---|
| 719 | OF THE LOWEST Q. :: |
---|
[3c3a440] | 720 | |
---|
[f60a8c2] | 721 | A B |
---|
[0997158f] | 722 | +-----------+--------+ <--- pixel size |
---|
[f60a8c2] | 723 | 0 1 |
---|
[0997158f] | 724 | Q_0 -------- Q ----- Q_1 <--- equivalent Q range |
---|
[76e2369] | 725 | if Q_1 > Q_0, A is returned |
---|
| 726 | if Q_1 < Q_0, B is returned |
---|
| 727 | if Q is outside the range of [Q_0, Q_1], None is returned |
---|
[3c3a440] | 728 | |
---|
[76e2369] | 729 | """ |
---|
| 730 | if q_1 > q_0: |
---|
[3c3a440] | 731 | if q > q_0 and q <= q_1: |
---|
[f60a8c2] | 732 | return (q - q_0) / (q_1 - q_0) |
---|
[76e2369] | 733 | else: |
---|
[3c3a440] | 734 | if q > q_1 and q <= q_0: |
---|
[f60a8c2] | 735 | return (q - q_1) / (q_0 - q_1) |
---|
[76e2369] | 736 | return None |
---|
[c8a6c3d7] | 737 | |
---|
| 738 | |
---|
[3c3a440] | 739 | class _Sector(object): |
---|
[fb198a9] | 740 | """ |
---|
[0997158f] | 741 | Defines a sector region on a 2D data set. |
---|
| 742 | The sector is defined by r_min, r_max, phi_min, phi_max, |
---|
[f60a8c2] | 743 | and the position of the center of the ring |
---|
[a7a5886] | 744 | where phi_min and phi_max are defined by the right |
---|
| 745 | and left lines wrt central line |
---|
[f60a8c2] | 746 | and phi_max could be less than phi_min. |
---|
[3c3a440] | 747 | |
---|
[f60a8c2] | 748 | Phi is defined between 0 and 2*pi in anti-clockwise |
---|
[a7a5886] | 749 | starting from the x- axis on the left-hand side |
---|
[fb198a9] | 750 | """ |
---|
[c8a6c3d7] | 751 | def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20): |
---|
[fb198a9] | 752 | self.r_min = r_min |
---|
| 753 | self.r_max = r_max |
---|
| 754 | self.phi_min = phi_min |
---|
| 755 | self.phi_max = phi_max |
---|
| 756 | self.nbins = nbins |
---|
[c8a6c3d7] | 757 | |
---|
[fb198a9] | 758 | def _agv(self, data2D, run='phi'): |
---|
| 759 | """ |
---|
[0997158f] | 760 | Perform sector averaging. |
---|
[3c3a440] | 761 | |
---|
[0997158f] | 762 | :param data2D: Data2D object |
---|
| 763 | :param run: define the varying parameter ('phi' , 'q' , or 'q2') |
---|
[3c3a440] | 764 | |
---|
[0997158f] | 765 | :return: Data1D object |
---|
[fb198a9] | 766 | """ |
---|
| 767 | if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: |
---|
| 768 | raise RuntimeError, "Ring averaging only take plottable_2D objects" |
---|
[095ab1b] | 769 | Pi = math.pi |
---|
[c6f95bb] | 770 | |
---|
[095ab1b] | 771 | # Get the all data & info |
---|
[9a5097c] | 772 | data = data2D.data[np.isfinite(data2D.data)] |
---|
| 773 | q_data = data2D.q_data[np.isfinite(data2D.data)] |
---|
| 774 | err_data = data2D.err_data[np.isfinite(data2D.data)] |
---|
| 775 | qx_data = data2D.qx_data[np.isfinite(data2D.data)] |
---|
| 776 | qy_data = data2D.qy_data[np.isfinite(data2D.data)] |
---|
[342a506] | 777 | dq_data = None |
---|
[c8a6c3d7] | 778 | |
---|
[729bcf6] | 779 | # Get the dq for resolution averaging |
---|
[342a506] | 780 | if data2D.dqx_data != None and data2D.dqy_data != None: |
---|
[f60a8c2] | 781 | # The pinholes and det. pix contribution present |
---|
[729bcf6] | 782 | # in both direction of the 2D which must be subtracted when |
---|
| 783 | # converting to 1D: dq_overlap should calculated ideally at |
---|
[f60a8c2] | 784 | # q = 0. |
---|
[729bcf6] | 785 | # Extrapolate dqy(perp) at q = 0 |
---|
| 786 | z_max = max(data2D.q_data) |
---|
| 787 | z_min = min(data2D.q_data) |
---|
| 788 | x_max = data2D.dqx_data[data2D.q_data[z_max]] |
---|
[f60a8c2] | 789 | x_min = data2D.dqx_data[data2D.q_data[z_min]] |
---|
[729bcf6] | 790 | y_max = data2D.dqy_data[data2D.q_data[z_max]] |
---|
[f60a8c2] | 791 | y_min = data2D.dqy_data[data2D.q_data[z_min]] |
---|
[729bcf6] | 792 | # Find qdx at q = 0 |
---|
| 793 | dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) |
---|
| 794 | # when extrapolation goes wrong |
---|
| 795 | if dq_overlap_x > min(data2D.dqx_data): |
---|
| 796 | dq_overlap_x = min(data2D.dqx_data) |
---|
[f60a8c2] | 797 | dq_overlap_x *= dq_overlap_x |
---|
[729bcf6] | 798 | # Find qdx at q = 0 |
---|
| 799 | dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) |
---|
| 800 | # when extrapolation goes wrong |
---|
| 801 | if dq_overlap_y > min(data2D.dqy_data): |
---|
| 802 | dq_overlap_y = min(data2D.dqy_data) |
---|
| 803 | # get dq at q=0. |
---|
| 804 | dq_overlap_y *= dq_overlap_y |
---|
| 805 | |
---|
[9a5097c] | 806 | dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) |
---|
[729bcf6] | 807 | if dq_overlap < 0: |
---|
| 808 | dq_overlap = y_min |
---|
[9a5097c] | 809 | dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] |
---|
| 810 | dqy_data = data2D.dqy_data[np.isfinite(data2D.data)] - dq_overlap |
---|
[729bcf6] | 811 | # def; dqx_data = dq_r dqy_data = dq_phi |
---|
| 812 | # Convert dq 2D to 1D here |
---|
[f60a8c2] | 813 | dqx = dqx_data * dqx_data |
---|
[729bcf6] | 814 | dqy = dqy_data * dqy_data |
---|
[9a5097c] | 815 | dq_data = np.add(dqx, dqy) |
---|
| 816 | dq_data = np.sqrt(dq_data) |
---|
[c8a6c3d7] | 817 | |
---|
[095ab1b] | 818 | #set space for 1d outputs |
---|
[9a5097c] | 819 | x = np.zeros(self.nbins) |
---|
| 820 | y = np.zeros(self.nbins) |
---|
| 821 | y_err = np.zeros(self.nbins) |
---|
| 822 | x_err = np.zeros(self.nbins) |
---|
| 823 | y_counts = np.zeros(self.nbins) |
---|
[c8a6c3d7] | 824 | |
---|
[095ab1b] | 825 | # Get the min and max into the region: 0 <= phi < 2Pi |
---|
| 826 | phi_min = flip_phi(self.phi_min) |
---|
| 827 | phi_max = flip_phi(self.phi_max) |
---|
[c8a6c3d7] | 828 | |
---|
[f60a8c2] | 829 | for n in range(len(data)): |
---|
[a7a5886] | 830 | frac = 0 |
---|
[c8a6c3d7] | 831 | |
---|
[a7a5886] | 832 | # q-value at the pixel (j,i) |
---|
| 833 | q_value = q_data[n] |
---|
| 834 | data_n = data[n] |
---|
[c8a6c3d7] | 835 | |
---|
[a7a5886] | 836 | # Is pixel within range? |
---|
| 837 | is_in = False |
---|
[c8a6c3d7] | 838 | |
---|
[a7a5886] | 839 | # phi-value of the pixel (j,i) |
---|
[f60a8c2] | 840 | phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi |
---|
[c8a6c3d7] | 841 | |
---|
[a7a5886] | 842 | ## No need to calculate the frac when all data are within range |
---|
| 843 | if self.r_min <= q_value and q_value <= self.r_max: |
---|
[f60a8c2] | 844 | frac = 1 |
---|
[a7a5886] | 845 | if frac == 0: |
---|
| 846 | continue |
---|
| 847 | #In case of two ROIs (symmetric major and minor regions)(for 'q2') |
---|
[f60a8c2] | 848 | if run.lower() == 'q2': |
---|
| 849 | ## For minor sector wing |
---|
[a7a5886] | 850 | # Calculate the minor wing phis |
---|
| 851 | phi_min_minor = flip_phi(phi_min - Pi) |
---|
| 852 | phi_max_minor = flip_phi(phi_max - Pi) |
---|
| 853 | # Check if phis of the minor ring is within 0 to 2pi |
---|
| 854 | if phi_min_minor > phi_max_minor: |
---|
| 855 | is_in = (phi_value > phi_min_minor or \ |
---|
| 856 | phi_value < phi_max_minor) |
---|
| 857 | else: |
---|
| 858 | is_in = (phi_value > phi_min_minor and \ |
---|
| 859 | phi_value < phi_max_minor) |
---|
[3c67340] | 860 | |
---|
[f60a8c2] | 861 | #For all cases(i.e.,for 'q', 'q2', and 'phi') |
---|
| 862 | #Find pixels within ROI |
---|
| 863 | if phi_min > phi_max: |
---|
[a7a5886] | 864 | is_in = is_in or (phi_value > phi_min or \ |
---|
[f60a8c2] | 865 | phi_value < phi_max) |
---|
[a7a5886] | 866 | else: |
---|
| 867 | is_in = is_in or (phi_value >= phi_min and \ |
---|
| 868 | phi_value < phi_max) |
---|
[c8a6c3d7] | 869 | |
---|
[a7a5886] | 870 | if not is_in: |
---|
[f60a8c2] | 871 | frac = 0 |
---|
[a7a5886] | 872 | if frac == 0: |
---|
| 873 | continue |
---|
| 874 | # Check which type of averaging we need |
---|
[f60a8c2] | 875 | if run.lower() == 'phi': |
---|
[a7a5886] | 876 | temp_x = (self.nbins) * (phi_value - self.phi_min) |
---|
| 877 | temp_y = (self.phi_max - self.phi_min) |
---|
| 878 | i_bin = int(math.floor(temp_x / temp_y)) |
---|
| 879 | else: |
---|
| 880 | temp_x = (self.nbins) * (q_value - self.r_min) |
---|
[ec3959ab] | 881 | temp_y = (self.r_max - self.r_min) |
---|
[a7a5886] | 882 | i_bin = int(math.floor(temp_x / temp_y)) |
---|
[bb0b12c] | 883 | |
---|
[f60a8c2] | 884 | # Take care of the edge case at phi = 2pi. |
---|
| 885 | if i_bin == self.nbins: |
---|
| 886 | i_bin = self.nbins - 1 |
---|
[c8a6c3d7] | 887 | |
---|
[f60a8c2] | 888 | ## Get the total y |
---|
[a7a5886] | 889 | y[i_bin] += frac * data_n |
---|
[729bcf6] | 890 | x[i_bin] += frac * q_value |
---|
[342a506] | 891 | if err_data[n] == None or err_data[n] == 0.0: |
---|
[a7a5886] | 892 | if data_n < 0: |
---|
| 893 | data_n = -data_n |
---|
| 894 | y_err[i_bin] += frac * frac * data_n |
---|
| 895 | else: |
---|
| 896 | y_err[i_bin] += frac * frac * err_data[n] * err_data[n] |
---|
[c8a6c3d7] | 897 | |
---|
[342a506] | 898 | if dq_data != None: |
---|
[f60a8c2] | 899 | # To be consistent with dq calculation in 1d reduction, |
---|
| 900 | # we need just the averages (not quadratures) because |
---|
| 901 | # it should not depend on the number of the q points |
---|
[729bcf6] | 902 | # in the qr bins. |
---|
| 903 | x_err[i_bin] += frac * dq_data[n] |
---|
[342a506] | 904 | else: |
---|
| 905 | x_err = None |
---|
[a7a5886] | 906 | y_counts[i_bin] += frac |
---|
[c8a6c3d7] | 907 | |
---|
[095ab1b] | 908 | # Organize the results |
---|
[fb198a9] | 909 | for i in range(self.nbins): |
---|
| 910 | y[i] = y[i] / y_counts[i] |
---|
| 911 | y_err[i] = math.sqrt(y_err[i]) / y_counts[i] |
---|
[729bcf6] | 912 | |
---|
[095ab1b] | 913 | # The type of averaging: phi,q2, or q |
---|
| 914 | # Calculate x[i]should be at the center of the bin |
---|
[f60a8c2] | 915 | if run.lower() == 'phi': |
---|
[12c5b87] | 916 | x[i] = (self.phi_max - self.phi_min) / self.nbins * \ |
---|
| 917 | (1.0 * i + 0.5) + self.phi_min |
---|
[095ab1b] | 918 | else: |
---|
[f60a8c2] | 919 | # We take the center of ring area, not radius. |
---|
[342a506] | 920 | # This is more accurate than taking the radial center of ring. |
---|
[729bcf6] | 921 | #delta_r = (self.r_max - self.r_min) / self.nbins |
---|
| 922 | #r_inner = self.r_min + delta_r * i |
---|
| 923 | #r_outer = r_inner + delta_r |
---|
| 924 | #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) |
---|
| 925 | x[i] = x[i] / y_counts[i] |
---|
[9a5097c] | 926 | y_err[y_err == 0] = np.average(y_err) |
---|
| 927 | idx = (np.isfinite(y) & np.isfinite(y_err)) |
---|
[342a506] | 928 | if x_err != None: |
---|
[729bcf6] | 929 | d_x = x_err[idx] / y_counts[idx] |
---|
[342a506] | 930 | else: |
---|
| 931 | d_x = None |
---|
[a7a5886] | 932 | if not idx.any(): |
---|
[f60a8c2] | 933 | msg = "Average Error: No points inside sector of ROI to average..." |
---|
[a7a5886] | 934 | raise ValueError, msg |
---|
| 935 | #elif len(y[idx])!= self.nbins: |
---|
| 936 | # print "resulted",self.nbins- len(y[idx]), |
---|
| 937 | #"empty bin(s) due to tight binning..." |
---|
[342a506] | 938 | return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x) |
---|
[c8a6c3d7] | 939 | |
---|
| 940 | |
---|
[2e83ff3] | 941 | class SectorPhi(_Sector): |
---|
| 942 | """ |
---|
[0997158f] | 943 | Sector average as a function of phi. |
---|
| 944 | I(phi) is return and the data is averaged over Q. |
---|
[3c3a440] | 945 | |
---|
[0997158f] | 946 | A sector is defined by r_min, r_max, phi_min, phi_max. |
---|
| 947 | The number of bin in phi also has to be defined. |
---|
[2e83ff3] | 948 | """ |
---|
| 949 | def __call__(self, data2D): |
---|
| 950 | """ |
---|
[0997158f] | 951 | Perform sector average and return I(phi). |
---|
[3c3a440] | 952 | |
---|
[0997158f] | 953 | :param data2D: Data2D object |
---|
| 954 | :return: Data1D object |
---|
[2e83ff3] | 955 | """ |
---|
| 956 | return self._agv(data2D, 'phi') |
---|
[c8a6c3d7] | 957 | |
---|
| 958 | |
---|
[fb198a9] | 959 | class SectorQ(_Sector): |
---|
| 960 | """ |
---|
[0997158f] | 961 | Sector average as a function of Q for both symatric wings. |
---|
| 962 | I(Q) is return and the data is averaged over phi. |
---|
[3c3a440] | 963 | |
---|
[0997158f] | 964 | A sector is defined by r_min, r_max, phi_min, phi_max. |
---|
[f60a8c2] | 965 | r_min, r_max, phi_min, phi_max >0. |
---|
[0997158f] | 966 | The number of bin in Q also has to be defined. |
---|
[fb198a9] | 967 | """ |
---|
| 968 | def __call__(self, data2D): |
---|
| 969 | """ |
---|
[0997158f] | 970 | Perform sector average and return I(Q). |
---|
[3c3a440] | 971 | |
---|
[0997158f] | 972 | :param data2D: Data2D object |
---|
[3c3a440] | 973 | |
---|
[0997158f] | 974 | :return: Data1D object |
---|
[fb198a9] | 975 | """ |
---|
| 976 | return self._agv(data2D, 'q2') |
---|
[c6f95bb] | 977 | |
---|
[f60a8c2] | 978 | |
---|
[f265927] | 979 | class Ringcut(object): |
---|
| 980 | """ |
---|
[0997158f] | 981 | Defines a ring on a 2D data set. |
---|
| 982 | The ring is defined by r_min, r_max, and |
---|
| 983 | the position of the center of the ring. |
---|
[3c3a440] | 984 | |
---|
[0997158f] | 985 | The data returned is the region inside the ring |
---|
[3c3a440] | 986 | |
---|
[f60a8c2] | 987 | Phi_min and phi_max should be defined between 0 and 2*pi |
---|
[0997158f] | 988 | in anti-clockwise starting from the x- axis on the left-hand side |
---|
[f265927] | 989 | """ |
---|
[f60a8c2] | 990 | def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0): |
---|
[f265927] | 991 | # Minimum radius |
---|
| 992 | self.r_min = r_min |
---|
| 993 | # Maximum radius |
---|
| 994 | self.r_max = r_max |
---|
| 995 | # Center of the ring in x |
---|
| 996 | self.center_x = center_x |
---|
| 997 | # Center of the ring in y |
---|
| 998 | self.center_y = center_y |
---|
| 999 | |
---|
| 1000 | def __call__(self, data2D): |
---|
| 1001 | """ |
---|
[0997158f] | 1002 | Apply the ring to the data set. |
---|
| 1003 | Returns the angular distribution for a given q range |
---|
[3c3a440] | 1004 | |
---|
[0997158f] | 1005 | :param data2D: Data2D object |
---|
[3c3a440] | 1006 | |
---|
[0997158f] | 1007 | :return: index array in the range |
---|
[f265927] | 1008 | """ |
---|
| 1009 | if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: |
---|
| 1010 | raise RuntimeError, "Ring cut only take plottable_2D objects" |
---|
| 1011 | |
---|
| 1012 | # Get data |
---|
[f60a8c2] | 1013 | qx_data = data2D.qx_data |
---|
[f265927] | 1014 | qy_data = data2D.qy_data |
---|
[9a5097c] | 1015 | q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) |
---|
[f265927] | 1016 | |
---|
| 1017 | # check whether or not the data point is inside ROI |
---|
| 1018 | out = (self.r_min <= q_data) & (self.r_max >= q_data) |
---|
[3c3a440] | 1019 | return out |
---|
[c8a6c3d7] | 1020 | |
---|
[f265927] | 1021 | |
---|
[c6f95bb] | 1022 | class Boxcut(object): |
---|
| 1023 | """ |
---|
[0997158f] | 1024 | Find a rectangular 2D region of interest. |
---|
[c6f95bb] | 1025 | """ |
---|
| 1026 | def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): |
---|
| 1027 | # Minimum Qx value [A-1] |
---|
| 1028 | self.x_min = x_min |
---|
| 1029 | # Maximum Qx value [A-1] |
---|
| 1030 | self.x_max = x_max |
---|
| 1031 | # Minimum Qy value [A-1] |
---|
| 1032 | self.y_min = y_min |
---|
| 1033 | # Maximum Qy value [A-1] |
---|
| 1034 | self.y_max = y_max |
---|
| 1035 | |
---|
| 1036 | def __call__(self, data2D): |
---|
| 1037 | """ |
---|
[0997158f] | 1038 | Find a rectangular 2D region of interest. |
---|
[3c3a440] | 1039 | |
---|
[0997158f] | 1040 | :param data2D: Data2D object |
---|
[f60a8c2] | 1041 | :return: mask, 1d array (len = len(data)) |
---|
[0997158f] | 1042 | with Trues where the data points are inside ROI, otherwise False |
---|
[c6f95bb] | 1043 | """ |
---|
| 1044 | mask = self._find(data2D) |
---|
[c8a6c3d7] | 1045 | |
---|
[c6f95bb] | 1046 | return mask |
---|
[c8a6c3d7] | 1047 | |
---|
[c6f95bb] | 1048 | def _find(self, data2D): |
---|
| 1049 | """ |
---|
[f60a8c2] | 1050 | Find a rectangular 2D region of interest. |
---|
[3c3a440] | 1051 | |
---|
[0997158f] | 1052 | :param data2D: Data2D object |
---|
[3c3a440] | 1053 | |
---|
[f60a8c2] | 1054 | :return: out, 1d array (length = len(data)) |
---|
[0997158f] | 1055 | with Trues where the data points are inside ROI, otherwise Falses |
---|
[c6f95bb] | 1056 | """ |
---|
| 1057 | if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: |
---|
| 1058 | raise RuntimeError, "Boxcut take only plottable_2D objects" |
---|
[f60a8c2] | 1059 | # Get qx_ and qy_data |
---|
[c6f95bb] | 1060 | qx_data = data2D.qx_data |
---|
| 1061 | qy_data = data2D.qy_data |
---|
[c8a6c3d7] | 1062 | |
---|
[c6f95bb] | 1063 | # check whether or not the data point is inside ROI |
---|
[f265927] | 1064 | outx = (self.x_min <= qx_data) & (self.x_max > qx_data) |
---|
| 1065 | outy = (self.y_min <= qy_data) & (self.y_max > qy_data) |
---|
[c6f95bb] | 1066 | |
---|
[3c3a440] | 1067 | return outx & outy |
---|
[c6f95bb] | 1068 | |
---|
[f60a8c2] | 1069 | |
---|
[c6f95bb] | 1070 | class Sectorcut(object): |
---|
| 1071 | """ |
---|
[0997158f] | 1072 | Defines a sector (major + minor) region on a 2D data set. |
---|
| 1073 | The sector is defined by phi_min, phi_max, |
---|
[f60a8c2] | 1074 | where phi_min and phi_max are defined by the right |
---|
| 1075 | and left lines wrt central line. |
---|
[3c3a440] | 1076 | |
---|
[f60a8c2] | 1077 | Phi_min and phi_max are given in units of radian |
---|
[0997158f] | 1078 | and (phi_max-phi_min) should not be larger than pi |
---|
[c6f95bb] | 1079 | """ |
---|
[a7a5886] | 1080 | def __init__(self, phi_min=0, phi_max=math.pi): |
---|
[c6f95bb] | 1081 | self.phi_min = phi_min |
---|
| 1082 | self.phi_max = phi_max |
---|
[c8a6c3d7] | 1083 | |
---|
[c6f95bb] | 1084 | def __call__(self, data2D): |
---|
| 1085 | """ |
---|
[0997158f] | 1086 | Find a rectangular 2D region of interest. |
---|
[3c3a440] | 1087 | |
---|
[0997158f] | 1088 | :param data2D: Data2D object |
---|
[3c3a440] | 1089 | |
---|
[f60a8c2] | 1090 | :return: mask, 1d array (len = len(data)) |
---|
[3c3a440] | 1091 | |
---|
[0997158f] | 1092 | with Trues where the data points are inside ROI, otherwise False |
---|
[c6f95bb] | 1093 | """ |
---|
| 1094 | mask = self._find(data2D) |
---|
[c8a6c3d7] | 1095 | |
---|
[c6f95bb] | 1096 | return mask |
---|
[c8a6c3d7] | 1097 | |
---|
[c6f95bb] | 1098 | def _find(self, data2D): |
---|
| 1099 | """ |
---|
[f60a8c2] | 1100 | Find a rectangular 2D region of interest. |
---|
[3c3a440] | 1101 | |
---|
[0997158f] | 1102 | :param data2D: Data2D object |
---|
[3c3a440] | 1103 | |
---|
[f60a8c2] | 1104 | :return: out, 1d array (length = len(data)) |
---|
[3c3a440] | 1105 | |
---|
[0997158f] | 1106 | with Trues where the data points are inside ROI, otherwise Falses |
---|
[c6f95bb] | 1107 | """ |
---|
| 1108 | if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: |
---|
[f60a8c2] | 1109 | raise RuntimeError, "Sectorcut take only plottable_2D objects" |
---|
[c6f95bb] | 1110 | Pi = math.pi |
---|
[3c3a440] | 1111 | # Get data |
---|
[c6f95bb] | 1112 | qx_data = data2D.qx_data |
---|
[f60a8c2] | 1113 | qy_data = data2D.qy_data |
---|
[c6f95bb] | 1114 | |
---|
| 1115 | # get phi from data |
---|
[9a5097c] | 1116 | phi_data = np.arctan2(qy_data, qx_data) |
---|
[c8a6c3d7] | 1117 | |
---|
[f265927] | 1118 | # Get the min and max into the region: -pi <= phi < Pi |
---|
[a7a5886] | 1119 | phi_min_major = flip_phi(self.phi_min + Pi) - Pi |
---|
[f60a8c2] | 1120 | phi_max_major = flip_phi(self.phi_max + Pi) - Pi |
---|
[c6f95bb] | 1121 | # check for major sector |
---|
[f265927] | 1122 | if phi_min_major > phi_max_major: |
---|
| 1123 | out_major = (phi_min_major <= phi_data) + (phi_max_major > phi_data) |
---|
[c6f95bb] | 1124 | else: |
---|
[f265927] | 1125 | out_major = (phi_min_major <= phi_data) & (phi_max_major > phi_data) |
---|
[c8a6c3d7] | 1126 | |
---|
[c6f95bb] | 1127 | # minor sector |
---|
| 1128 | # Get the min and max into the region: -pi <= phi < Pi |
---|
[a7a5886] | 1129 | phi_min_minor = flip_phi(self.phi_min) - Pi |
---|
| 1130 | phi_max_minor = flip_phi(self.phi_max) - Pi |
---|
[c8a6c3d7] | 1131 | |
---|
[c6f95bb] | 1132 | # check for minor sector |
---|
| 1133 | if phi_min_minor > phi_max_minor: |
---|
[a7a5886] | 1134 | out_minor = (phi_min_minor <= phi_data) + \ |
---|
[f60a8c2] | 1135 | (phi_max_minor >= phi_data) |
---|
[c6f95bb] | 1136 | else: |
---|
[a7a5886] | 1137 | out_minor = (phi_min_minor <= phi_data) & \ |
---|
[f60a8c2] | 1138 | (phi_max_minor >= phi_data) |
---|
[c6f95bb] | 1139 | out = out_major + out_minor |
---|
[c8a6c3d7] | 1140 | |
---|
[c6f95bb] | 1141 | return out |
---|