Changeset bb0b12c in sasview
- Timestamp:
- Jan 25, 2010 2:58:07 PM (14 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:
- 2f569b3
- Parents:
- 1702180
- File:
-
- 1 edited
-
DataLoader/manipulations.py (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
DataLoader/manipulations.py
r8ba103f rbb0b12c 104 104 y = numpy.zeros(nbins) 105 105 err_y = numpy.zeros(nbins) 106 x_counts = numpy.zeros(nbins) 106 107 y_counts = numpy.zeros(nbins) 107 108 … … 160 161 continue 161 162 162 x[i_q] = q_value163 x[i_q] += q_value 163 164 y[i_q] += frac * data2D.data[j][i] 164 165 if data2D.err_data == None or data2D.err_data[j][i]==0.0: … … 166 167 else: 167 168 err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 169 x_counts[i_q] += 1 168 170 y_counts[i_q] += frac 169 171 170 172 # Average the sums 171 173 for i in range(nbins): 172 if y_counts[i]>0:174 if x_counts[i]>0 and y_counts[i]>0: 173 175 err_y[i] = math.sqrt(err_y[i])/y_counts[i] 176 x[i] = x[i]/x_counts[i] 174 177 y[i] = y[i]/y_counts[i] 175 178 … … 356 359 as a function of Q 357 360 """ 358 def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.00 1):361 def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 359 362 # Minimum radius included in the average [A-1] 360 363 self.r_min = r_min … … 378 381 det_dist = data2D.detector[0].distance 379 382 wavelength = data2D.source.wavelength 380 center_x = data2D.detector[0].beam_center.x/pixel_width_x +0.5381 center_y = data2D.detector[0].beam_center.y/pixel_width_y +0.5383 center_x = data2D.detector[0].beam_center.x/pixel_width_x 384 center_y = data2D.detector[0].beam_center.y/pixel_width_y 382 385 383 386 # Find out the maximum Q range … … 393 396 394 397 qmax = get_q(dx_max, dy_max, det_dist, wavelength) 395 398 396 399 # Build array of Q intervals 397 400 nbins = int(math.ceil((qmax-self.r_min)/self.bin_width)) … … 399 402 400 403 x = numpy.zeros(nbins) 404 x_counts = numpy.zeros(nbins) 401 405 y = numpy.zeros(nbins) 402 406 err_y = numpy.zeros(nbins) 403 407 y_counts = numpy.zeros(nbins) 404 408 405 409 for i in range(numpy.size(data2D.data,1)): 406 dx = pixel_width_x*(i+0.5 - center_x)410 407 411 408 412 # Min and max x-value for the pixel … … 411 415 412 416 for j in range(numpy.size(data2D.data,0)): 413 dy = pixel_width_y*(j+0.5 - center_y) 414 415 q_value = get_q(dx, dy, det_dist, wavelength) 417 416 418 417 419 # Min and max y-value for the pixel … … 426 428 q_11 = get_q(maxx, maxy, det_dist, wavelength) 427 429 428 # Look for intercept between each side of the pixel 429 # and the constant-q ring for qmax 430 frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11) 431 432 # Look for intercept between each side of the pixel 433 # and the constant-q ring for qmin 434 frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11) 435 436 # We are interested in the region between qmin and qmax 437 # therefore the fraction of the surface of the pixel 438 # that we will use to calculate the number of counts to 439 # include is given by: 440 frac = frac_max - frac_min 441 442 i_q = int(math.ceil((q_value-self.r_min)/self.bin_width)) - 1 443 if q_value > qmax or q_value < self.r_min: 444 continue 445 446 x[i_q] = q_value 447 y[i_q] += frac * data2D.data[j][i] 448 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 449 err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i]) 450 else: 451 err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 452 y_counts[i_q] += frac 453 430 #number of sub-bin: default = 1 (ie., 1 bin: no sub-bin) 431 nsubbins = 2 432 433 #if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \ 434 # or is_intercept(self.r_min, q_00, q_01, q_10, q_11): 435 # # Number of sub-bins near the boundary of ROI in x-direction 436 # nsubbins = 3 437 438 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 439 for x_subbins in range(nsubbins): 440 for y_subbins in range(nsubbins): 441 442 dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 443 dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 444 445 q_value = get_q(dx, dy, det_dist, wavelength) 446 # Min and max x-value for the subpixel 447 minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x) 448 maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x) 449 # Min and max y-value for the subpixel 450 miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y) 451 maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y) 452 # Look for intercept between each side of the pixel 453 # and the constant-q ring for qmax 454 # Calculate the q-value for each corner of sub pixel 455 # q_[x min or max][y min or max] 456 q_00 = get_q(minx, miny, det_dist, wavelength) 457 q_01 = get_q(minx, maxy, det_dist, wavelength) 458 q_10 = get_q(maxx, miny, det_dist, wavelength) 459 q_11 = get_q(maxx, maxy, det_dist, wavelength) 460 frac_max = get_pixel_fraction(self.r_max, q_00, q_01, q_10, q_11)/(nsubbins*nsubbins) 461 462 # Look for intercept between each side of the pixel 463 # and the constant-q ring for qmin 464 frac_min = get_pixel_fraction(self.r_min, q_00, q_01, q_10, q_11)/(nsubbins*nsubbins) 465 466 # We are interested in the region between qmin and qmax 467 # therefore the fraction of the surface of the pixel 468 # that we will use to calculate the number of counts to 469 # include is given by: 470 frac = (frac_max - frac_min) #??? Todo: What if frac_max - fracmin <0? 471 472 473 i_q = int(math.floor((q_value-self.r_min)/self.bin_width)) #- 1 474 if q_value > qmax or q_value < self.r_min: 475 continue 476 #print q_value 477 478 x[i_q] += q_value 479 x_counts[i_q] += 1 480 481 y[i_q] += frac * data2D.data[j][i] 482 483 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 484 err_y[i_q] += frac * frac * math.fabs(data2D.data[j][i]) 485 else: 486 err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 487 y_counts[i_q] += frac 488 454 489 # Average the sums 455 490 nzero = 0 491 456 492 for i in range(nbins): 457 if y_counts[i]>0: 493 #print x_counts[i],x[i],y[i] 494 if y_counts[i]>0 and x_counts[i]>0: 458 495 err_y[i] = math.sqrt(err_y[i])/y_counts[i] 459 496 y[i] = y[i]/y_counts[i] 497 x[i] = (qmax-self.r_min)/nbins*(1.0*i + 0.5)+ self.r_min#x[i]/x_counts[i] 498 460 499 else: 461 500 nzero += 1 … … 466 505 j=0 467 506 for i in range(nbins): 468 if err_y[i] != 0 : 469 tx[j] = x[i] 470 ty[j] = y[i] 471 terr_y[j] = err_y[i] 472 j+=1 473 507 if err_y[i] != 0 and y_counts[i]>0 and x_counts[i]>0: 508 #Make sure wintin range 509 if x[i] <= self.r_max and x[i] >= self.r_min: 510 tx[j] = x[i] 511 ty[j] = y[i] 512 terr_y[j] = err_y[i] 513 j+=1 474 514 return Data1D(x=tx, y=ty, dy=terr_y) 475 515 … … 597 637 | | 598 638 y=0 +--------------+ 599 q_00 q_ 01639 q_00 q_10 600 640 601 641 x=0 x=1 … … 651 691 elif (q_00+q_01+q_10+q_11)/4.0 < qmax: 652 692 frac_max = 1.0 653 693 654 694 return frac_max 655 695 … … 683 723 return (q-q_1)/(q_0 - q_1) 684 724 return None 685 686 #This class can be removed. 687 class _Sectorold: 725 726 class _Sector: 688 727 """ 689 728 Defines a sector region on a 2D data set. 690 729 The sector is defined by r_min, r_max, phi_min, phi_max, 691 and the position of the center of the ring. 730 and the position of the center of the ring 731 where phi_min and phi_max are defined by the right and left lines wrt central line 732 and phi_max could be less than phi_min. 733 692 734 Phi is defined between 0 and 2pi 693 735 """ … … 726 768 y_counts = numpy.zeros(self.nbins) 727 769 x = numpy.zeros(self.nbins) 770 x_counts = numpy.zeros(self.nbins) 728 771 y_err = numpy.zeros(self.nbins) 729 772 730 for i in range(numpy.size(data,1)): 731 dx = pixel_width_x*(i+0.5 - center_x) 732 733 # Min and max x-value for the pixel 734 minx = pixel_width_x*(i - center_x) 735 maxx = pixel_width_x*(i+1.0 - center_x) 736 737 for j in range(numpy.size(data,0)): 738 dy = pixel_width_y*(j+0.5 - center_y) 739 740 q_value = get_q(dx, dy, det_dist, wavelength) 741 742 # Min and max y-value for the pixel 743 miny = pixel_width_y*(j - center_y) 744 maxy = pixel_width_y*(j+1.0 - center_y) 745 746 # Calculate the q-value for each corner 747 # q_[x min or max][y min or max] 748 q_00 = get_q(minx, miny, det_dist, wavelength) 749 q_01 = get_q(minx, maxy, det_dist, wavelength) 750 q_10 = get_q(maxx, miny, det_dist, wavelength) 751 q_11 = get_q(maxx, maxy, det_dist, wavelength) 752 753 # Look for intercept between each side of the pixel 754 # and the constant-q ring for qmax 755 frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 756 757 # Look for intercept between each side of the pixel 758 # and the constant-q ring for qmin 759 frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 760 761 # We are interested in the region between qmin and qmax 762 # therefore the fraction of the surface of the pixel 763 # that we will use to calculate the number of counts to 764 # include is given by: 765 766 frac = frac_max - frac_min 767 768 # Compute phi and check whether it's within the limits 769 phi_value=math.atan2(dy,dx)+math.pi 770 # if phi_value<self.phi_min or phi_value>self.phi_max: 771 if phi_value<self.phi_min or phi_value>self.phi_max: 772 continue 773 774 # Check which type of averaging we need 775 if run.lower()=='phi': 776 i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 777 else: 778 # If we don't need this pixel, skip the rest of the work 779 #TODO: an improvement here would be to compute the average 780 # Q for the pixel from the part that is covered by 781 # the ring defined by q_min/q_max rather than the complete 782 # pixel 783 if q_value<self.r_min or q_value>self.r_max: 784 continue 785 i_bin = int(math.ceil(self.nbins*(q_value-self.r_min)/(self.r_max-self.r_min))) - 1 786 787 try: 788 y[i_bin] += frac * data[j][i] 789 except: 790 import sys 791 print sys.exc_value 792 print i_bin, frac 793 794 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 795 y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 796 else: 797 y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 798 y_counts[i_bin] += frac 799 800 for i in range(self.nbins): 801 y[i] = y[i] / y_counts[i] 802 y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 803 # Check which type of averaging we need 804 if run.lower()=='phi': 805 x[i] = (self.phi_max-self.phi_min)/self.nbins*(1.0*i + 0.5)+self.phi_min 806 else: 807 x[i] = (self.r_max-self.r_min)/self.nbins*(1.0*i + 0.5)+self.r_min 808 809 return Data1D(x=x, y=y, dy=y_err) 810 811 class _Sector: 812 """ 813 Defines a sector region on a 2D data set. 814 The sector is defined by r_min, r_max, phi_min, phi_max, 815 and the position of the center of the ring 816 where phi_min and phi_max are defined by the right and left lines wrt central line 817 and phi_max could be less than phi_min. 818 819 Phi is defined between 0 and 2pi 820 """ 821 def __init__(self, r_min, r_max, phi_min, phi_max,nbins=20): 822 self.r_min = r_min 823 self.r_max = r_max 824 self.phi_min = phi_min 825 self.phi_max = phi_max 826 self.nbins = nbins 827 828 def _agv(self, data2D, run='phi'): 829 """ 830 Perform sector averaging. 831 832 @param data2D: Data2D object 833 @param run: define the varying parameter ('phi' or 'q') 834 @return: Data1D object 835 """ 836 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 837 raise RuntimeError, "Ring averaging only take plottable_2D objects" 838 839 data = data2D.data 840 qmax = self.r_max 841 qmin = self.r_min 842 843 if len(data2D.detector) != 1: 844 raise RuntimeError, "Ring averaging: invalid number of detectors: %g" % len(data2D.detector) 845 pixel_width_x = data2D.detector[0].pixel_size.x 846 pixel_width_y = data2D.detector[0].pixel_size.y 847 det_dist = data2D.detector[0].distance 848 wavelength = data2D.source.wavelength 849 center_x = data2D.detector[0].beam_center.x/pixel_width_x 850 center_y = data2D.detector[0].beam_center.y/pixel_width_y 851 852 y = numpy.zeros(self.nbins) 853 y_counts = numpy.zeros(self.nbins) 854 x = numpy.zeros(self.nbins) 855 y_err = numpy.zeros(self.nbins) 773 774 center_x=center_x 775 center_y=center_y 856 776 857 777 # This If finds qmax within ROI defined by sector lines … … 859 779 temp0=1000000 860 780 temp1=0 861 for i in range(numpy.size(data,1)): 862 dx = pixel_width_x*(i+0.5 - center_x) 781 for i in range(numpy.size(data,1)): 863 782 for j in range(numpy.size(data,0)): 864 865 dy = pixel_width_y*(j+0.5 - center_y) 866 q_value = get_q(dx, dy, det_dist, wavelength) 867 # Compute phi and check whether it's within the limits 868 phi_value=math.atan2(dy,dx)+math.pi 869 if self.phi_max>2*math.pi: 870 self.phi_max=self.phi_max-2*math.pi 871 if self.phi_min<0: 872 self.phi_max=self.phi_max+2*math.pi 873 874 #In case of two ROI (symmetric major and minor wings)(for 'q2') 875 if run.lower()=='q2': 876 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 877 temp_max=self.phi_max+math.pi 878 temp_min=self.phi_min+math.pi 879 else: 880 temp_max=self.phi_max 881 temp_min=self.phi_min 882 883 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 884 if (phi_value<temp_min or phi_value>temp_max): 885 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 886 continue 887 if (self.phi_max<self.phi_min): 888 tmp_max=self.phi_max+math.pi 889 tmp_min=self.phi_min-math.pi 890 else: 891 tmp_max=self.phi_max 892 tmp_min=self.phi_min 893 if (tmp_min<math.pi and tmp_max>math.pi): 894 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 895 continue 896 #In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'. 897 elif run.lower()=='q': 898 if (self.phi_max>=self.phi_min): 899 if (phi_value<self.phi_min or phi_value>self.phi_max): 900 continue 901 else: 902 if (phi_value<self.phi_min and phi_value>self.phi_max): 903 continue 904 if q_value<qmin or q_value>qmax: 905 continue 906 907 if run.lower()=='phi': 908 if temp1<phi_value: 909 temp1=phi_value 910 if temp0>phi_value: 911 temp0=phi_value 912 913 elif temp<q_value: 914 temp=q_value 783 #number of sub-bin: default = 2 (ie., 1 bin: no sub-bin) 784 nsubbins = 2 785 786 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 787 for x_subbins in range(nsubbins): 788 for y_subbins in range(nsubbins): 789 790 dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 791 dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 792 793 q_value = get_q(dx, dy, det_dist, wavelength) 794 795 # Compute phi and check whether it's within the limits 796 phi_value=math.atan2(dy,dx)+math.pi 797 if self.phi_max>2*math.pi: 798 self.phi_max=self.phi_max-2*math.pi 799 if self.phi_min<0: 800 self.phi_max=self.phi_max+2*math.pi 801 802 #In case of two ROI (symmetric major and minor wings)(for 'q2') 803 if run.lower()=='q2': 804 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 805 temp_max=self.phi_max+math.pi 806 temp_min=self.phi_min+math.pi 807 else: 808 temp_max=self.phi_max 809 temp_min=self.phi_min 810 811 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 812 if (phi_value<temp_min or phi_value>temp_max): 813 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 814 continue 815 if (self.phi_max<self.phi_min): 816 tmp_max=self.phi_max+math.pi 817 tmp_min=self.phi_min-math.pi 818 else: 819 tmp_max=self.phi_max 820 tmp_min=self.phi_min 821 if (tmp_min<math.pi and tmp_max>math.pi): 822 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 823 continue 824 #In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'. 825 elif run.lower()=='q': 826 if (self.phi_max>=self.phi_min): 827 if (phi_value<self.phi_min or phi_value>self.phi_max): 828 continue 829 else: 830 if (phi_value<self.phi_min and phi_value>self.phi_max): 831 continue 832 if q_value<qmin or q_value>qmax: 833 continue 834 835 if run.lower()=='phi': 836 if temp1<phi_value: 837 temp1=phi_value 838 if temp0>phi_value: 839 temp0=phi_value 840 841 elif temp<q_value: 842 temp=q_value 915 843 916 844 if run.lower()=='phi': … … 921 849 #Beam center is already corrected, but the calculation below assumed it was not. 922 850 # Thus Beam center shifted back to uncorrected value. ToDo: cleanup the mess. 923 center_x=center_x+0.5 924 center_y=center_y+0.5 925 for i in range(numpy.size(data,1)): 926 dx = pixel_width_x*(i+0.5 - center_x) 927 928 # Min and max x-value for the pixel 929 minx = pixel_width_x*(i - center_x) 930 maxx = pixel_width_x*(i+1.0 - center_x) 931 851 #center_x=center_x+0.5 852 #center_y=center_y+0.5 853 for i in range(numpy.size(data,1)): 932 854 for j in range(numpy.size(data,0)): 933 dy = pixel_width_y*(j+0.5 - center_y) 934 935 q_value = get_q(dx, dy, det_dist, wavelength) 936 937 # Min and max y-value for the pixel 938 miny = pixel_width_y*(j - center_y) 939 maxy = pixel_width_y*(j+1.0 - center_y) 940 941 # Calculate the q-value for each corner 942 # q_[x min or max][y min or max] 943 q_00 = get_q(minx, miny, det_dist, wavelength) 944 q_01 = get_q(minx, maxy, det_dist, wavelength) 945 q_10 = get_q(maxx, miny, det_dist, wavelength) 946 q_11 = get_q(maxx, maxy, det_dist, wavelength) 947 948 # Compute phi and check whether it's within the limits 949 phi_value=math.atan2(dy,dx)+math.pi 950 if self.phi_max>2*math.pi: 951 self.phi_max=self.phi_max-2*math.pi 952 if self.phi_min<0: 953 self.phi_max=self.phi_max+2*math.pi 954 955 # Look for intercept between each side of the pixel 956 # and the constant-q ring for qmax 957 frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 958 959 # Look for intercept between each side of the pixel 960 # and the constant-q ring for qmin 961 frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 962 963 # We are interested in the region between qmin and qmax 964 # therefore the fraction of the surface of the pixel 965 # that we will use to calculate the number of counts to 966 # include is given by: 967 968 frac = frac_max - frac_min 969 970 #In case of two ROI (symmetric major and minor regions)(for 'q2') 971 if run.lower()=='q2': 972 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 973 temp_max=self.phi_max+math.pi 974 temp_min=self.phi_min+math.pi 975 else: 976 temp_max=self.phi_max 977 temp_min=self.phi_min 978 979 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 980 if (phi_value<temp_min or phi_value>temp_max): 981 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 982 continue 983 if (self.phi_max<self.phi_min): 984 tmp_max=self.phi_max+math.pi 985 tmp_min=self.phi_min-math.pi 986 else: 987 tmp_max=self.phi_max 988 tmp_min=self.phi_min 989 if (tmp_min<math.pi and tmp_max>math.pi): 990 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 855 #number of sub-bin: default = 1 (ie., 1 bin: no sub-bin) 856 nsubbins = 2 857 858 #if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \ 859 # or is_intercept(self.r_min, q_00, q_01, q_10, q_11): 860 # # Number of sub-bins near the boundary of ROI in x-direction 861 # nsubbins = 3 862 863 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 864 for x_subbins in range(nsubbins): 865 for y_subbins in range(nsubbins): 866 867 dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 868 dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 869 870 q_value = get_q(dx, dy, det_dist, wavelength) 871 # Min and max x-value for the subpixel 872 minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x) 873 maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x) 874 # Min and max y-value for the subpixel 875 miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y) 876 maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y) 877 878 879 # Calculate the q-value for each corner 880 # q_[x min or max][y min or max] 881 q_00 = get_q(minx, miny, det_dist, wavelength) 882 q_01 = get_q(minx, maxy, det_dist, wavelength) 883 q_10 = get_q(maxx, miny, det_dist, wavelength) 884 q_11 = get_q(maxx, maxy, det_dist, wavelength) 885 886 # Compute phi and check whether it's within the limits 887 phi_value=math.atan2(dy,dx)+math.pi 888 if self.phi_max>2*math.pi: 889 self.phi_max=self.phi_max-2*math.pi 890 if self.phi_min<0: 891 self.phi_max=self.phi_max+2*math.pi 892 893 # Look for intercept between each side of the pixel 894 # and the constant-q ring for qmax 895 frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 896 897 # Look for intercept between each side of the pixel 898 # and the constant-q ring for qmin 899 frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 900 901 # We are interested in the region between qmin and qmax 902 # therefore the fraction of the surface of the pixel 903 # that we will use to calculate the number of counts to 904 # include is given by: 905 906 frac = frac_max - frac_min 907 908 #In case of two ROI (symmetric major and minor regions)(for 'q2') 909 if run.lower()=='q2': 910 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 911 temp_max=self.phi_max+math.pi 912 temp_min=self.phi_min+math.pi 913 else: 914 temp_max=self.phi_max 915 temp_min=self.phi_min 916 917 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 918 if (phi_value<temp_min or phi_value>temp_max): 919 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 920 continue 921 if (self.phi_max<self.phi_min): 922 tmp_max=self.phi_max+math.pi 923 tmp_min=self.phi_min-math.pi 924 else: 925 tmp_max=self.phi_max 926 tmp_min=self.phi_min 927 if (tmp_min<math.pi and tmp_max>math.pi): 928 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 929 continue 930 #In case of one ROI (major only)(i.e.,for 'q' and 'phi') 931 else: 932 if (self.phi_max>=self.phi_min): 933 if (phi_value<self.phi_min or phi_value>self.phi_max): 934 continue 935 936 # Check which type of averaging we need 937 if run.lower()=='phi': 938 i_bin = int(math.floor(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) 939 else: 940 # If we don't need this pixel, skip the rest of the work 941 #TODO: an improvement here would be to compute the average 942 # Q for the pixel from the part that is covered by 943 # the ring defined by q_min/q_max rather than the complete 944 # pixel 945 i_bin = int(math.floor(self.nbins*(q_value-qmin)/(qmax-qmin))) 946 947 else: 948 if (phi_value<self.phi_min and phi_value>self.phi_max): 949 continue 950 #print "qmax=",qmax,qmin 951 952 if q_value<qmin or q_value>qmax: 991 953 continue 992 #In case of one ROI (major only)(i.e.,for 'q' and 'phi') 993 else: 994 if (self.phi_max>=self.phi_min): 995 if (phi_value<self.phi_min or phi_value>self.phi_max): 996 continue 997 998 else: 999 if (phi_value<self.phi_min and phi_value>self.phi_max): 1000 continue 1001 #print "qmax=",qmax,qmin 1002 1003 if q_value<qmin or q_value>qmax: 1004 continue 1005 1006 # Check which type of averaging we need 1007 if run.lower()=='phi': 1008 i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 1009 else: 1010 # If we don't need this pixel, skip the rest of the work 1011 #TODO: an improvement here would be to compute the average 1012 # Q for the pixel from the part that is covered by 1013 # the ring defined by q_min/q_max rather than the complete 1014 # pixel 1015 i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1 1016 1017 try: 1018 y[i_bin] += frac * data[j][i] 1019 except: 1020 import sys 1021 print sys.exc_value 1022 print i_bin, frac 1023 1024 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 1025 y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 1026 else: 1027 y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 1028 y_counts[i_bin] += frac 1029 954 955 # Check which type of averaging we need 956 if run.lower()=='phi': 957 i_bin = int(math.ceil(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) - 1 958 else: 959 # If we don't need this pixel, skip the rest of the work 960 #TODO: an improvement here would be to compute the average 961 # Q for the pixel from the part that is covered by 962 # the ring defined by q_min/q_max rather than the complete 963 # pixel 964 i_bin = int(math.ceil(self.nbins*(q_value-qmin)/(qmax-qmin))) - 1 965 966 try: 967 y[i_bin] += frac * data[j][i] 968 except: 969 import sys 970 print sys.exc_value 971 print i_bin, frac 972 973 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 974 y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 975 else: 976 y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 977 y_counts[i_bin] += frac 978 979 1030 980 for i in range(self.nbins): 1031 981 y[i] = y[i] / y_counts[i] … … 1058 1008 """ 1059 1009 return self._agv(data2D, 'phi') 1060 1061 class SectorQold(_Sector):1062 """1063 Sector average as a function of Q.1064 I(Q) is return and the data is averaged over phi.1065 1066 A sector is defined by r_min, r_max, phi_min, phi_max.1067 The number of bin in Q also has to be defined.1068 """1069 def __call__(self, data2D):1070 """1071 Perform sector average and return I(Q).1072 1073 @param data2D: Data2D object1074 @return: Data1D object1075 """1076 return self._agv(data2D, 'q')1077 1010 1078 1011 class SectorQ(_Sector):
Note: See TracChangeset
for help on using the changeset viewer.
