- Timestamp:
- Mar 3, 2015 8:28:02 AM (10 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 3c3a440
- Parents:
- f06d7fc
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/dataloader/manipulations.py
rddc192a rc8a6c3d7 24 24 :param dx: x-distance from beam center [mm] 25 25 :param dy: y-distance from beam center [mm] 26 27 26 :return: q-value at the given position 28 27 """ 29 28 # Distance from beam center in the plane of detector 30 plane_dist = math.sqrt(dx *dx + dy*dy)29 plane_dist = math.sqrt(dx * dx + dy * dy) 31 30 # Half of the scattering angle 32 theta = 0.5 * math.atan(plane_dist /det_dist)33 return (4.0 * math.pi /wavelength) * math.sin(theta)31 theta = 0.5 * math.atan(plane_dist / det_dist) 32 return (4.0 * math.pi / wavelength) * math.sin(theta) 34 33 35 34 … … 45 44 angle_xy = math.pi 46 45 else: 47 angle_xy = math.atan(dx /dy)48 46 angle_xy = math.atan(dx / dy) 47 49 48 if compo == "x": 50 49 out = get_q(dx, dy, det_dist, wavelength) * math.cos(angle_xy) … … 59 58 """ 60 59 Correct phi to within the 0 <= to <= 2pi range 61 60 62 61 :return: phi in >=0 and <=2Pi 63 62 """ … … 76 75 convert old 2d format opened by IhorReader or danse_reader 77 76 to new Data2D format 78 77 79 78 :param data2d: 2d array of Data2D object 80 81 79 :return: 1d arrays of Data2D object 82 80 83 81 """ 84 82 if data2d.data == None or data2d.x_bins == None or data2d.y_bins == None: 85 83 raise ValueError, "Can't convert this data: data=None..." 86 87 from sas.dataloader.data_info import Data2D88 89 84 new_x = numpy.tile(data2d.x_bins, (len(data2d.y_bins), 1)) 90 85 new_y = numpy.tile(data2d.y_bins, (len(data2d.x_bins), 1)) … … 94 89 qx_data = new_x.flatten() 95 90 qy_data = new_y.flatten() 96 q_data = numpy.sqrt(qx_data *qx_data + qy_data*qy_data)91 q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data) 97 92 if data2d.err_data == None or numpy.any(data2d.err_data <= 0): 98 93 new_err_data = numpy.sqrt(numpy.abs(new_data)) … … 102 97 103 98 #TODO: make sense of the following two lines... 104 output = Data2D() 99 #from sas.dataloader.data_info import Data2D 100 #output = Data2D() 105 101 output = data2d 106 102 output.data = new_data … … 133 129 # negative q-values are allowed 134 130 self.fold = False 135 131 136 132 def __call__(self, data2D): 137 133 return NotImplemented 138 134 139 135 def _avg(self, data2D, maj): 140 136 """ … … 142 138 The major axis is defined as the axis of Q_maj. 143 139 The minor axis is the axis that we average over. 144 140 145 141 :param data2D: Data2D object 146 142 :param maj_min: min value on the major axis 147 148 143 :return: Data1D object 149 144 """ … … 152 147 msg += " detectors: %g" % len(data2D.detector) 153 148 raise RuntimeError, msg 154 149 155 150 # Get data 156 151 data = data2D.data[numpy.isfinite(data2D.data)] 157 q_data = data2D.q_data[numpy.isfinite(data2D.data)]158 152 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 159 153 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 160 154 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 161 155 162 156 # Build array of Q intervals 163 157 if maj == 'x': … … 167 161 x_min = self.x_min 168 162 nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) 169 qbins = self.bin_width * numpy.arange(nbins) + x_min170 163 elif maj == 'y': 171 164 if self.fold: … … 173 166 else: 174 167 y_min = self.y_min 175 nbins = int(math.ceil((self.y_max - y_min)/self.bin_width)) 176 qbins = self.bin_width * numpy.arange(nbins) + y_min 168 nbins = int(math.ceil((self.y_max - y_min) / self.bin_width)) 177 169 else: 178 170 raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj) 179 171 180 172 x = numpy.zeros(nbins) 181 173 y = numpy.zeros(nbins) … … 194 186 frac_y = 1 195 187 frac = frac_x * frac_y 196 188 197 189 if frac == 0: 198 190 continue … … 200 192 if maj == 'x': 201 193 q_value = qx_data[npts] 202 min = x_min194 min_value = x_min 203 195 if maj == 'y': 204 196 q_value = qy_data[npts] 205 min = y_min197 min_value = y_min 206 198 if self.fold and q_value < 0: 207 199 q_value = -q_value 208 200 # bin 209 i_q = int(math.ceil((q_value - min ) / self.bin_width)) - 1210 201 i_q = int(math.ceil((q_value - min_value) / self.bin_width)) - 1 202 211 203 # skip outside of max bins 212 204 if i_q < 0 or i_q >= nbins: 213 205 continue 214 206 215 207 #TODO: find better definition of x[i_q] based on q_data 216 x[i_q] += frac * q_value # min + (i_q + 1) * self.bin_width / 2.0208 x[i_q] += frac * q_value # min_value + (i_q + 1) * self.bin_width / 2.0 217 209 y[i_q] += frac * data[npts] 218 210 219 211 if err_data == None or err_data[npts] == 0.0: 220 212 if data[npts] < 0: … … 224 216 err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] 225 217 y_counts[i_q] += frac 226 218 227 219 # Average the sums 228 220 for n in range(nbins): 229 221 err_y[n] = math.sqrt(err_y[n]) 230 222 231 223 err_y = err_y / y_counts 232 224 y = y / y_counts 233 225 x = x / y_counts 234 226 idx = (numpy.isfinite(y) & numpy.isfinite(x)) 235 236 if not idx.any(): 227 228 if not idx.any(): 237 229 msg = "Average Error: No points inside ROI to average..." 238 230 raise ValueError, msg 239 #elif len(y[idx])!= nbins:240 # msg = "empty bin(s) due to tight binning..."241 # print "resulted",nbins- len(y[idx]), msg242 231 return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) 243 244 232 233 245 234 class SlabY(_Slab): 246 235 """ … … 250 239 """ 251 240 Compute average I(Qy) for a region of interest 252 253 :param data2D: Data2D object 254 241 242 :param data2D: Data2D object 255 243 :return: Data1D object 256 244 """ 257 245 return self._avg(data2D, 'y') 258 259 246 247 260 248 class SlabX(_Slab): 261 249 """ … … 265 253 """ 266 254 Compute average I(Qx) for a region of interest 267 268 :param data2D: Data2D object 269 255 :param data2D: Data2D object 270 256 :return: Data1D object 271 272 257 """ 273 258 return self._avg(data2D, 'x') … … 291 276 """ 292 277 Perform the sum in the region of interest 293 294 :param data2D: Data2D object 295 278 279 :param data2D: Data2D object 296 280 :return: number of counts, error on number of counts, 297 281 number of points summed 298 299 282 """ 300 283 y, err_y, y_counts = self._sum(data2D) 301 284 302 285 # Average the sums 303 286 counts = 0 if y_counts == 0 else y 304 287 error = 0 if y_counts == 0 else math.sqrt(err_y) 305 288 306 289 # Added y_counts to return, SMK & PDB, 04/03/2013 307 290 return counts, error, y_counts 308 291 309 292 def _sum(self, data2D): 310 293 """ 311 294 Perform the sum in the region of interest 312 313 :param data2D: Data2D object 314 295 296 :param data2D: Data2D object 315 297 :return: number of counts, 316 298 error on number of counts, number of entries summed 317 318 299 """ 319 300 if len(data2D.detector) != 1: … … 323 304 # Get data 324 305 data = data2D.data[numpy.isfinite(data2D.data)] 325 q_data = data2D.q_data[numpy.isfinite(data2D.data)]326 306 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 327 307 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 328 308 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 329 309 330 310 y = 0.0 331 311 err_y = 0.0 … … 337 317 frac_x = 0 338 318 frac_y = 0 339 319 340 320 # get min and max at each points 341 321 qx = qx_data[npts] 342 322 qy = qy_data[npts] 343 323 344 324 # get the ROI 345 325 if self.x_min <= qx and self.x_max > qx: … … 368 348 def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 369 349 super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, 370 350 y_min=y_min, y_max=y_max) 371 351 372 352 def __call__(self, data2D): 373 353 """ 374 354 Perform the sum in the region of interest 375 376 :param data2D: Data2D object 377 355 356 :param data2D: Data2D object 378 357 :return: average counts, error on average counts 379 358 380 359 """ 381 360 y, err_y, y_counts = self._sum(data2D) 382 361 383 362 # Average the sums 384 363 counts = 0 if y_counts == 0 else y / y_counts 385 364 error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts 386 365 387 366 return counts, error 388 389 367 368 390 369 def get_pixel_fraction_square(x, xmin, xmax): 391 370 """ 392 371 Return the fraction of the length 393 372 from xmin to x.:: 394 395 373 396 374 A B 397 375 +-----------+---------+ 398 376 xmin x xmax 399 377 400 378 :param x: x-value 401 379 :param xmin: minimum x for the length considered 402 380 :param xmax: minimum x for the length considered 403 404 381 :return: (x-xmin)/(xmax-xmin) when xmin < x < xmax 405 382 406 383 """ 407 384 if x <= xmin: … … 416 393 """ 417 394 Perform circular averaging on 2D data 418 395 419 396 The data returned is the distribution of counts 420 397 as a function of Q … … 431 408 """ 432 409 Perform circular averaging on the data 433 434 :param data2D: Data2D object 435 410 411 :param data2D: Data2D object 436 412 :return: Data1D object 437 413 """ … … 439 415 data = data2D.data[numpy.isfinite(data2D.data)] 440 416 q_data = data2D.q_data[numpy.isfinite(data2D.data)] 441 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)]442 417 err_data = data2D.err_data[numpy.isfinite(data2D.data)] 443 418 mask_data = data2D.mask[numpy.isfinite(data2D.data)] 444 419 445 420 dq_data = None 446 421 447 422 # Get the dq for resolution averaging 448 423 if data2D.dqx_data != None and data2D.dqy_data != None: … … 484 459 dq_data = numpy.add(dqx, dqy) 485 460 dq_data = numpy.sqrt(dq_data) 486 461 487 462 #q_data_max = numpy.max(q_data) 488 463 if len(data2D.q_data) == None: … … 492 467 # Build array of Q intervals 493 468 nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) 494 qbins = self.bin_width * numpy.arange(nbins) + self.r_min495 469 496 470 x = numpy.zeros(nbins) … … 501 475 502 476 for npt in range(len(data)): 503 477 504 478 if ismask and not mask_data[npt]: 505 479 continue 506 480 507 481 frac = 0 508 482 509 483 # q-value at the pixel (j,i) 510 484 q_value = q_data[npt] 511 485 data_n = data[npt] 512 486 513 487 ## No need to calculate the frac when all data are within range 514 488 if self.r_min >= self.r_max: 515 489 raise ValueError, "Limit Error: min > max" 516 490 517 491 if self.r_min <= q_value and q_value <= self.r_max: 518 492 frac = 1 519 493 if frac == 0: 520 continue 494 continue 521 495 i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 522 496 … … 542 516 err_x = None 543 517 y_counts[i_q] += frac 544 518 545 519 # Average the sums 546 520 for n in range(nbins): … … 550 524 #if err_x != None: 551 525 # err_x[n] = math.sqrt(err_x[n]) 552 526 553 527 err_y = err_y / y_counts 554 528 err_y[err_y == 0] = numpy.average(err_y) … … 556 530 x = x / y_counts 557 531 idx = (numpy.isfinite(y)) & (numpy.isfinite(x)) 558 532 559 533 if err_x != None: 560 534 d_x = err_x[idx] / y_counts[idx] … … 565 539 msg = "Average Error: No points inside ROI to average..." 566 540 raise ValueError, msg 567 541 568 542 return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) 569 543 570 544 571 545 class Ring(object): … … 574 548 The ring is defined by r_min, r_max, and 575 549 the position of the center of the ring. 576 550 577 551 The data returned is the distribution of counts 578 552 around the ring as a function of phi. 579 553 580 554 Phi_min and phi_max should be defined between 0 and 2*pi 581 555 in anti-clockwise starting from the x- axis on the left-hand side … … 594 568 self.nbins_phi = nbins 595 569 596 570 597 571 def __call__(self, data2D): 598 572 """ … … 606 580 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 607 581 raise RuntimeError, "Ring averaging only take plottable_2D objects" 608 582 609 583 Pi = math.pi 610 584 611 585 # Get data 612 586 data = data2D.data[numpy.isfinite(data2D.data)] … … 615 589 qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 616 590 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 617 618 q_data_max = numpy.max(q_data) 619 591 620 592 # Set space for 1d outputs 621 phi_bins 593 phi_bins = numpy.zeros(self.nbins_phi) 622 594 phi_counts = numpy.zeros(self.nbins_phi) 623 595 phi_values = numpy.zeros(self.nbins_phi) 624 phi_err 625 596 phi_err = numpy.zeros(self.nbins_phi) 597 626 598 # Shift to apply to calculated phi values in order to center first bin at zero 627 599 phi_shift = Pi / self.nbins_phi … … 632 604 q_value = q_data[npt] 633 605 data_n = data[npt] 634 606 635 607 # phi-value at the point (npt) 636 608 phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi 637 609 638 610 if self.r_min <= q_value and q_value <= self.r_max: 639 611 frac = 1 … … 641 613 continue 642 614 # binning 643 i_phi = int(math.floor((self.nbins_phi) * (phi_value +phi_shift) / (2 * Pi)))644 615 i_phi = int(math.floor((self.nbins_phi) * (phi_value + phi_shift) / (2 * Pi))) 616 645 617 # Take care of the edge case at phi = 2pi. 646 618 if i_phi >= self.nbins_phi: 647 i_phi = 619 i_phi = 0 648 620 phi_bins[i_phi] += frac * data[npt] 649 621 650 622 if err_data == None or err_data[npt] == 0.0: 651 623 if data_n < 0: … … 655 627 phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt] 656 628 phi_counts[i_phi] += frac 657 629 658 630 for i in range(self.nbins_phi): 659 631 phi_bins[i] = phi_bins[i] / phi_counts[i] 660 632 phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i] 661 633 phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) 662 634 663 635 idx = (numpy.isfinite(phi_bins)) 664 636 … … 670 642 #,"empty bin(s) due to tight binning..." 671 643 return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) 672 673 644 645 674 646 def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11): 675 647 """ … … 693 665 # y side for x = maxx 694 666 x_1 = get_intercept(qmax, q_10, q_11) 695 667 696 668 # x side for y = miny 697 669 y_0 = get_intercept(qmax, q_00, q_10) 698 670 # x side for y = maxy 699 671 y_1 = get_intercept(qmax, q_01, q_11) 700 672 701 673 # surface fraction for a 1x1 pixel 702 674 frac_max = 0 703 675 704 676 if x_0 and x_1: 705 677 frac_max = (x_0 + x_1) / 2.0 … … 726 698 else: 727 699 frac_max = (1.0 - x_1) * (1.0 - y_1) / 2.0 728 700 729 701 # If we make it here, there is no intercept between 730 702 # this pixel and the constant-q ring. We only need 731 703 # to know if we have to include it or exclude it. 732 elif (q_00 + q_01 + q_10 + q_11) / 4.0 < 704 elif (q_00 + q_01 + q_10 + q_11) / 4.0 < qmax: 733 705 frac_max = 1.0 734 706 735 707 return frac_max 736 737 708 709 738 710 def get_intercept(q, q_0, q_1): 739 711 """ … … 760 732 return (q - q_1) / (q_0 - q_1) 761 733 return None 762 763 734 735 764 736 class _Sector: 765 737 """ … … 774 746 starting from the x- axis on the left-hand side 775 747 """ 776 def __init__(self, r_min, r_max, phi_min=0, phi_max=2 *math.pi, nbins=20):748 def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20): 777 749 self.r_min = r_min 778 750 self.r_max = r_max … … 780 752 self.phi_max = phi_max 781 753 self.nbins = nbins 782 754 783 755 def _agv(self, data2D, run='phi'): 784 756 """ … … 801 773 qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 802 774 dq_data = None 803 775 804 776 # Get the dq for resolution averaging 805 777 if data2D.dqx_data != None and data2D.dqy_data != None: … … 840 812 dq_data = numpy.add(dqx, dqy) 841 813 dq_data = numpy.sqrt(dq_data) 842 814 843 815 #set space for 1d outputs 844 x 845 y 846 y_err 847 x_err 816 x = numpy.zeros(self.nbins) 817 y = numpy.zeros(self.nbins) 818 y_err = numpy.zeros(self.nbins) 819 x_err = numpy.zeros(self.nbins) 848 820 y_counts = numpy.zeros(self.nbins) 849 821 850 822 # Get the min and max into the region: 0 <= phi < 2Pi 851 823 phi_min = flip_phi(self.phi_min) 852 824 phi_max = flip_phi(self.phi_max) 853 854 q_data_max = numpy.max(q_data) 855 825 856 826 for n in range(len(data)): 857 827 frac = 0 858 828 859 829 # q-value at the pixel (j,i) 860 830 q_value = q_data[n] 861 831 data_n = data[n] 862 832 863 833 # Is pixel within range? 864 834 is_in = False 865 835 866 836 # phi-value of the pixel (j,i) 867 837 phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 868 838 869 839 ## No need to calculate the frac when all data are within range 870 840 if self.r_min <= q_value and q_value <= self.r_max: … … 894 864 is_in = is_in or (phi_value >= phi_min and \ 895 865 phi_value < phi_max) 896 866 897 867 if not is_in: 898 868 frac = 0 … … 912 882 if i_bin == self.nbins: 913 883 i_bin = self.nbins - 1 914 884 915 885 ## Get the total y 916 886 y[i_bin] += frac * data_n … … 922 892 else: 923 893 y_err[i_bin] += frac * frac * err_data[n] * err_data[n] 924 894 925 895 if dq_data != None: 926 896 # To be consistent with dq calculation in 1d reduction, … … 932 902 x_err = None 933 903 y_counts[i_bin] += frac 934 904 935 905 # Organize the results 936 906 for i in range(self.nbins): … … 964 934 #"empty bin(s) due to tight binning..." 965 935 return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x) 966 967 936 937 968 938 class SectorPhi(_Sector): 969 939 """ … … 982 952 """ 983 953 return self._agv(data2D, 'phi') 984 985 954 955 986 956 class SectorQ(_Sector): 987 957 """ … … 1040 1010 qx_data = data2D.qx_data 1041 1011 qy_data = data2D.qy_data 1042 mask = data2D.mask1043 1012 q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data) 1044 1013 … … 1047 1016 1048 1017 return (out) 1049 1018 1050 1019 1051 1020 class Boxcut(object): … … 1072 1041 """ 1073 1042 mask = self._find(data2D) 1074 1043 1075 1044 return mask 1076 1045 1077 1046 def _find(self, data2D): 1078 1047 """ … … 1089 1058 qx_data = data2D.qx_data 1090 1059 qy_data = data2D.qy_data 1091 mask = data2D.mask 1092 1060 1093 1061 # check whether or not the data point is inside ROI 1094 1062 outx = (self.x_min <= qx_data) & (self.x_max > qx_data) … … 1111 1079 self.phi_min = phi_min 1112 1080 self.phi_max = phi_max 1113 1081 1114 1082 def __call__(self, data2D): 1115 1083 """ … … 1123 1091 """ 1124 1092 mask = self._find(data2D) 1125 1093 1126 1094 return mask 1127 1095 1128 1096 def _find(self, data2D): 1129 1097 """ … … 1142 1110 qx_data = data2D.qx_data 1143 1111 qy_data = data2D.qy_data 1144 phi_data = numpy.zeros(len(qx_data))1145 1112 1146 1113 # get phi from data 1147 1114 phi_data = numpy.arctan2(qy_data, qx_data) 1148 1115 1149 1116 # Get the min and max into the region: -pi <= phi < Pi 1150 1117 phi_min_major = flip_phi(self.phi_min + Pi) - Pi … … 1155 1122 else: 1156 1123 out_major = (phi_min_major <= phi_data) & (phi_max_major > phi_data) 1157 1124 1158 1125 # minor sector 1159 1126 # Get the min and max into the region: -pi <= phi < Pi 1160 1127 phi_min_minor = flip_phi(self.phi_min) - Pi 1161 1128 phi_max_minor = flip_phi(self.phi_max) - Pi 1162 1129 1163 1130 # check for minor sector 1164 1131 if phi_min_minor > phi_max_minor: … … 1169 1136 (phi_max_minor >= phi_data) 1170 1137 out = out_major + out_minor 1171 1138 1172 1139 return out
Note: See TracChangeset
for help on using the changeset viewer.