Changeset 8a621ac in sasview for sanscalculator/src/sans/calculator/resolution_calculator.py
- Timestamp:
- Apr 27, 2012 9:22:40 AM (12 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:
- 10bfeb3
- Parents:
- 34f3ad0
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sanscalculator/src/sans/calculator/resolution_calculator.py
r5f3164c r8a621ac 1 1 """ 2 2 This object is a small tool to allow user to quickly 3 determine the variance in q from the 3 determine the variance in q from the 4 4 instrumental parameters. 5 5 """ … … 11 11 from math import pi 12 12 from math import sqrt 13 import math 14 import scipy 13 import math 15 14 import numpy 16 15 … … 19 18 #Gravitational acc. in cgs unit 20 19 _GRAVITY = 981.0 20 21 21 22 22 class ResolutionCalculator(object): … … 86 86 self.qyrange = [] 87 87 88 def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 89 qy_min, qy_max, coord ='cartesian'):88 def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, 89 qy_min, qy_max, coord='cartesian'): 90 90 """ 91 91 Compute the resolution … … 98 98 # wavelength etc. 99 99 lamda_list, dlamb_list = self.get_wave_list() 100 intens_list = []#self.get_intensity_list() 101 100 intens_list = [] 102 101 sig1_list = [] 103 102 sig2_list = [] … … 120 119 self.compute(lam, dlam, qx_value, qy_value, coord, tof) 121 120 # make image 122 image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 121 image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, 123 122 sigma_r, qx_min, qx_max, qy_min, qy_max, 124 123 coord, False) 125 124 126 # Non tof mode to be speed up 125 # Non tof mode to be speed up 127 126 #if num_lamda < 2: 128 127 # return self.plot_image(image) … … 140 139 self.qxrange = [qx_min, qx_max] 141 140 self.qyrange = [qy_min, qy_max] 142 #print qy_max+qy_min,qy_max,qy_min143 141 sig1_list.append(sigma_1) 144 142 sig2_list.append(sigma_2) 145 143 sigr_list.append(sigma_r) 146 144 sigma1d_list.append(sigma1d) 147 # redraw image in global 2d q-space. 145 # redraw image in global 2d q-space. 148 146 self.image_lam = [] 149 147 total_intensity = 0 … … 156 154 dlam = dlamb_list[ind] 157 155 intens = self.setup_tof(lam, dlam) 158 out = self.get_image(qx_value, qy_value, sig1_list[ind], 159 sig2_list[ind], sigr_list[ind], 156 out = self.get_image(qx_value, qy_value, sig1_list[ind], 157 sig2_list[ind], sigr_list[ind], 160 158 qx_min, qx_max, qy_min, qy_max, coord) 161 159 # this is the case of q being outside the detector … … 164 162 image = out 165 163 # set variance as sigmas 166 sigma_1 += sig1_list[ind] * 164 sigma_1 += sig1_list[ind] * sig1_list[ind] * self.intensity 167 165 sigma_r += sigr_list[ind] * sigr_list[ind] * self.intensity 168 166 sigma_2 += sig2_list[ind] * sig2_list[ind] * self.intensity … … 182 180 self.sigma_2 = sqrt(sigma_2) 183 181 self.sigma_1d = sqrt(sigma1d) 184 # rescale 185 max_im_val = 1 #image_out.max()182 # rescale 183 max_im_val = 1 186 184 if max_im_val > 0: 187 185 image_out /= max_im_val … … 199 197 200 198 # plot image 201 return self.plot_image(self.image) 199 return self.plot_image(self.image) 202 200 203 201 def setup_tof(self, wavelength, wavelength_spread): … … 214 212 215 213 if wavelength == 0: 216 msg = "Can't compute the resolution: the wavelength is zero..." 214 msg = "Can't compute the resolution: the wavelength is zero..." 217 215 raise RuntimeError, msg 218 216 return self.intensity 219 217 220 def compute(self, wavelength, wavelength_spread, qx_value, qy_value, 221 coord ='cartesian', tof=False):218 def compute(self, wavelength, wavelength_spread, qx_value, qy_value, 219 coord='cartesian', tof=False): 222 220 """ 223 221 Compute the Q resoltuion in || and + direction of 2D … … 232 230 if tof: 233 231 # rectangular 234 tof_factor = 232 tof_factor = 2 235 233 else: 236 234 # triangular 237 tof_factor = 235 tof_factor = 1 238 236 # Find polar values 239 237 qr_value, phi = self._get_polar_value(qx_value, qy_value) … … 285 283 sigma_1 += self.get_variance(rthree, l_two, phi, comp1) 286 284 # for gravity term for 1d 287 sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread,288 phi, comp1, 'on') / tof_factor285 sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, 286 lamb_spread, phi, comp1, 'on') / tof_factor 289 287 # for wavelength spread 290 288 # reserve for 1d calculation 291 289 A_value = self._cal_A_value(lamb, l_ssa, l_sad) 292 sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value, 293 radius, l_two, lamb_spread, 290 sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value, 291 radius, l_two, lamb_spread, 294 292 phi, 'radial', 'on') 295 293 sigma_wave_1 /= tof_factor 296 sigma_wave_1_1d /= 294 sigma_wave_1_1d /= tof_factor 297 295 # for 1d 298 variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d296 variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d 299 297 # normalize 300 298 variance_1d_1 = knot * knot * variance_1d_1 / 12 … … 303 301 #sigma_1 += sigma_wave_1 304 302 # normalize 305 sigma_1 = knot *sqrt(sigma_1 / 12)306 sigma_r = knot *sqrt(sigma_wave_1 / (tof_factor *12))303 sigma_1 = knot * sqrt(sigma_1 / 12) 304 sigma_r = knot * sqrt(sigma_wave_1 / (tof_factor *12)) 307 305 # sigma in the phi/y direction 308 306 # for source apperture … … 316 314 317 315 # for gravity term for 1d 318 sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, 319 phi, comp2, 'on') / tof_factor 320 321 316 sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, 317 lamb_spread, phi, comp2, 'on') / tof_factor 318 322 319 # for wavelength spread 323 320 # reserve for 1d calculation 324 sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value, 325 radius, l_two, lamb_spread, 326 phi, 'phi', 'on') 327 sigma_wave_2 /= 328 sigma_wave_2_1d /= 321 sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value, 322 radius, l_two, lamb_spread, 323 phi, 'phi', 'on') 324 sigma_wave_2 /= tof_factor 325 sigma_wave_2_1d /= tof_factor 329 326 # for 1d 330 327 variance_1d_2 = (sigma_2 + sigma_2grav1d) / 2 + sigma_wave_2_1d 331 328 # normalize 332 variance_1d_2 = knot *knot*variance_1d_2 / 12329 variance_1d_2 = knot * knot * variance_1d_2 / 12 333 330 334 331 # for 2d … … 336 333 #sigma_2 += sigma_wave_2 337 334 # normalize 338 sigma_2 = 335 sigma_2 = knot * sqrt(sigma_2 / 12) 339 336 sigma1d = sqrt(variance_1d_1 + variance_1d_2) 340 337 # set sigmas … … 345 342 return qr_value, phi, sigma_1, sigma_2, sigma_r, sigma1d 346 343 347 def _within_detector_range(self, qx_value, qy_value):344 def _within_detector_range(self, qx_value, qy_value): 348 345 """ 349 346 check if qvalues are within detector range … … 369 366 370 367 def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r, 371 qx_min, qx_max, qy_min, qy_max, 372 coord = 'cartesian', full_cal=True):368 qx_min, qx_max, qy_min, qy_max, 369 coord='cartesian', full_cal=True): 373 370 """ 374 371 Get the resolution in polar coordinate ready to plot … … 380 377 """ 381 378 # Get qx_max and qy_max... 382 output =self._get_detector_qxqy_pixels()379 self._get_detector_qxqy_pixels() 383 380 384 381 qr_value, phi = self._get_polar_value(qx_value, qy_value) 385 382 386 383 # Check whether the q value is within the detector range 387 msg = "Invalid input: Q value out of the detector range..."384 #msg = "Invalid input: Q value out of the detector range..." 388 385 if qx_min < self.qx_min: 389 386 self.qx_min = qx_min … … 416 413 qc_2 = 0.0 417 414 # Calculate the 2D Gaussian distribution image 418 image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, 415 image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, 419 416 sigma_1, sigma_2, sigma_r) 420 417 else: … … 426 423 427 424 # Calculate the 2D Gaussian distribution image 428 image = self._gaussian2d(q_1, q_2, qc_1, qc_2, 425 image = self._gaussian2d(q_1, q_2, qc_1, qc_2, 429 426 sigma_1, sigma_2, sigma_r) 430 427 # out side of detector … … 445 442 """ 446 443 Plot image using pyplot 447 : image: 2d resolution image 444 : image: 2d resolution image 448 445 449 446 : return plt: pylab object … … 455 452 plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$') 456 453 # Max value of the image 457 max = numpy.max(image)454 # max = numpy.max(image) 458 455 qx_min, qx_max, qy_min, qy_max = self.get_detector_qrange() 459 456 460 457 # Image 461 im = plt.imshow(image, 462 extent =[qx_min, qx_max, qy_min, qy_max])458 im = plt.imshow(image, 459 extent=[qx_min, qx_max, qy_min, qy_max]) 463 460 464 461 # bilinear interpolation to make it smoother … … 473 470 self.image = [] 474 471 475 def get_variance(self, size = [], distance = 0, phi = 0, comp ='radial'):472 def get_variance(self, size=[], distance=0, phi=0, comp='radial'): 476 473 """ 477 474 Get the variance when the slit/pinhole size is given 478 : size: list that can be one(diameter for circular) 475 : size: list that can be one(diameter for circular) 479 476 or two components(lengths for rectangular) 480 477 : distance: [z, x] where z along the incident beam, x // qx_value … … 482 479 483 480 : return variance: sigma^2 484 """ 481 """ 485 482 # check the length of size (list) 486 483 len_size = len(size) … … 509 506 # for rectangular slit 510 507 elif len_size == 2: 511 x_comp = size[0] * phi_x 508 x_comp = size[0] * phi_x 512 509 y_comp = size[1] * phi_y 513 510 # otherwise 514 511 else: 515 512 raise ValueError, " Improper input..." 516 # get them squared 517 sigma = x_comp * x_comp 513 # get them squared 514 sigma = x_comp * x_comp 518 515 sigma += y_comp * y_comp 519 516 # normalize by distance … … 522 519 return sigma 523 520 524 def get_variance_wave(self, A_value, radius, distance, spread, phi, 525 comp = 'radial', switch ='on'):521 def get_variance_wave(self, A_value, radius, distance, spread, phi, 522 comp='radial', switch='on'): 526 523 """ 527 524 Get the variance when the wavelength spread is given … … 547 544 sigma1d *= (math.sin(phi)*math.sin(phi)) 548 545 else: 549 sigma1d *= 1 550 # sigma^2 for 2d 546 sigma1d *= 1 547 # sigma^2 for 2d 551 548 # shift the coordinate due to the gravitational shift 552 549 rad_x = radius * math.cos(phi) 553 550 rad_y = A_value - radius * math.sin(phi) 554 551 radius = math.sqrt(rad_x * rad_x + rad_y * rad_y) 555 # new phi 552 # new phi 556 553 phi = math.atan2(-rad_y, rad_x) 557 self.gravity_phi = phi 554 self.gravity_phi = phi 558 555 # calculate sigma^2 559 556 sigma = 2 * math.pow(radius/distance*spread, 2) … … 563 560 sigma *= (math.sin(phi)*math.sin(phi)) 564 561 else: 565 sigma *= 1 562 sigma *= 1 566 563 567 564 return sigma, sigma1d 568 565 569 def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, 570 phi, comp = 'radial', switch ='on'):566 def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, 567 phi, comp='radial', switch='on'): 571 568 """ 572 569 Get the variance from gravity when the wavelength spread is given … … 603 600 return sigma 604 601 605 def _cal_A_value(self, lamda, s_distance, d_distance): 602 def _cal_A_value(self, lamda, s_distance, d_distance): 606 603 """ 607 604 Calculate A value for gravity … … 617 614 gravy = _GRAVITY 618 615 # m/h 619 m_over_h = self.mass / h_constant616 m_over_h = self.mass / h_constant 620 617 # A value 621 618 a_value = d_distance * (s_distance + d_distance) … … 641 638 return self.wave.wavelength 642 639 640 #TODO: why was this method duplicated? 641 #def get_spectrum(self): 642 # """ 643 # Get spectrum 644 # """ 645 # return self.wave.spectrum 646 647 def get_default_spectrum(self): 648 """ 649 Get default_spectrum 650 """ 651 return self.wave.get_default_spectrum() 652 643 653 def get_spectrum(self): 644 654 """ 645 Get spectrum646 """647 return self.wave.spectrum648 649 def get_default_spectrum(self):650 """651 Get default_spectrum652 """653 return self.wave.get_default_spectrum()654 655 def get_spectrum(self):656 """657 655 Get _spectrum 658 656 """ 659 return self.wave.get_spectrum() 657 return self.wave.get_spectrum() 660 658 661 659 def get_wavelength_spread(self): … … 693 691 Get detector size 694 692 """ 695 return self.detector.size 693 return self.detector.size 696 694 697 695 def get_source2sample_distance(self): … … 755 753 """ 756 754 self.spectrum = spectrum 757 self.wave.set_spectrum(spectrum) 755 self.wave.set_spectrum(spectrum) 758 756 759 757 def set_wavelength_spread(self, wavelength_spread): … … 787 785 788 786 : param size: [dia_value] or [x_value, y_value] 789 """ 787 """ 790 788 if len(size) < 1 or len(size) > 2: 791 789 raise RuntimeError, "The length of the size must be one or two." … … 884 882 return qx_min, qx_max, qy_min, qy_max 885 883 886 def _rotate_z(self, x_value, y_value, theta= 884 def _rotate_z(self, x_value, y_value, theta=0.0): 887 885 """ 888 886 Rotate x-y cordinate around z-axis by theta … … 895 893 # rotate by theta 896 894 x_prime = x_value * math.cos(theta) + y_value * math.sin(theta) 897 y_prime = 895 y_prime = -x_value * math.sin(theta) + y_value * math.cos(theta) 898 896 899 897 return x_prime, y_prime 900 898 901 def _gaussian2d(self, x_val, y_val, x0_val, y0_val, 899 def _gaussian2d(self, x_val, y_val, x0_val, y0_val, 902 900 sigma_x, sigma_y, sigma_r): 903 901 """ … … 927 925 y_p = -x_value * sin_phi + y_value * cos_phi 928 926 929 new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x 930 new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y 927 new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x) + 1) 928 new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y) + 1) 931 929 new_x = x_p * cos_phi / new_sig_x - y_p * sin_phi 932 930 new_x /= sigma_x … … 934 932 new_y /= sigma_y 935 933 936 nu_value = -0.5 * (new_x * new_x + new_y * new_y)934 nu_value = -0.5 * (new_x * new_x + new_y * new_y) 937 935 938 936 gaussian = numpy.exp(nu_value) … … 942 940 return gaussian 943 941 944 def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, 942 def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, 945 943 sigma_x, sigma_y, sigma_r): 946 944 """ 947 Calculate 2D Gaussian distribution for polar coodinate 945 Calculate 2D Gaussian distribution for polar coodinate 948 946 : x_val: x value 949 947 : y_val: y value … … 957 955 """ 958 956 sigma_x = sqrt(sigma_x * sigma_x + sigma_r * sigma_r) 959 # call gaussian1d 957 # call gaussian1d 960 958 gaussian = self._gaussian1d(x_val, x0_val, sigma_x) 961 959 gaussian *= self._gaussian1d(y_val, y0_val, sigma_y) … … 971 969 : value: value 972 970 : mean: mean value 973 : sigma: variance 971 : sigma: variance 974 972 975 973 : return: gaussian (value) … … 1028 1026 wavelength = self.wave.wavelength 1029 1027 # Gavity correction 1030 delta_y = self._get_beamcenter_drop() # in cm1028 delta_y = self._get_beamcenter_drop() # in cm 1031 1029 1032 1030 # detector_pix size … … 1042 1040 else: 1043 1041 raise ValueError, " Input value format error..." 1044 # Sample to detector distance = sample slit to detector 1042 # Sample to detector distance = sample slit to detector 1045 1043 # minus sample offset 1046 1044 sample2detector_distance = self.sample2detector_distance[0] - \ … … 1098 1096 1099 1097 # qx_value and qy_value values in array 1100 qx_value = qx_value.repeat(detector_pix_nums_y) 1098 qx_value = qx_value.repeat(detector_pix_nums_y) 1101 1099 qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y) 1102 qy_value = qy_value.repeat(detector_pix_nums_x) 1100 qy_value = qy_value.repeat(detector_pix_nums_x) 1103 1101 qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x) 1104 1102 qy_value = qy_value.transpose() … … 1110 1108 self.qy_max = numpy.max(qy_value) 1111 1109 1112 # Appr. min and max values of the detector display limits 1110 # Appr. min and max values of the detector display limits 1113 1111 # i.e., edges of the last pixels. 1114 self.qy_min += self._get_qx(-0.5 * pix_y_size, 1112 self.qy_min += self._get_qx(-0.5 * pix_y_size, 1115 1113 sample2detector_distance, wavelength) 1116 self.qy_max += self._get_qx(0.5 * pix_y_size, 1114 self.qy_max += self._get_qx(0.5 * pix_y_size, 1117 1115 sample2detector_distance, wavelength) 1118 1116 #if self.qx_min == self.qx_max: 1119 self.qx_min += self._get_qx(-0.5 * pix_x_size, 1117 self.qx_min += self._get_qx(-0.5 * pix_x_size, 1120 1118 sample2detector_distance, wavelength) 1121 self.qx_max += self._get_qx(0.5 * pix_x_size, 1119 self.qx_max += self._get_qx(0.5 * pix_x_size, 1122 1120 sample2detector_distance, wavelength) 1123 1121 … … 1128 1126 self.detector_qy_max = self.qy_max 1129 1127 1130 1131 1132 1128 # try to set it as a Data2D otherwise pass (not required for now) 1133 1129 try: … … 1135 1131 output = Data2D() 1136 1132 inten = numpy.zeros_like(qx_value) 1137 output.data 1138 output.qx_data 1139 output.qy_data = qy_value1133 output.data = inten 1134 output.qx_data = qx_value 1135 output.qy_data = qy_value 1140 1136 except: 1141 1137 pass 1142 1138 1143 return output #qx_value,qy_value1139 return output 1144 1140 1145 1141 def _get_qx(self, dx_size, det_dist, wavelength): … … 1153 1149 plane_dist = dx_size 1154 1150 # full scattering angle on the x-axis 1155 theta 1156 qx_value 1157 return qx_value 1151 theta = numpy.arctan(plane_dist / det_dist) 1152 qx_value = (2.0 * pi / wavelength) * numpy.sin(theta) 1153 return qx_value 1158 1154 1159 1155 def _get_polar_value(self, qx_value, qy_value): … … 1187 1183 pos_y -= offset_y 1188 1184 1189 return pos_x, pos_y 1185 return pos_x, pos_y 1190 1186 1191 1187 def _get_beamcenter_drop(self): … … 1195 1191 :return delta y: the beam center drop in cm 1196 1192 """ 1197 # Check if mass == 0 (X-ray). 1193 # Check if mass == 0 (X-ray). 1198 1194 if self.mass == 0: 1199 1195 return 0 … … 1211 1207 delta_y /= (velocity * velocity) 1212 1208 1213 return delta_y 1214 1215 1209 return delta_y
Note: See TracChangeset
for help on using the changeset viewer.
