Changeset 2e83ff3 in sasview for DataLoader/manipulations.py
- Timestamp:
- Sep 22, 2008 3:53:09 PM (16 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:
- 1374174f
- Parents:
- 55fd102
- File:
-
- 1 edited
-
DataLoader/manipulations.py (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
DataLoader/manipulations.py
r70975f3 r2e83ff3 15 15 """ 16 16 #TODO: copy the meta data from the 2D object to the resulting 1D object 17 #TODO: Don't assume square pixels18 17 19 18 from data_info import plottable_2D, Data1D … … 222 221 raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector) 223 222 224 pixel_width = data2D.detector[0].pixel_size.x 223 pixel_width_x = data2D.detector[0].pixel_size.x 224 pixel_width_y = data2D.detector[0].pixel_size.y 225 225 det_dist = data2D.detector[0].distance 226 226 wavelength = data2D.source.wavelength 227 center_x = data2D.detector[0].beam_center.x/pixel_width 228 center_y = data2D.detector[0].beam_center.y/pixel_width 227 center_x = data2D.detector[0].beam_center.x/pixel_width_x 228 center_y = data2D.detector[0].beam_center.y/pixel_width_y 229 229 230 230 y = 0.0 … … 232 232 y_counts = 0.0 233 233 234 for i in range( len(data2D.data)):234 for i in range(numpy.size(data2D.data,1)): 235 235 # Min and max x-value for the pixel 236 minx = pixel_width *(i - center_x)237 maxx = pixel_width *(i+1.0 - center_x)236 minx = pixel_width_x*(i - center_x) 237 maxx = pixel_width_x*(i+1.0 - center_x) 238 238 239 239 qxmin = get_q(minx, 0.0, det_dist, wavelength) … … 245 245 frac_x = frac_max - frac_min 246 246 247 for j in range( len(data2D.data)):247 for j in range(numpy.size(data2D.data,0)): 248 248 # Min and max y-value for the pixel 249 miny = pixel_width *(j - center_y)250 maxy = pixel_width *(j+1.0 - center_y)249 miny = pixel_width_y*(j - center_y) 250 maxy = pixel_width_y*(j+1.0 - center_y) 251 251 252 252 qymin = get_q(0.0, miny, det_dist, wavelength) … … 339 339 raise RuntimeError, "Circular averaging: invalid number of detectors: %g" % len(data2D.detector) 340 340 341 pixel_width = data2D.detector[0].pixel_size.x 341 pixel_width_x = data2D.detector[0].pixel_size.x 342 pixel_width_y = data2D.detector[0].pixel_size.y 342 343 det_dist = data2D.detector[0].distance 343 344 wavelength = data2D.source.wavelength 344 center_x = data2D.detector[0].beam_center.x/pixel_width 345 center_y = data2D.detector[0].beam_center.y/pixel_width 345 center_x = data2D.detector[0].beam_center.x/pixel_width_x 346 center_y = data2D.detector[0].beam_center.y/pixel_width_y 346 347 347 348 # Find out the maximum Q range 348 xwidth = numpy.size(data2D.data,1)*pixel_width 349 xwidth = numpy.size(data2D.data,1)*pixel_width_x 349 350 dx_max = xwidth - data2D.detector[0].beam_center.x 350 351 if xwidth-dx_max>dx_max: 351 352 dx_max = xwidth-dx_max 352 353 353 ywidth = numpy.size(data2D.data,0)*pixel_width 354 ywidth = numpy.size(data2D.data,0)*pixel_width_y 354 355 dy_max = ywidth - data2D.detector[0].beam_center.y 355 356 if ywidth-dy_max>dy_max: … … 367 368 y_counts = numpy.zeros(nbins) 368 369 369 for i in range( len(data2D.data)):370 dx = pixel_width *(i+0.5 - center_x)370 for i in range(numpy.size(data2D.data,1)): 371 dx = pixel_width_x*(i+0.5 - center_x) 371 372 372 373 # Min and max x-value for the pixel 373 minx = pixel_width *(i - center_x)374 maxx = pixel_width *(i+1.0 - center_x)375 376 for j in range( len(data2D.data)):377 dy = pixel_width *(j+0.5 - center_y)374 minx = pixel_width_x*(i - center_x) 375 maxx = pixel_width_x*(i+1.0 - center_x) 376 377 for j in range(numpy.size(data2D.data,0)): 378 dy = pixel_width_y*(j+0.5 - center_y) 378 379 379 380 q_value = get_q(dx, dy, det_dist, wavelength) 380 381 381 382 # Min and max y-value for the pixel 382 miny = pixel_width *(j - center_y)383 maxy = pixel_width *(j+1.0 - center_y)383 miny = pixel_width_y*(j - center_y) 384 maxy = pixel_width_y*(j+1.0 - center_y) 384 385 385 386 # Calculate the q-value for each corner … … 463 464 if len(data2D.detector) != 1: 464 465 raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector) 465 pixel_width = data2D.detector[0].pixel_size.x 466 pixel_width_x = data2D.detector[0].pixel_size.x 467 pixel_width_y = data2D.detector[0].pixel_size.y 466 468 det_dist = data2D.detector[0].distance 467 469 wavelength = data2D.source.wavelength 468 center_x = self.center_x/pixel_width 469 center_y = self.center_y/pixel_width 470 #center_x = self.center_x/pixel_width_x 471 #center_y = self.center_y/pixel_width_y 472 center_x = data2D.detector[0].beam_center.x/pixel_width_x 473 center_y = data2D.detector[0].beam_center.y/pixel_width_y 474 470 475 471 476 phi_bins = numpy.zeros(self.nbins_phi) … … 474 479 phi_err = numpy.zeros(self.nbins_phi) 475 480 476 for i in range( len(data)):477 dx = pixel_width *(i+0.5 - center_x)481 for i in range(numpy.size(data,1)): 482 dx = pixel_width_x*(i+0.5 - center_x) 478 483 479 484 # Min and max x-value for the pixel 480 minx = pixel_width *(i - center_x)481 maxx = pixel_width *(i+1.0 - center_x)482 483 for j in range( len(data)):484 dy = pixel_width *(j+0.5 - center_y)485 minx = pixel_width_x*(i - center_x) 486 maxx = pixel_width_x*(i+1.0 - center_x) 487 488 for j in range(numpy.size(data,0)): 489 dy = pixel_width_y*(j+0.5 - center_y) 485 490 486 491 q_value = get_q(dx, dy, det_dist, wavelength) 487 492 488 493 # Min and max y-value for the pixel 489 miny = pixel_width *(j - center_y)490 maxy = pixel_width *(j+1.0 - center_y)494 miny = pixel_width_y*(j - center_y) 495 maxy = pixel_width_y*(j+1.0 - center_y) 491 496 492 497 # Calculate the q-value for each corner … … 512 517 frac = frac_max - frac_min 513 518 514 i_phi = int(math.ceil(self.nbins_phi*(math.atan2(dy, dx)+math.pi)/(2.0*math.pi)) - 1)519 i_phi = int(math.ceil(self.nbins_phi*(math.atan2(dy, dx)+math.pi)/(2.0*math.pi))) - 1 515 520 516 521 phi_bins[i_phi] += frac * data[j][i] … … 629 634 630 635 631 class Sector:636 class _Sector: 632 637 """ 633 638 Defines a sector region on a 2D data set. 634 639 The sector is defined by r_min, r_max, phi_min, phi_max, 635 640 and the position of the center of the ring. 636 """ 637 pass 641 642 Phi is defined between 0 and 2pi 643 """ 644 def __init__(self, r_min, r_max, phi_min, phi_max): 645 self.r_min = r_min 646 self.r_max = r_max 647 self.phi_min = phi_min 648 self.phi_max = phi_max 649 self.nbins = 20 650 651 def _agv(self, data2D, run='phi'): 652 """ 653 Perform sector averaging. 654 655 @param data2D: Data2D object 656 @param run: define the varying parameter ('phi' or 'q') 657 @return: Data1D object 658 """ 659 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 660 raise RuntimeError, "Ring averaging only take plottable_2D objects" 661 662 data = data2D.data 663 qmin = self.r_min 664 qmax = self.r_max 665 666 if len(data2D.detector) != 1: 667 raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector) 668 pixel_width_x = data2D.detector[0].pixel_size.x 669 pixel_width_y = data2D.detector[0].pixel_size.y 670 det_dist = data2D.detector[0].distance 671 wavelength = data2D.source.wavelength 672 center_x = data2D.detector[0].beam_center.x/pixel_width_x 673 center_y = data2D.detector[0].beam_center.y/pixel_width_y 674 675 y = numpy.zeros(self.nbins) 676 y_counts = numpy.zeros(self.nbins) 677 x = numpy.zeros(self.nbins) 678 y_err = numpy.zeros(self.nbins) 679 680 for i in range(numpy.size(data,1)): 681 dx = pixel_width_x*(i+0.5 - center_x) 682 683 # Min and max x-value for the pixel 684 minx = pixel_width_x*(i - center_x) 685 maxx = pixel_width_x*(i+1.0 - center_x) 686 687 for j in range(numpy.size(data,0)): 688 dy = pixel_width_y*(j+0.5 - center_y) 689 690 q_value = get_q(dx, dy, det_dist, wavelength) 691 692 # Min and max y-value for the pixel 693 miny = pixel_width_y*(j - center_y) 694 maxy = pixel_width_y*(j+1.0 - center_y) 695 696 # Calculate the q-value for each corner 697 # q_[x min or max][y min or max] 698 q_00 = get_q(minx, miny, det_dist, wavelength) 699 q_01 = get_q(minx, maxy, det_dist, wavelength) 700 q_10 = get_q(maxx, miny, det_dist, wavelength) 701 q_11 = get_q(maxx, maxy, det_dist, wavelength) 702 703 # Look for intercept between each side of the pixel 704 # and the constant-q ring for qmax 705 frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 706 707 # Look for intercept between each side of the pixel 708 # and the constant-q ring for qmin 709 frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 710 711 # We are interested in the region between qmin and qmax 712 # therefore the fraction of the surface of the pixel 713 # that we will use to calculate the number of counts to 714 # include is given by: 715 716 frac = frac_max - frac_min 717 718 # Compute phi and check whether it's within the limits 719 phi_value = math.atan2(dy, dx)+math.pi 720 if phi_value<self.phi_min or phi_value>self.phi_max: 721 continue 722 723 # Check which type of averaging we need 724 if run.lower()=='phi': 725 i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 726 else: 727 # If we don't need this pixel, skip the rest of the work 728 #TODO: an improvement here would be to compute the average 729 # Q for the pixel from the part that is covered by 730 # the ring defined by q_min/q_max rather than the complete 731 # pixel 732 if q_value<self.r_min or q_value>self.r_max: 733 continue 734 i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1 735 736 try: 737 y[i_bin] += frac * data[j][i] 738 except: 739 import sys 740 print sys.exc_value 741 print i_bin, frac 742 743 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 744 y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 745 else: 746 y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 747 y_counts[i_bin] += frac 748 749 for i in range(self.nbins): 750 y[i] = y[i] / y_counts[i] 751 y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 752 # Check which type of averaging we need 753 if run.lower()=='phi': 754 x[i] = (self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min 755 else: 756 x[i] = (self.r_max-self.r_min)/self.nbins*(1.0*i + 0.5)+self.r_min 757 758 return Data1D(x=x, y=y, dy=y_err) 759 760 761 class SectorPhi(_Sector): 762 """ 763 Sector average as a function of phi. 764 I(phi) is return and the data is averaged over Q. 765 766 A sector is defined by r_min, r_max, phi_min, phi_max. 767 The number of bin in phi also has to be defined. 768 """ 769 def __call__(self, data2D): 770 """ 771 Perform sector average and return I(phi). 772 773 @param data2D: Data2D object 774 @return: Data1D object 775 """ 776 return self._agv(data2D, 'phi') 777 778 class SectorQ(_Sector): 779 """ 780 Sector average as a function of Q. 781 I(Q) is return and the data is averaged over phi. 782 783 A sector is defined by r_min, r_max, phi_min, phi_max. 784 The number of bin in Q also has to be defined. 785 """ 786 def __call__(self, data2D): 787 """ 788 Perform sector average and return I(Q). 789 790 @param data2D: Data2D object 791 @return: Data1D object 792 """ 793 return self._agv(data2D, 'q') 638 794 639 795 if __name__ == "__main__": … … 646 802 647 803 648 #r = Boxsum(x_min=.2, x_max=.4, y_min=0.2, y_max=0.4) 649 r = SlabX(x_min=-.01, x_max=.01, y_min=-0.0002, y_max=0.0002, bin_width=0.0004) 650 r.fold = False 804 r = SectorQ(r_min=.005, r_max=.01, phi_min=0.0, phi_max=math.pi/2.0) 651 805 o = r(d) 652 806 653 #s = SlabY(x_min=-.01, x_max=.01, y_min=-0.0002, y_max=0.0002, bin_width=0.0004) 654 #s.fold = False 655 #p = s(d) 807 s = Ring(r_min=.005, r_max=.01) 808 p = s(d) 656 809 657 810 for i in range(len(o.x)):
Note: See TracChangeset
for help on using the changeset viewer.
