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