Changeset 3c67340 in sasview
- Timestamp:
- Jan 25, 2010 6:22:03 PM (15 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:
- 393f0f3
- Parents:
- 2f569b3
- File:
-
- 1 edited
-
DataLoader/manipulations.py (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
DataLoader/manipulations.py
r2f569b3 r3c67340 104 104 y = numpy.zeros(nbins) 105 105 err_y = numpy.zeros(nbins) 106 x_counts = numpy.zeros(nbins)107 106 y_counts = numpy.zeros(nbins) 108 107 … … 161 160 continue 162 161 163 x[i_q] += q_value162 x[i_q] = q_value 164 163 y[i_q] += frac * data2D.data[j][i] 165 164 if data2D.err_data == None or data2D.err_data[j][i]==0.0: … … 167 166 else: 168 167 err_y[i_q] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 169 x_counts[i_q] += 1170 168 y_counts[i_q] += frac 171 169 172 170 # Average the sums 173 171 for i in range(nbins): 174 if x_counts[i]>0 andy_counts[i]>0:172 if y_counts[i]>0: 175 173 err_y[i] = math.sqrt(err_y[i])/y_counts[i] 176 x[i] = x[i]/x_counts[i]177 174 y[i] = y[i]/y_counts[i] 178 175 … … 359 356 as a function of Q 360 357 """ 361 def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.00 05):358 def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.001): 362 359 # Minimum radius included in the average [A-1] 363 360 self.r_min = r_min … … 381 378 det_dist = data2D.detector[0].distance 382 379 wavelength = data2D.source.wavelength 383 center_x = data2D.detector[0].beam_center.x/pixel_width_x 384 center_y = data2D.detector[0].beam_center.y/pixel_width_y 380 center_x = data2D.detector[0].beam_center.x/pixel_width_x+0.5 381 center_y = data2D.detector[0].beam_center.y/pixel_width_y+0.5 385 382 386 383 # Find out the maximum Q range … … 396 393 397 394 qmax = get_q(dx_max, dy_max, det_dist, wavelength) 398 395 399 396 # Build array of Q intervals 400 397 nbins = int(math.ceil((qmax-self.r_min)/self.bin_width)) … … 402 399 403 400 x = numpy.zeros(nbins) 404 x_counts = numpy.zeros(nbins)405 401 y = numpy.zeros(nbins) 406 402 err_y = numpy.zeros(nbins) 407 403 y_counts = numpy.zeros(nbins) 408 404 409 405 for i in range(numpy.size(data2D.data,1)): 410 406 dx = pixel_width_x*(i+0.5 - center_x) 411 407 412 408 # Min and max x-value for the pixel … … 415 411 416 412 for j in range(numpy.size(data2D.data,0)): 417 413 dy = pixel_width_y*(j+0.5 - center_y) 414 415 q_value = get_q(dx, dy, det_dist, wavelength) 418 416 419 417 # Min and max y-value for the pixel … … 428 426 q_11 = get_q(maxx, maxy, det_dist, wavelength) 429 427 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 475 if q_value > qmax or q_value < self.r_min: 476 continue 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 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 489 454 # Average the sums 490 455 nzero = 0 491 492 456 for i in range(nbins): 493 #print x_counts[i],x[i],y[i] 494 if y_counts[i]>0 and x_counts[i]>0: 457 if y_counts[i]>0: 495 458 err_y[i] = math.sqrt(err_y[i])/y_counts[i] 496 459 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 499 460 else: 500 461 nzero += 1 … … 505 466 j=0 506 467 for i in range(nbins): 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 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 514 474 return Data1D(x=tx, y=ty, dy=terr_y) 515 475 … … 691 651 elif (q_00+q_01+q_10+q_11)/4.0 < qmax: 692 652 frac_max = 1.0 693 653 694 654 return frac_max 695 655 … … 723 683 return (q-q_1)/(q_0 - q_1) 724 684 return None 725 726 class _Sector: 685 686 #This class can be removed. 687 class _Sectorold: 727 688 """ 728 689 Defines a sector region on a 2D data set. 729 690 The sector is defined by r_min, r_max, phi_min, phi_max, 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 691 and the position of the center of the ring. 734 692 Phi is defined between 0 and 2pi 735 693 """ … … 768 726 y_counts = numpy.zeros(self.nbins) 769 727 x = numpy.zeros(self.nbins) 770 x_counts = numpy.zeros(self.nbins) 771 y_err = numpy.zeros(self.nbins) 728 y_err = numpy.zeros(self.nbins) 729 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) 772 856 773 857 # This If finds qmax within ROI defined by sector lines … … 775 859 temp0=1000000 776 860 temp1=0 777 for i in range(numpy.size(data,1)): 861 for i in range(numpy.size(data,1)): 862 dx = pixel_width_x*(i+0.5 - center_x) 778 863 for j in range(numpy.size(data,0)): 779 #number of sub-bin: default = 2 (ie., 1 bin: no sub-bin) 780 nsubbins = 1 781 782 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 783 for x_subbins in range(nsubbins): 784 for y_subbins in range(nsubbins): 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 785 906 786 dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 787 dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 788 789 q_value = get_q(dx, dy, det_dist, wavelength) 790 791 # Compute phi and check whether it's within the limits 792 phi_value=math.atan2(dy,dx)+math.pi 793 if self.phi_max>2*math.pi: 794 self.phi_max=self.phi_max-2*math.pi 795 if self.phi_min<0: 796 self.phi_max=self.phi_max+2*math.pi 797 798 #In case of two ROI (symmetric major and minor wings)(for 'q2') 799 if run.lower()=='q2': 800 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 801 temp_max=self.phi_max+math.pi 802 temp_min=self.phi_min+math.pi 803 else: 804 temp_max=self.phi_max 805 temp_min=self.phi_min 806 807 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 808 if (phi_value<temp_min or phi_value>temp_max): 809 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 810 continue 811 if (self.phi_max<self.phi_min): 812 tmp_max=self.phi_max+math.pi 813 tmp_min=self.phi_min-math.pi 814 else: 815 tmp_max=self.phi_max 816 tmp_min=self.phi_min 817 if (tmp_min<math.pi and tmp_max>math.pi): 818 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 819 continue 820 #In case of one ROI (major only)(i.e.,for 'q')and do nothing for 'phi'. 821 elif run.lower()=='q': 822 if (self.phi_max>=self.phi_min): 823 if (phi_value<self.phi_min or phi_value>self.phi_max): 824 continue 825 else: 826 if (phi_value<self.phi_min and phi_value>self.phi_max): 827 continue 828 if q_value<qmin or q_value>qmax: 829 continue 830 831 if run.lower()=='phi': 832 if temp1<phi_value: 833 temp1=phi_value 834 if temp0>phi_value: 835 temp0=phi_value 836 837 elif temp<q_value: 838 temp=q_value 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 839 915 840 916 if run.lower()=='phi': … … 845 921 #Beam center is already corrected, but the calculation below assumed it was not. 846 922 # Thus Beam center shifted back to uncorrected value. ToDo: cleanup the mess. 847 #center_x=center_x+0.5 848 #center_y=center_y+0.5 849 for i in range(numpy.size(data,1)): 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 850 932 for j in range(numpy.size(data,0)): 851 #number of sub-bin: default = 1 (ie., 1 bin: no sub-bin) 852 nsubbins = 1 853 854 #if is_intercept(self.r_max, q_00, q_01, q_10, q_11) \ 855 # or is_intercept(self.r_min, q_00, q_01, q_10, q_11): 856 # # Number of sub-bins near the boundary of ROI in x-direction 857 # nsubbins = 3 858 859 #Sub divide one pixel into nine sub-pixels only for the pixels where the qmax or qmin line crosses over. 860 for x_subbins in range(nsubbins): 861 for y_subbins in range(nsubbins): 862 863 dx = pixel_width_x*(i+(x_subbins)/nsubbins - center_x) 864 dy = pixel_width_y*(j+(y_subbins)/nsubbins - center_y) 865 866 q_value = get_q(dx, dy, det_dist, wavelength) 867 # Min and max x-value for the subpixel 868 minx = pixel_width_x*(i+(-0.5+x_subbins)/nsubbins - center_x) 869 maxx = pixel_width_x*(i+(0.5 +x_subbins)/nsubbins - center_x) 870 # Min and max y-value for the subpixel 871 miny = pixel_width_y*(j+(-0.5+y_subbins)/nsubbins - center_y) 872 maxy = pixel_width_y*(j+(0.5 +y_subbins)/nsubbins - center_y) 873 874 875 # Calculate the q-value for each corner 876 # q_[x min or max][y min or max] 877 q_00 = get_q(minx, miny, det_dist, wavelength) 878 q_01 = get_q(minx, maxy, det_dist, wavelength) 879 q_10 = get_q(maxx, miny, det_dist, wavelength) 880 q_11 = get_q(maxx, maxy, det_dist, wavelength) 881 882 # Compute phi and check whether it's within the limits 883 phi_value=math.atan2(dy,dx)+math.pi 884 if self.phi_max>2*math.pi: 885 self.phi_max=self.phi_max-2*math.pi 886 if self.phi_min<0: 887 self.phi_max=self.phi_max+2*math.pi 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)): 991 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 888 997 889 # Look for intercept between each side of the pixel 890 # and the constant-q ring for qmax 891 frac_max = get_pixel_fraction(qmax, q_00, q_01, q_10, q_11) 892 893 # Look for intercept between each side of the pixel 894 # and the constant-q ring for qmin 895 frac_min = get_pixel_fraction(qmin, q_00, q_01, q_10, q_11) 896 897 # We are interested in the region between qmin and qmax 898 # therefore the fraction of the surface of the pixel 899 # that we will use to calculate the number of counts to 900 # include is given by: 901 902 frac = frac_max - frac_min 903 904 #In case of two ROI (symmetric major and minor regions)(for 'q2') 905 if run.lower()=='q2': 906 if ((self.phi_max>=0 and self.phi_max<math.pi)and (self.phi_min>=0 and self.phi_min<math.pi)): 907 temp_max=self.phi_max+math.pi 908 temp_min=self.phi_min+math.pi 909 else: 910 temp_max=self.phi_max 911 temp_min=self.phi_min 912 913 if ((temp_max>=math.pi and temp_max<2*math.pi)and (temp_min>=math.pi and temp_min<2*math.pi)): 914 if (phi_value<temp_min or phi_value>temp_max): 915 if (phi_value<temp_min-math.pi or phi_value>temp_max-math.pi): 916 continue 917 if (self.phi_max<self.phi_min): 918 tmp_max=self.phi_max+math.pi 919 tmp_min=self.phi_min-math.pi 920 else: 921 tmp_max=self.phi_max 922 tmp_min=self.phi_min 923 if (tmp_min<math.pi and tmp_max>math.pi): 924 if((phi_value>tmp_max and phi_value<tmp_min+math.pi)or (phi_value>tmp_max-math.pi and phi_value<tmp_min)): 925 continue 926 #In case of one ROI (major only)(i.e.,for 'q' and 'phi') 927 else: 928 if (self.phi_max>=self.phi_min): 929 if (phi_value<self.phi_min or phi_value>self.phi_max): 930 continue 931 932 # Check which type of averaging we need 933 if run.lower()=='phi': 934 i_bin = int(math.floor(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) 935 else: 936 # If we don't need this pixel, skip the rest of the work 937 #TODO: an improvement here would be to compute the average 938 # Q for the pixel from the part that is covered by 939 # the ring defined by q_min/q_max rather than the complete 940 # pixel 941 i_bin = int(math.floor(self.nbins*(q_value-qmin)/(qmax-qmin))) 942 943 else: 944 if (phi_value<self.phi_min and phi_value>self.phi_max): 945 continue 946 #print "qmax=",qmax,qmin 947 948 if q_value<qmin or q_value>qmax: 998 else: 999 if (phi_value<self.phi_min and phi_value>self.phi_max): 949 1000 continue 950 951 # Check which type of averaging we need 952 if run.lower()=='phi': 953 i_bin = int(math.floor(self.nbins*(phi_value-self.phi_min)/(self.phi_max-self.phi_min))) 954 else: 955 # If we don't need this pixel, skip the rest of the work 956 #TODO: an improvement here would be to compute the average 957 # Q for the pixel from the part that is covered by 958 # the ring defined by q_min/q_max rather than the complete 959 # pixel 960 i_bin = int(math.floor(self.nbins*(q_value-qmin)/(qmax-qmin))) 961 962 #Make sure that no out of bound happens 963 #This is needed because of a numerical error in math.ceil(). 964 if i_bin >= self.nbins: 965 i_bin = self.nbins - 1 966 elif i_bin < 0: 967 i_bin = 0 968 969 try: 970 y[i_bin] += frac * data[j][i] 971 except: 972 import sys 973 print sys.exc_value 974 print i_bin, frac 975 976 if data2D.err_data == None or data2D.err_data[j][i]==0.0: 977 y_err[i_bin] += frac * frac * math.fabs(data2D.data[j][i]) 978 else: 979 y_err[i_bin] += frac * frac * data2D.err_data[j][i] * data2D.err_data[j][i] 980 y_counts[i_bin] += frac 981 982 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 983 1030 for i in range(self.nbins): 984 1031 y[i] = y[i] / y_counts[i] … … 1011 1058 """ 1012 1059 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 object 1074 @return: Data1D object 1075 """ 1076 return self._agv(data2D, 'q') 1013 1077 1014 1078 class SectorQ(_Sector):
Note: See TracChangeset
for help on using the changeset viewer.
