Changeset fb198a9 in sasview
- Timestamp:
- Jan 12, 2009 2:09:27 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:
- e23a20c
- Parents:
- db709e4
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
DataLoader/manipulations.py
rf2eee4a rfb198a9 634 634 635 635 636 class _Sector :636 class _Sectorold: 637 637 """ 638 638 Defines a sector region on a 2D data set. 639 639 The sector is defined by r_min, r_max, phi_min, phi_max, 640 and the position of the center of the ring. 641 640 and the position of the center of the ring. 642 641 Phi is defined between 0 and 2pi 643 642 """ … … 659 658 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 660 659 raise RuntimeError, "Ring averaging only take plottable_2D objects" 661 662 data = data2D.data 660 661 data = data2D.data 662 qmax = self.r_max 663 663 qmin = self.r_min 664 qmax = self.r_max665 664 666 665 if len(data2D.detector) != 1: … … 717 716 718 717 # Compute phi and check whether it's within the limits 719 phi_value = math.atan2(dy, dx)+math.pi 718 phi_value=math.atan2(dy,dx)+math.pi 719 # if phi_value<self.phi_min or phi_value>self.phi_max: 720 720 if phi_value<self.phi_min or phi_value>self.phi_max: 721 721 continue … … 758 758 return Data1D(x=x, y=y, dy=y_err) 759 759 760 760 class _Sector: 761 """ 762 Defines a sector region on a 2D data set. 763 The sector is defined by r_min, r_max, phi_min, phi_max, 764 and the position of the center of the ring 765 where phi_min and phi_max are defined by the right and left lines wrt central line 766 and phi_max could be less than phi_min. 767 768 Phi is defined between 0 and 2pi 769 """ 770 def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20): 771 self.r_min = r_min 772 self.r_max = r_max 773 self.phi_min = phi_min 774 self.phi_max = phi_max 775 self.nbins = nbins 776 777 def _agv(self, data2D, run='phi'): 778 """ 779 Perform sector averaging. 780 781 @param data2D: Data2D object 782 @param run: define the varying parameter ('phi' or 'q') 783 @return: Data1D object 784 """ 785 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 786 raise RuntimeError, "Ring averaging only take plottable_2D objects" 787 788 data = data2D.data 789 qmax = self.r_max 790 qmin = self.r_min 791 792 if len(data2D.detector) != 1: 793 raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector) 794 pixel_width_x = data2D.detector[0].pixel_size.x 795 pixel_width_y = data2D.detector[0].pixel_size.y 796 det_dist = data2D.detector[0].distance 797 wavelength = data2D.source.wavelength 798 center_x = data2D.detector[0].beam_center.x/pixel_width_x 799 center_y = data2D.detector[0].beam_center.y/pixel_width_y 800 801 y = numpy.zeros(self.nbins) 802 y_counts = numpy.zeros(self.nbins) 803 x = numpy.zeros(self.nbins) 804 y_err = numpy.zeros(self.nbins) 805 806 for i in range(numpy.size(data,1)): 807 dx = pixel_width_x*(i+0.5 - center_x) 808 809 # Min and max x-value for the pixel 810 minx = pixel_width_x*(i - center_x) 811 maxx = pixel_width_x*(i+1.0 - center_x) 812 813 for j in range(numpy.size(data,0)): 814 dy = pixel_width_y*(j+0.5 - center_y) 815 816 q_value = get_q(dx, dy, det_dist, wavelength) 817 818 # Min and max y-value for the pixel 819 miny = pixel_width_y*(j - center_y) 820 maxy = pixel_width_y*(j+1.0 - center_y) 821 822 # Calculate the q-value for each corner 823 # q_[x min or max][y min or max] 824 q_00 = get_q(minx, miny, det_dist, wavelength) 825 q_01 = get_q(minx, maxy, det_dist, wavelength) 826 q_10 = get_q(maxx, miny, det_dist, wavelength) 827 q_11 = get_q(maxx, maxy, det_dist, wavelength) 828 829 # Look for intercept between each side of the pixel 830 # and the constant-q ring for qmax 831 frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 832 833 # Look for intercept between each side of the pixel 834 # and the constant-q ring for qmin 835 frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 836 837 # We are interested in the region between qmin and qmax 838 # therefore the fraction of the surface of the pixel 839 # that we will use to calculate the number of counts to 840 # include is given by: 841 842 frac = frac_max - frac_min 843 844 # Compute phi and check whether it's within the limits 845 phi_value=math.atan2(dy,dx)+math.pi 846 if self.phi_max>2*math.pi: 847 self.phi_max=self.phi_max-2*math.pi 848 if self.phi_min<0: 849 self.phi_max=self.phi_max+2*math.pi 850 851 #In case of two ROI (symetric major and minor regions)(for 'q2') 852 if run.lower()=='q2': 853 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 854 temp_max=self.phi_max+math.pi 855 temp_min=self.phi_min+math.pi 856 else: 857 temp_max=self.phi_max 858 temp_min=self.phi_min 859 860 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 861 if (phi_value<temp_min or phi_value>temp_max): 862 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 863 continue 864 if (self.phi_max<self.phi_min): 865 tmp_max=self.phi_max+math.pi 866 tmp_min=self.phi_min-math.pi 867 else: 868 tmp_max=self.phi_max 869 tmp_min=self.phi_min 870 if (tmp_min<math.pi and tmp_max>math.pi): 871 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 872 continue 873 #In case of one ROI (major only)(for 'q' and 'phi') 874 else: 875 if (self.phi_max>=self.phi_min): 876 if (phi_value<self.phi_min or phi_value>self.phi_max): 877 continue 878 else: 879 if (phi_value<self.phi_min and phi_value>self.phi_max): 880 continue 881 882 # Check which type of averaging we need 883 if run.lower()=='phi': 884 i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 885 else: 886 # If we don't need this pixel, skip the rest of the work 887 #TODO: an improvement here would be to compute the average 888 # Q for the pixel from the part that is covered by 889 # the ring defined by q_min/q_max rather than the complete 890 # pixel 891 if q_value<self.r_min or q_value>self.r_max: 892 continue 893 i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1 894 895 try: 896 y[i_bin] += frac * data[j][i] 897 except: 898 import sys 899 print sys.exc_value 900 print i_bin, frac 901 902 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 903 y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 904 else: 905 y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 906 y_counts[i_bin] += frac 907 908 for i in range(self.nbins): 909 y[i] = y[i] / y_counts[i] 910 y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 911 # Check which type of averaging we need 912 if run.lower()=='phi': 913 x[i] = (self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min 914 else: 915 x[i] = (self.r_max-self.r_min)/self.nbins*(1.0*i + 0.5)+self.r_min 916 917 return Data1D(x=x, y=y, dy=y_err) 918 761 919 class SectorPhi(_Sector): 762 920 """ … … 776 934 return self._agv(data2D, 'phi') 777 935 778 class SectorQ (_Sector):936 class SectorQold(_Sector): 779 937 """ 780 938 Sector average as a function of Q. … … 792 950 """ 793 951 return self._agv(data2D, 'q') 794 952 953 class SectorQ(_Sector): 954 """ 955 Sector average as a function of Q for both symatric wings. 956 I(Q) is return and the data is averaged over phi. 957 958 A sector is defined by r_min, r_max, phi_min, phi_max. 959 r_min, r_max, phi_min, phi_max >0. 960 The number of bin in Q also has to be defined. 961 """ 962 def __call__(self, data2D): 963 """ 964 Perform sector average and return I(Q). 965 966 @param data2D: Data2D object 967 @return: Data1D object 968 """ 969 return self._agv(data2D, 'q2') 795 970 if __name__ == "__main__": 796 971
Note: See TracChangeset
for help on using the changeset viewer.