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