| 1 | """ |
|---|
| 2 | This object is a small tool to allow user to quickly |
|---|
| 3 | determine the variance in q from the |
|---|
| 4 | instrumental parameters. |
|---|
| 5 | """ |
|---|
| 6 | from instrument import Sample |
|---|
| 7 | from instrument import Detector |
|---|
| 8 | from instrument import TOF as Neutron |
|---|
| 9 | from instrument import Aperture |
|---|
| 10 | # import math stuffs |
|---|
| 11 | from math import pi |
|---|
| 12 | from math import sqrt |
|---|
| 13 | import math |
|---|
| 14 | import numpy as np |
|---|
| 15 | import sys |
|---|
| 16 | import logging |
|---|
| 17 | |
|---|
| 18 | logger = logging.getLogger(__name__) |
|---|
| 19 | |
|---|
| 20 | #Plank's constant in cgs unit |
|---|
| 21 | _PLANK_H = 6.62606896E-27 |
|---|
| 22 | #Gravitational acc. in cgs unit |
|---|
| 23 | _GRAVITY = 981.0 |
|---|
| 24 | |
|---|
| 25 | |
|---|
| 26 | class ResolutionCalculator(object): |
|---|
| 27 | """ |
|---|
| 28 | compute resolution in 2D |
|---|
| 29 | """ |
|---|
| 30 | def __init__(self): |
|---|
| 31 | |
|---|
| 32 | # wavelength |
|---|
| 33 | self.wave = Neutron() |
|---|
| 34 | # sample |
|---|
| 35 | self.sample = Sample() |
|---|
| 36 | # aperture |
|---|
| 37 | self.aperture = Aperture() |
|---|
| 38 | # detector |
|---|
| 39 | self.detector = Detector() |
|---|
| 40 | # 2d image of the resolution |
|---|
| 41 | self.image = [] |
|---|
| 42 | self.image_lam = [] |
|---|
| 43 | # resolutions |
|---|
| 44 | # lamda in r-direction |
|---|
| 45 | self.sigma_lamd = 0 |
|---|
| 46 | # x-dir (no lamda) |
|---|
| 47 | self.sigma_1 = 0 |
|---|
| 48 | #y-dir (no lamda) |
|---|
| 49 | self.sigma_2 = 0 |
|---|
| 50 | # 1D total |
|---|
| 51 | self.sigma_1d = 0 |
|---|
| 52 | self.gravity_phi = None |
|---|
| 53 | # q min and max |
|---|
| 54 | self.qx_min = -0.3 |
|---|
| 55 | self.qx_max = 0.3 |
|---|
| 56 | self.qy_min = -0.3 |
|---|
| 57 | self.qy_max = 0.3 |
|---|
| 58 | # q min and max of the detector |
|---|
| 59 | self.detector_qx_min = -0.3 |
|---|
| 60 | self.detector_qx_max = 0.3 |
|---|
| 61 | self.detector_qy_min = -0.3 |
|---|
| 62 | self.detector_qy_max = 0.3 |
|---|
| 63 | # possible max qrange |
|---|
| 64 | self.qxmin_limit = 0 |
|---|
| 65 | self.qxmax_limit = 0 |
|---|
| 66 | self.qymin_limit = 0 |
|---|
| 67 | self.qymax_limit = 0 |
|---|
| 68 | |
|---|
| 69 | # plots |
|---|
| 70 | self.plot = None |
|---|
| 71 | # instrumental params defaults |
|---|
| 72 | self.mass = 0 |
|---|
| 73 | self.intensity = 0 |
|---|
| 74 | self.wavelength = 0 |
|---|
| 75 | self.wavelength_spread = 0 |
|---|
| 76 | self.source_aperture_size = [] |
|---|
| 77 | self.source2sample_distance = [] |
|---|
| 78 | self.sample2sample_distance = [] |
|---|
| 79 | self.sample_aperture_size = [] |
|---|
| 80 | self.sample2detector_distance = [] |
|---|
| 81 | self.detector_pix_size = [] |
|---|
| 82 | self.detector_size = [] |
|---|
| 83 | self.get_all_instrument_params() |
|---|
| 84 | # max q range for all lambdas |
|---|
| 85 | self.qxrange = [] |
|---|
| 86 | self.qyrange = [] |
|---|
| 87 | |
|---|
| 88 | def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, |
|---|
| 89 | qy_min, qy_max, coord='cartesian'): |
|---|
| 90 | """ |
|---|
| 91 | Compute the resolution |
|---|
| 92 | : qx_value: x component of q |
|---|
| 93 | : qy_value: y component of q |
|---|
| 94 | """ |
|---|
| 95 | |
|---|
| 96 | # make sure to update all the variables need. |
|---|
| 97 | # except lambda, dlambda, and intensity |
|---|
| 98 | self.get_all_instrument_params() |
|---|
| 99 | # wavelength etc. |
|---|
| 100 | lamda_list, dlamb_list = self.get_wave_list() |
|---|
| 101 | intens_list = [] |
|---|
| 102 | sig1_list = [] |
|---|
| 103 | sig2_list = [] |
|---|
| 104 | sigr_list = [] |
|---|
| 105 | sigma1d_list = [] |
|---|
| 106 | num_lamda = len(lamda_list) |
|---|
| 107 | for num in range(num_lamda): |
|---|
| 108 | lam = lamda_list[num] |
|---|
| 109 | # wavelength spread |
|---|
| 110 | dlam = dlamb_list[num] |
|---|
| 111 | intens = self.setup_tof(lam, dlam) |
|---|
| 112 | intens_list.append(intens) |
|---|
| 113 | # cehck if tof |
|---|
| 114 | if num_lamda > 1: |
|---|
| 115 | tof = True |
|---|
| 116 | else: |
|---|
| 117 | tof = False |
|---|
| 118 | # compute 2d resolution |
|---|
| 119 | _, _, sigma_1, sigma_2, sigma_r, sigma1d = \ |
|---|
| 120 | self.compute(lam, dlam, qx_value, qy_value, coord, tof) |
|---|
| 121 | # make image |
|---|
| 122 | image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, |
|---|
| 123 | sigma_r, qx_min, qx_max, qy_min, qy_max, |
|---|
| 124 | coord, False) |
|---|
| 125 | if qx_min > self.qx_min: |
|---|
| 126 | qx_min = self.qx_min |
|---|
| 127 | if qx_max < self.qx_max: |
|---|
| 128 | qx_max = self.qx_max |
|---|
| 129 | if qy_min > self.qy_min: |
|---|
| 130 | qy_min = self.qy_min |
|---|
| 131 | if qy_max < self.qy_max: |
|---|
| 132 | qy_max = self.qy_max |
|---|
| 133 | |
|---|
| 134 | # set max qranges |
|---|
| 135 | self.qxrange = [qx_min, qx_max] |
|---|
| 136 | self.qyrange = [qy_min, qy_max] |
|---|
| 137 | sig1_list.append(sigma_1) |
|---|
| 138 | sig2_list.append(sigma_2) |
|---|
| 139 | sigr_list.append(sigma_r) |
|---|
| 140 | sigma1d_list.append(sigma1d) |
|---|
| 141 | # redraw image in global 2d q-space. |
|---|
| 142 | self.image_lam = [] |
|---|
| 143 | total_intensity = 0 |
|---|
| 144 | sigma_1 = 0 |
|---|
| 145 | sigma_r = 0 |
|---|
| 146 | sigma_2 = 0 |
|---|
| 147 | sigma1d = 0 |
|---|
| 148 | for ind in range(num_lamda): |
|---|
| 149 | lam = lamda_list[ind] |
|---|
| 150 | dlam = dlamb_list[ind] |
|---|
| 151 | intens = self.setup_tof(lam, dlam) |
|---|
| 152 | out = self.get_image(qx_value, qy_value, sig1_list[ind], |
|---|
| 153 | sig2_list[ind], sigr_list[ind], |
|---|
| 154 | qx_min, qx_max, qy_min, qy_max, coord) |
|---|
| 155 | image = out |
|---|
| 156 | # set variance as sigmas |
|---|
| 157 | sigma_1 += sig1_list[ind] * sig1_list[ind] * self.intensity |
|---|
| 158 | sigma_r += sigr_list[ind] * sigr_list[ind] * self.intensity |
|---|
| 159 | sigma_2 += sig2_list[ind] * sig2_list[ind] * self.intensity |
|---|
| 160 | sigma1d += sigma1d_list[ind] * sigma1d_list[ind] * self.intensity |
|---|
| 161 | total_intensity += self.intensity |
|---|
| 162 | |
|---|
| 163 | if total_intensity != 0: |
|---|
| 164 | # average variance |
|---|
| 165 | image_out = image / total_intensity |
|---|
| 166 | sigma_1 = sigma_1 / total_intensity |
|---|
| 167 | sigma_r = sigma_r / total_intensity |
|---|
| 168 | sigma_2 = sigma_2 / total_intensity |
|---|
| 169 | sigma1d = sigma1d / total_intensity |
|---|
| 170 | # set sigmas |
|---|
| 171 | self.sigma_1 = sqrt(sigma_1) |
|---|
| 172 | self.sigma_lamd = sqrt(sigma_r) |
|---|
| 173 | self.sigma_2 = sqrt(sigma_2) |
|---|
| 174 | self.sigma_1d = sqrt(sigma1d) |
|---|
| 175 | # rescale |
|---|
| 176 | max_im_val = 1 |
|---|
| 177 | if max_im_val > 0: |
|---|
| 178 | image_out /= max_im_val |
|---|
| 179 | else: |
|---|
| 180 | image_out = image * 0.0 |
|---|
| 181 | # Don't calculate sigmas nor set self.sigmas! |
|---|
| 182 | sigma_1 = 0 |
|---|
| 183 | sigma_r = 0 |
|---|
| 184 | sigma_2 = 0 |
|---|
| 185 | sigma1d = 0 |
|---|
| 186 | |
|---|
| 187 | if len(self.image) > 0: |
|---|
| 188 | self.image += image_out |
|---|
| 189 | else: |
|---|
| 190 | self.image = image_out |
|---|
| 191 | |
|---|
| 192 | # plot image |
|---|
| 193 | return self.image # self.plot_image(self.image) |
|---|
| 194 | |
|---|
| 195 | def setup_tof(self, wavelength, wavelength_spread): |
|---|
| 196 | """ |
|---|
| 197 | Setup all parameters in instrument |
|---|
| 198 | |
|---|
| 199 | : param ind: index of lambda, etc |
|---|
| 200 | """ |
|---|
| 201 | |
|---|
| 202 | # set wave.wavelength |
|---|
| 203 | self.set_wavelength(wavelength) |
|---|
| 204 | self.set_wavelength_spread(wavelength_spread) |
|---|
| 205 | self.intensity = self.wave.get_intensity() |
|---|
| 206 | |
|---|
| 207 | if wavelength == 0: |
|---|
| 208 | msg = "Can't compute the resolution: the wavelength is zero..." |
|---|
| 209 | raise RuntimeError, msg |
|---|
| 210 | return self.intensity |
|---|
| 211 | |
|---|
| 212 | def compute(self, wavelength, wavelength_spread, qx_value, qy_value, |
|---|
| 213 | coord='cartesian', tof=False): |
|---|
| 214 | """ |
|---|
| 215 | Compute the Q resolution in || and + direction of 2D |
|---|
| 216 | : qx_value: x component of q |
|---|
| 217 | : qy_value: y component of q |
|---|
| 218 | """ |
|---|
| 219 | coord = 'cartesian' |
|---|
| 220 | lamb = wavelength |
|---|
| 221 | lamb_spread = wavelength_spread |
|---|
| 222 | # the shape of wavelength distribution |
|---|
| 223 | |
|---|
| 224 | if tof: |
|---|
| 225 | # rectangular |
|---|
| 226 | tof_factor = 2 |
|---|
| 227 | else: |
|---|
| 228 | # triangular |
|---|
| 229 | tof_factor = 1 |
|---|
| 230 | # Find polar values |
|---|
| 231 | qr_value, phi = self._get_polar_value(qx_value, qy_value) |
|---|
| 232 | # vacuum wave transfer |
|---|
| 233 | knot = 2*pi/lamb |
|---|
| 234 | # scattering angle theta; always true for plane detector |
|---|
| 235 | # aligned vertically to the ko direction |
|---|
| 236 | if qr_value > knot: |
|---|
| 237 | theta = pi/2 |
|---|
| 238 | else: |
|---|
| 239 | theta = math.asin(qr_value/knot) |
|---|
| 240 | # source aperture size |
|---|
| 241 | rone = self.source_aperture_size |
|---|
| 242 | # sample aperture size |
|---|
| 243 | rtwo = self.sample_aperture_size |
|---|
| 244 | # detector pixel size |
|---|
| 245 | rthree = self.detector_pix_size |
|---|
| 246 | # source to sample(aperture) distance |
|---|
| 247 | l_ssa = self.source2sample_distance[0] |
|---|
| 248 | # sample(aperture) to detector distance |
|---|
| 249 | l_sad = self.sample2detector_distance[0] |
|---|
| 250 | # sample (aperture) to sample distance |
|---|
| 251 | l_sas = self.sample2sample_distance[0] |
|---|
| 252 | # source to sample distance |
|---|
| 253 | l_one = l_ssa + l_sas |
|---|
| 254 | # sample to detector distance |
|---|
| 255 | l_two = l_sad - l_sas |
|---|
| 256 | |
|---|
| 257 | # Sample offset correction for l_one and Lp on variance calculation |
|---|
| 258 | l1_cor = (l_ssa * l_two) / (l_sas + l_two) |
|---|
| 259 | lp_cor = (l_ssa * l_two) / (l_one + l_two) |
|---|
| 260 | # the radial distance to the pixel from the center of the detector |
|---|
| 261 | radius = math.tan(theta) * l_two |
|---|
| 262 | #Lp = l_one*l_two/(l_one+l_two) |
|---|
| 263 | # default polar coordinate |
|---|
| 264 | comp1 = 'radial' |
|---|
| 265 | comp2 = 'phi' |
|---|
| 266 | # in the case of the cartesian coordinate |
|---|
| 267 | if coord == 'cartesian': |
|---|
| 268 | comp1 = 'x' |
|---|
| 269 | comp2 = 'y' |
|---|
| 270 | |
|---|
| 271 | # sigma in the radial/x direction |
|---|
| 272 | # for source aperture |
|---|
| 273 | sigma_1 = self.get_variance(rone, l1_cor, phi, comp1) |
|---|
| 274 | # for sample apperture |
|---|
| 275 | sigma_1 += self.get_variance(rtwo, lp_cor, phi, comp1) |
|---|
| 276 | # for detector pix |
|---|
| 277 | sigma_1 += self.get_variance(rthree, l_two, phi, comp1) |
|---|
| 278 | # for gravity term for 1d |
|---|
| 279 | sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, |
|---|
| 280 | lamb_spread, phi, comp1, 'on') / tof_factor |
|---|
| 281 | # for wavelength spread |
|---|
| 282 | # reserve for 1d calculation |
|---|
| 283 | A_value = self._cal_A_value(lamb, l_ssa, l_sad) |
|---|
| 284 | sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value, |
|---|
| 285 | radius, l_two, lamb_spread, |
|---|
| 286 | phi, 'radial', 'on') |
|---|
| 287 | sigma_wave_1 /= tof_factor |
|---|
| 288 | sigma_wave_1_1d /= tof_factor |
|---|
| 289 | # for 1d |
|---|
| 290 | variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d |
|---|
| 291 | # normalize |
|---|
| 292 | variance_1d_1 = knot * knot * variance_1d_1 / 12 |
|---|
| 293 | |
|---|
| 294 | # for 2d |
|---|
| 295 | #sigma_1 += sigma_wave_1 |
|---|
| 296 | # normalize |
|---|
| 297 | sigma_1 = knot * sqrt(sigma_1 / 12) |
|---|
| 298 | sigma_r = knot * sqrt(sigma_wave_1 / (tof_factor *12)) |
|---|
| 299 | # sigma in the phi/y direction |
|---|
| 300 | # for source apperture |
|---|
| 301 | sigma_2 = self.get_variance(rone, l1_cor, phi, comp2) |
|---|
| 302 | |
|---|
| 303 | # for sample apperture |
|---|
| 304 | sigma_2 += self.get_variance(rtwo, lp_cor, phi, comp2) |
|---|
| 305 | |
|---|
| 306 | # for detector pix |
|---|
| 307 | sigma_2 += self.get_variance(rthree, l_two, phi, comp2) |
|---|
| 308 | |
|---|
| 309 | # for gravity term for 1d |
|---|
| 310 | sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, |
|---|
| 311 | lamb_spread, phi, comp2, 'on') / tof_factor |
|---|
| 312 | |
|---|
| 313 | # for wavelength spread |
|---|
| 314 | # reserve for 1d calculation |
|---|
| 315 | sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value, |
|---|
| 316 | radius, l_two, lamb_spread, |
|---|
| 317 | phi, 'phi', 'on') |
|---|
| 318 | sigma_wave_2 /= tof_factor |
|---|
| 319 | sigma_wave_2_1d /= tof_factor |
|---|
| 320 | # for 1d |
|---|
| 321 | variance_1d_2 = (sigma_2 + sigma_2grav1d) / 2 + sigma_wave_2_1d |
|---|
| 322 | # normalize |
|---|
| 323 | variance_1d_2 = knot * knot * variance_1d_2 / 12 |
|---|
| 324 | |
|---|
| 325 | # for 2d |
|---|
| 326 | #sigma_2 = knot*sqrt(sigma_2/12) |
|---|
| 327 | #sigma_2 += sigma_wave_2 |
|---|
| 328 | # normalize |
|---|
| 329 | sigma_2 = knot * sqrt(sigma_2 / 12) |
|---|
| 330 | sigma1d = sqrt(variance_1d_1 + variance_1d_2) |
|---|
| 331 | # set sigmas |
|---|
| 332 | self.sigma_1 = sigma_1 |
|---|
| 333 | self.sigma_lamd = sigma_r |
|---|
| 334 | self.sigma_2 = sigma_2 |
|---|
| 335 | self.sigma_1d = sigma1d |
|---|
| 336 | return qr_value, phi, sigma_1, sigma_2, sigma_r, sigma1d |
|---|
| 337 | |
|---|
| 338 | def _within_detector_range(self, qx_value, qy_value): |
|---|
| 339 | """ |
|---|
| 340 | check if qvalues are within detector range |
|---|
| 341 | """ |
|---|
| 342 | # detector range |
|---|
| 343 | detector_qx_min = self.detector_qx_min |
|---|
| 344 | detector_qx_max = self.detector_qx_max |
|---|
| 345 | detector_qy_min = self.detector_qy_min |
|---|
| 346 | detector_qy_max = self.detector_qy_max |
|---|
| 347 | if self.qxmin_limit > detector_qx_min: |
|---|
| 348 | self.qxmin_limit = detector_qx_min |
|---|
| 349 | if self.qxmax_limit < detector_qx_max: |
|---|
| 350 | self.qxmax_limit = detector_qx_max |
|---|
| 351 | if self.qymin_limit > detector_qy_min: |
|---|
| 352 | self.qymin_limit = detector_qy_min |
|---|
| 353 | if self.qymax_limit < detector_qy_max: |
|---|
| 354 | self.qymax_limit = detector_qy_max |
|---|
| 355 | if qx_value < detector_qx_min or qx_value > detector_qx_max: |
|---|
| 356 | return False |
|---|
| 357 | if qy_value < detector_qy_min or qy_value > detector_qy_max: |
|---|
| 358 | return False |
|---|
| 359 | return True |
|---|
| 360 | |
|---|
| 361 | def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r, |
|---|
| 362 | qx_min, qx_max, qy_min, qy_max, |
|---|
| 363 | coord='cartesian', full_cal=True): |
|---|
| 364 | """ |
|---|
| 365 | Get the resolution in polar coordinate ready to plot |
|---|
| 366 | : qx_value: qx_value value |
|---|
| 367 | : qy_value: qy_value value |
|---|
| 368 | : sigma_1: variance in r direction |
|---|
| 369 | : sigma_2: variance in phi direction |
|---|
| 370 | : coord: coordinate system of image, 'polar' or 'cartesian' |
|---|
| 371 | """ |
|---|
| 372 | # Get qx_max and qy_max... |
|---|
| 373 | self._get_detector_qxqy_pixels() |
|---|
| 374 | |
|---|
| 375 | qr_value, phi = self._get_polar_value(qx_value, qy_value) |
|---|
| 376 | |
|---|
| 377 | # Check whether the q value is within the detector range |
|---|
| 378 | if qx_min < self.qx_min: |
|---|
| 379 | self.qx_min = qx_min |
|---|
| 380 | #raise ValueError, msg |
|---|
| 381 | if qx_max > self.qx_max: |
|---|
| 382 | self.qx_max = qx_max |
|---|
| 383 | #raise ValueError, msg |
|---|
| 384 | if qy_min < self.qy_min: |
|---|
| 385 | self.qy_min = qy_min |
|---|
| 386 | #raise ValueError, msg |
|---|
| 387 | if qy_max > self.qy_max: |
|---|
| 388 | self.qy_max = qy_max |
|---|
| 389 | #raise ValueError, msg |
|---|
| 390 | if not full_cal: |
|---|
| 391 | return None |
|---|
| 392 | |
|---|
| 393 | # Make an empty graph in the detector scale |
|---|
| 394 | dx_size = (self.qx_max - self.qx_min) / (1000 - 1) |
|---|
| 395 | dy_size = (self.qy_max - self.qy_min) / (1000 - 1) |
|---|
| 396 | x_val = np.arange(self.qx_min, self.qx_max, dx_size) |
|---|
| 397 | y_val = np.arange(self.qy_max, self.qy_min, -dy_size) |
|---|
| 398 | q_1, q_2 = np.meshgrid(x_val, y_val) |
|---|
| 399 | #q_phi = numpy.arctan(q_1,q_2) |
|---|
| 400 | # check whether polar or cartesian |
|---|
| 401 | if coord == 'polar': |
|---|
| 402 | # Find polar values |
|---|
| 403 | qr_value, phi = self._get_polar_value(qx_value, qy_value) |
|---|
| 404 | q_1, q_2 = self._rotate_z(q_1, q_2, phi) |
|---|
| 405 | qc_1 = qr_value |
|---|
| 406 | qc_2 = 0.0 |
|---|
| 407 | # Calculate the 2D Gaussian distribution image |
|---|
| 408 | image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, |
|---|
| 409 | sigma_1, sigma_2, sigma_r) |
|---|
| 410 | else: |
|---|
| 411 | # catesian coordinate |
|---|
| 412 | # qx_center |
|---|
| 413 | qc_1 = qx_value |
|---|
| 414 | # qy_center |
|---|
| 415 | qc_2 = qy_value |
|---|
| 416 | |
|---|
| 417 | # Calculate the 2D Gaussian distribution image |
|---|
| 418 | image = self._gaussian2d(q_1, q_2, qc_1, qc_2, |
|---|
| 419 | sigma_1, sigma_2, sigma_r) |
|---|
| 420 | # out side of detector |
|---|
| 421 | if not self._within_detector_range(qx_value, qy_value): |
|---|
| 422 | image *= 0.0 |
|---|
| 423 | self.intensity = 0.0 |
|---|
| 424 | #return self.image |
|---|
| 425 | |
|---|
| 426 | # Add it if there are more than one inputs. |
|---|
| 427 | if len(self.image_lam) > 0: |
|---|
| 428 | self.image_lam += image * self.intensity |
|---|
| 429 | else: |
|---|
| 430 | self.image_lam = image * self.intensity |
|---|
| 431 | |
|---|
| 432 | return self.image_lam |
|---|
| 433 | |
|---|
| 434 | def plot_image(self, image): |
|---|
| 435 | """ |
|---|
| 436 | Plot image using pyplot |
|---|
| 437 | : image: 2d resolution image |
|---|
| 438 | |
|---|
| 439 | : return plt: pylab object |
|---|
| 440 | """ |
|---|
| 441 | import matplotlib.pyplot as plt |
|---|
| 442 | |
|---|
| 443 | self.plot = plt |
|---|
| 444 | plt.xlabel('$\\rm{Q}_{x} [A^{-1}]$') |
|---|
| 445 | plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$') |
|---|
| 446 | # Max value of the image |
|---|
| 447 | # max = numpy.max(image) |
|---|
| 448 | qx_min, qx_max, qy_min, qy_max = self.get_detector_qrange() |
|---|
| 449 | |
|---|
| 450 | # Image |
|---|
| 451 | im = plt.imshow(image, |
|---|
| 452 | extent=[qx_min, qx_max, qy_min, qy_max]) |
|---|
| 453 | |
|---|
| 454 | # bilinear interpolation to make it smoother |
|---|
| 455 | im.set_interpolation('bilinear') |
|---|
| 456 | |
|---|
| 457 | return plt |
|---|
| 458 | |
|---|
| 459 | def reset_image(self): |
|---|
| 460 | """ |
|---|
| 461 | Reset image to default (=[]) |
|---|
| 462 | """ |
|---|
| 463 | self.image = [] |
|---|
| 464 | |
|---|
| 465 | def get_variance(self, size=[], distance=0, phi=0, comp='radial'): |
|---|
| 466 | """ |
|---|
| 467 | Get the variance when the slit/pinhole size is given |
|---|
| 468 | : size: list that can be one(diameter for circular) or two components(lengths for rectangular) |
|---|
| 469 | : distance: [z, x] where z along the incident beam, x // qx_value |
|---|
| 470 | : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial' |
|---|
| 471 | |
|---|
| 472 | : return variance: sigma^2 |
|---|
| 473 | """ |
|---|
| 474 | # check the length of size (list) |
|---|
| 475 | len_size = len(size) |
|---|
| 476 | |
|---|
| 477 | # define sigma component direction |
|---|
| 478 | if comp == 'radial': |
|---|
| 479 | phi_x = math.cos(phi) |
|---|
| 480 | phi_y = math.sin(phi) |
|---|
| 481 | elif comp == 'phi': |
|---|
| 482 | phi_x = math.sin(phi) |
|---|
| 483 | phi_y = math.cos(phi) |
|---|
| 484 | elif comp == 'x': |
|---|
| 485 | phi_x = 1 |
|---|
| 486 | phi_y = 0 |
|---|
| 487 | elif comp == 'y': |
|---|
| 488 | phi_x = 0 |
|---|
| 489 | phi_y = 1 |
|---|
| 490 | else: |
|---|
| 491 | phi_x = 0 |
|---|
| 492 | phi_y = 0 |
|---|
| 493 | # calculate each component |
|---|
| 494 | # for pinhole w/ radius = size[0]/2 |
|---|
| 495 | if len_size == 1: |
|---|
| 496 | x_comp = (0.5 * size[0]) * sqrt(3) |
|---|
| 497 | y_comp = 0 |
|---|
| 498 | # for rectangular slit |
|---|
| 499 | elif len_size == 2: |
|---|
| 500 | x_comp = size[0] * phi_x |
|---|
| 501 | y_comp = size[1] * phi_y |
|---|
| 502 | # otherwise |
|---|
| 503 | else: |
|---|
| 504 | raise ValueError, " Improper input..." |
|---|
| 505 | # get them squared |
|---|
| 506 | sigma = x_comp * x_comp |
|---|
| 507 | sigma += y_comp * y_comp |
|---|
| 508 | # normalize by distance |
|---|
| 509 | sigma /= (distance * distance) |
|---|
| 510 | |
|---|
| 511 | return sigma |
|---|
| 512 | |
|---|
| 513 | def get_variance_wave(self, A_value, radius, distance, spread, phi, |
|---|
| 514 | comp='radial', switch='on'): |
|---|
| 515 | """ |
|---|
| 516 | Get the variance when the wavelength spread is given |
|---|
| 517 | |
|---|
| 518 | : radius: the radial distance from the beam center to the pix of q |
|---|
| 519 | : distance: sample to detector distance |
|---|
| 520 | : spread: wavelength spread (ratio) |
|---|
| 521 | : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial' |
|---|
| 522 | |
|---|
| 523 | : return variance: sigma^2 for 2d, sigma^2 for 1d [tuple] |
|---|
| 524 | """ |
|---|
| 525 | if switch.lower() == 'off': |
|---|
| 526 | return 0, 0 |
|---|
| 527 | # check the singular point |
|---|
| 528 | if distance == 0 or comp == 'phi': |
|---|
| 529 | return 0, 0 |
|---|
| 530 | else: |
|---|
| 531 | # calculate sigma^2 for 1d |
|---|
| 532 | sigma1d = 2 * math.pow(radius/distance*spread, 2) |
|---|
| 533 | if comp == 'x': |
|---|
| 534 | sigma1d *= (math.cos(phi)*math.cos(phi)) |
|---|
| 535 | elif comp == 'y': |
|---|
| 536 | sigma1d *= (math.sin(phi)*math.sin(phi)) |
|---|
| 537 | else: |
|---|
| 538 | sigma1d *= 1 |
|---|
| 539 | # sigma^2 for 2d |
|---|
| 540 | # shift the coordinate due to the gravitational shift |
|---|
| 541 | rad_x = radius * math.cos(phi) |
|---|
| 542 | rad_y = A_value - radius * math.sin(phi) |
|---|
| 543 | radius = math.sqrt(rad_x * rad_x + rad_y * rad_y) |
|---|
| 544 | # new phi |
|---|
| 545 | phi = math.atan2(-rad_y, rad_x) |
|---|
| 546 | self.gravity_phi = phi |
|---|
| 547 | # calculate sigma^2 |
|---|
| 548 | sigma = 2 * math.pow(radius/distance*spread, 2) |
|---|
| 549 | if comp == 'x': |
|---|
| 550 | sigma *= (math.cos(phi)*math.cos(phi)) |
|---|
| 551 | elif comp == 'y': |
|---|
| 552 | sigma *= (math.sin(phi)*math.sin(phi)) |
|---|
| 553 | else: |
|---|
| 554 | sigma *= 1 |
|---|
| 555 | |
|---|
| 556 | return sigma, sigma1d |
|---|
| 557 | |
|---|
| 558 | def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, |
|---|
| 559 | phi, comp='radial', switch='on'): |
|---|
| 560 | """ |
|---|
| 561 | Get the variance from gravity when the wavelength spread is given |
|---|
| 562 | |
|---|
| 563 | : s_distance: source to sample distance |
|---|
| 564 | : d_distance: sample to detector distance |
|---|
| 565 | : wavelength: wavelength |
|---|
| 566 | : spread: wavelength spread (ratio) |
|---|
| 567 | : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial' |
|---|
| 568 | |
|---|
| 569 | : return variance: sigma^2 |
|---|
| 570 | """ |
|---|
| 571 | if switch.lower() == 'off': |
|---|
| 572 | return 0 |
|---|
| 573 | if self.mass == 0.0: |
|---|
| 574 | return 0 |
|---|
| 575 | # check the singular point |
|---|
| 576 | if d_distance == 0 or comp == 'x': |
|---|
| 577 | return 0 |
|---|
| 578 | else: |
|---|
| 579 | a_value = self._cal_A_value(None, s_distance, d_distance) |
|---|
| 580 | # calculate sigma^2 |
|---|
| 581 | sigma = math.pow(a_value / d_distance, 2) |
|---|
| 582 | sigma *= math.pow(wavelength, 4) |
|---|
| 583 | sigma *= math.pow(spread, 2) |
|---|
| 584 | sigma *= 8 |
|---|
| 585 | return sigma |
|---|
| 586 | |
|---|
| 587 | def _cal_A_value(self, lamda, s_distance, d_distance): |
|---|
| 588 | """ |
|---|
| 589 | Calculate A value for gravity |
|---|
| 590 | |
|---|
| 591 | : s_distance: source to sample distance |
|---|
| 592 | : d_distance: sample to detector distance |
|---|
| 593 | """ |
|---|
| 594 | # neutron mass in cgs unit |
|---|
| 595 | self.mass = self.get_neutron_mass() |
|---|
| 596 | # plank constant in cgs unit |
|---|
| 597 | h_constant = _PLANK_H |
|---|
| 598 | # gravity in cgs unit |
|---|
| 599 | gravy = _GRAVITY |
|---|
| 600 | # m/h |
|---|
| 601 | m_over_h = self.mass / h_constant |
|---|
| 602 | # A value |
|---|
| 603 | a_value = d_distance * (s_distance + d_distance) |
|---|
| 604 | a_value *= math.pow(m_over_h / 2, 2) |
|---|
| 605 | a_value *= gravy |
|---|
| 606 | # unit correction (1/cm to 1/A) for A and d_distance below |
|---|
| 607 | a_value *= 1.0E-16 |
|---|
| 608 | # if lamda is give (broad meanning of A) return 2* lamda^2 * A |
|---|
| 609 | if lamda is not None: |
|---|
| 610 | a_value *= (4 * lamda * lamda) |
|---|
| 611 | return a_value |
|---|
| 612 | |
|---|
| 613 | def get_intensity(self): |
|---|
| 614 | """ |
|---|
| 615 | Get intensity |
|---|
| 616 | """ |
|---|
| 617 | return self.wave.intensity |
|---|
| 618 | |
|---|
| 619 | def get_wavelength(self): |
|---|
| 620 | """ |
|---|
| 621 | Get wavelength |
|---|
| 622 | """ |
|---|
| 623 | return self.wave.wavelength |
|---|
| 624 | |
|---|
| 625 | def get_default_spectrum(self): |
|---|
| 626 | """ |
|---|
| 627 | Get default_spectrum |
|---|
| 628 | """ |
|---|
| 629 | return self.wave.get_default_spectrum() |
|---|
| 630 | |
|---|
| 631 | def get_spectrum(self): |
|---|
| 632 | """ |
|---|
| 633 | Get _spectrum |
|---|
| 634 | """ |
|---|
| 635 | return self.wave.get_spectrum() |
|---|
| 636 | |
|---|
| 637 | def get_wavelength_spread(self): |
|---|
| 638 | """ |
|---|
| 639 | Get wavelength spread |
|---|
| 640 | """ |
|---|
| 641 | return self.wave.wavelength_spread |
|---|
| 642 | |
|---|
| 643 | def get_neutron_mass(self): |
|---|
| 644 | """ |
|---|
| 645 | Get Neutron mass |
|---|
| 646 | """ |
|---|
| 647 | return self.wave.mass |
|---|
| 648 | |
|---|
| 649 | def get_source_aperture_size(self): |
|---|
| 650 | """ |
|---|
| 651 | Get source aperture size |
|---|
| 652 | """ |
|---|
| 653 | return self.aperture.source_size |
|---|
| 654 | |
|---|
| 655 | def get_sample_aperture_size(self): |
|---|
| 656 | """ |
|---|
| 657 | Get sample aperture size |
|---|
| 658 | """ |
|---|
| 659 | return self.aperture.sample_size |
|---|
| 660 | |
|---|
| 661 | def get_detector_pix_size(self): |
|---|
| 662 | """ |
|---|
| 663 | Get detector pixel size |
|---|
| 664 | """ |
|---|
| 665 | return self.detector.pix_size |
|---|
| 666 | |
|---|
| 667 | def get_detector_size(self): |
|---|
| 668 | """ |
|---|
| 669 | Get detector size |
|---|
| 670 | """ |
|---|
| 671 | return self.detector.size |
|---|
| 672 | |
|---|
| 673 | def get_source2sample_distance(self): |
|---|
| 674 | """ |
|---|
| 675 | Get detector source2sample_distance |
|---|
| 676 | """ |
|---|
| 677 | return self.aperture.sample_distance |
|---|
| 678 | |
|---|
| 679 | def get_sample2sample_distance(self): |
|---|
| 680 | """ |
|---|
| 681 | Get detector sampleslitsample_distance |
|---|
| 682 | """ |
|---|
| 683 | return self.sample.distance |
|---|
| 684 | |
|---|
| 685 | def get_sample2detector_distance(self): |
|---|
| 686 | """ |
|---|
| 687 | Get detector sample2detector_distance |
|---|
| 688 | """ |
|---|
| 689 | return self.detector.distance |
|---|
| 690 | |
|---|
| 691 | def set_intensity(self, intensity): |
|---|
| 692 | """ |
|---|
| 693 | Set intensity |
|---|
| 694 | """ |
|---|
| 695 | self.wave.set_intensity(intensity) |
|---|
| 696 | |
|---|
| 697 | def set_wave(self, wavelength): |
|---|
| 698 | """ |
|---|
| 699 | Set wavelength list or wavelength |
|---|
| 700 | """ |
|---|
| 701 | if wavelength.__class__.__name__ == 'list': |
|---|
| 702 | self.wave.set_wave_list(wavelength) |
|---|
| 703 | elif wavelength.__class__.__name__ == 'float': |
|---|
| 704 | self.wave.set_wave_list([wavelength]) |
|---|
| 705 | #self.set_wavelength(wavelength) |
|---|
| 706 | else: |
|---|
| 707 | raise |
|---|
| 708 | |
|---|
| 709 | def set_wave_spread(self, wavelength_spread): |
|---|
| 710 | """ |
|---|
| 711 | Set wavelength spread or wavelength spread |
|---|
| 712 | """ |
|---|
| 713 | if wavelength_spread.__class__.__name__ == 'list': |
|---|
| 714 | self.wave.set_wave_spread_list(wavelength_spread) |
|---|
| 715 | elif wavelength_spread.__class__.__name__ == 'float': |
|---|
| 716 | self.wave.set_wave_spread_list([wavelength_spread]) |
|---|
| 717 | else: |
|---|
| 718 | raise |
|---|
| 719 | |
|---|
| 720 | def set_wavelength(self, wavelength): |
|---|
| 721 | """ |
|---|
| 722 | Set wavelength |
|---|
| 723 | """ |
|---|
| 724 | self.wavelength = wavelength |
|---|
| 725 | self.wave.set_wavelength(wavelength) |
|---|
| 726 | |
|---|
| 727 | def set_spectrum(self, spectrum): |
|---|
| 728 | """ |
|---|
| 729 | Set spectrum |
|---|
| 730 | """ |
|---|
| 731 | self.spectrum = spectrum |
|---|
| 732 | self.wave.set_spectrum(spectrum) |
|---|
| 733 | |
|---|
| 734 | def set_wavelength_spread(self, wavelength_spread): |
|---|
| 735 | """ |
|---|
| 736 | Set wavelength spread |
|---|
| 737 | """ |
|---|
| 738 | self.wavelength_spread = wavelength_spread |
|---|
| 739 | self.wave.set_wavelength_spread(wavelength_spread) |
|---|
| 740 | |
|---|
| 741 | def set_wave_list(self, wavelength_list, wavelengthspread_list): |
|---|
| 742 | """ |
|---|
| 743 | Set wavelength and its spread list |
|---|
| 744 | """ |
|---|
| 745 | self.wave.set_wave_list(wavelength_list) |
|---|
| 746 | self.wave.set_wave_spread_list(wavelengthspread_list) |
|---|
| 747 | |
|---|
| 748 | def get_wave_list(self): |
|---|
| 749 | """ |
|---|
| 750 | Set wavelength spread |
|---|
| 751 | """ |
|---|
| 752 | return self.wave.get_wave_list() |
|---|
| 753 | |
|---|
| 754 | def get_intensity_list(self): |
|---|
| 755 | """ |
|---|
| 756 | Set wavelength spread |
|---|
| 757 | """ |
|---|
| 758 | return self.wave.get_intensity_list() |
|---|
| 759 | |
|---|
| 760 | def set_source_aperture_size(self, size): |
|---|
| 761 | """ |
|---|
| 762 | Set source aperture size |
|---|
| 763 | |
|---|
| 764 | : param size: [dia_value] or [x_value, y_value] |
|---|
| 765 | """ |
|---|
| 766 | if len(size) < 1 or len(size) > 2: |
|---|
| 767 | raise RuntimeError, "The length of the size must be one or two." |
|---|
| 768 | self.aperture.set_source_size(size) |
|---|
| 769 | |
|---|
| 770 | def set_neutron_mass(self, mass): |
|---|
| 771 | """ |
|---|
| 772 | Set Neutron mass |
|---|
| 773 | """ |
|---|
| 774 | self.wave.set_mass(mass) |
|---|
| 775 | self.mass = mass |
|---|
| 776 | |
|---|
| 777 | def set_sample_aperture_size(self, size): |
|---|
| 778 | """ |
|---|
| 779 | Set sample aperture size |
|---|
| 780 | |
|---|
| 781 | : param size: [dia_value] or [xheight_value, yheight_value] |
|---|
| 782 | """ |
|---|
| 783 | if len(size) < 1 or len(size) > 2: |
|---|
| 784 | raise RuntimeError, "The length of the size must be one or two." |
|---|
| 785 | self.aperture.set_sample_size(size) |
|---|
| 786 | |
|---|
| 787 | def set_detector_pix_size(self, size): |
|---|
| 788 | """ |
|---|
| 789 | Set detector pixel size |
|---|
| 790 | """ |
|---|
| 791 | self.detector.set_pix_size(size) |
|---|
| 792 | |
|---|
| 793 | def set_detector_size(self, size): |
|---|
| 794 | """ |
|---|
| 795 | Set detector size in number of pixels |
|---|
| 796 | : param size: [pixel_nums] or [x_pix_num, yx_pix_num] |
|---|
| 797 | """ |
|---|
| 798 | self.detector.set_size(size) |
|---|
| 799 | |
|---|
| 800 | def set_source2sample_distance(self, distance): |
|---|
| 801 | """ |
|---|
| 802 | Set detector source2sample_distance |
|---|
| 803 | |
|---|
| 804 | : param distance: [distance, x_offset] |
|---|
| 805 | """ |
|---|
| 806 | if len(distance) < 1 or len(distance) > 2: |
|---|
| 807 | raise RuntimeError, "The length of the size must be one or two." |
|---|
| 808 | self.aperture.set_sample_distance(distance) |
|---|
| 809 | |
|---|
| 810 | def set_sample2sample_distance(self, distance): |
|---|
| 811 | """ |
|---|
| 812 | Set detector sample_slit2sample_distance |
|---|
| 813 | |
|---|
| 814 | : param distance: [distance, x_offset] |
|---|
| 815 | """ |
|---|
| 816 | if len(distance) < 1 or len(distance) > 2: |
|---|
| 817 | raise RuntimeError, "The length of the size must be one or two." |
|---|
| 818 | self.sample.set_distance(distance) |
|---|
| 819 | |
|---|
| 820 | def set_sample2detector_distance(self, distance): |
|---|
| 821 | """ |
|---|
| 822 | Set detector sample2detector_distance |
|---|
| 823 | |
|---|
| 824 | : param distance: [distance, x_offset] |
|---|
| 825 | """ |
|---|
| 826 | if len(distance) < 1 or len(distance) > 2: |
|---|
| 827 | raise RuntimeError, "The length of the size must be one or two." |
|---|
| 828 | self.detector.set_distance(distance) |
|---|
| 829 | |
|---|
| 830 | def get_all_instrument_params(self): |
|---|
| 831 | """ |
|---|
| 832 | Get all instrumental parameters |
|---|
| 833 | """ |
|---|
| 834 | self.mass = self.get_neutron_mass() |
|---|
| 835 | self.spectrum = self.get_spectrum() |
|---|
| 836 | self.source_aperture_size = self.get_source_aperture_size() |
|---|
| 837 | self.sample_aperture_size = self.get_sample_aperture_size() |
|---|
| 838 | self.detector_pix_size = self.get_detector_pix_size() |
|---|
| 839 | self.detector_size = self.get_detector_size() |
|---|
| 840 | self.source2sample_distance = self.get_source2sample_distance() |
|---|
| 841 | self.sample2sample_distance = self.get_sample2sample_distance() |
|---|
| 842 | self.sample2detector_distance = self.get_sample2detector_distance() |
|---|
| 843 | |
|---|
| 844 | def get_detector_qrange(self): |
|---|
| 845 | """ |
|---|
| 846 | get max detector q ranges |
|---|
| 847 | |
|---|
| 848 | : return: qx_min, qx_max, qy_min, qy_max tuple |
|---|
| 849 | """ |
|---|
| 850 | if len(self.qxrange) != 2 or len(self.qyrange) != 2: |
|---|
| 851 | return None |
|---|
| 852 | qx_min = self.qxrange[0] |
|---|
| 853 | qx_max = self.qxrange[1] |
|---|
| 854 | qy_min = self.qyrange[0] |
|---|
| 855 | qy_max = self.qyrange[1] |
|---|
| 856 | |
|---|
| 857 | return qx_min, qx_max, qy_min, qy_max |
|---|
| 858 | |
|---|
| 859 | def _rotate_z(self, x_value, y_value, theta=0.0): |
|---|
| 860 | """ |
|---|
| 861 | Rotate x-y cordinate around z-axis by theta |
|---|
| 862 | : x_value: numpy array of x values |
|---|
| 863 | : y_value: numpy array of y values |
|---|
| 864 | : theta: angle to rotate by in rad |
|---|
| 865 | |
|---|
| 866 | :return: x_prime, y-prime |
|---|
| 867 | """ |
|---|
| 868 | # rotate by theta |
|---|
| 869 | x_prime = x_value * math.cos(theta) + y_value * math.sin(theta) |
|---|
| 870 | y_prime = -x_value * math.sin(theta) + y_value * math.cos(theta) |
|---|
| 871 | |
|---|
| 872 | return x_prime, y_prime |
|---|
| 873 | |
|---|
| 874 | def _gaussian2d(self, x_val, y_val, x0_val, y0_val, |
|---|
| 875 | sigma_x, sigma_y, sigma_r): |
|---|
| 876 | """ |
|---|
| 877 | Calculate 2D Gaussian distribution |
|---|
| 878 | : x_val: x value |
|---|
| 879 | : y_val: y value |
|---|
| 880 | : x0_val: mean value in x-axis |
|---|
| 881 | : y0_val: mean value in y-axis |
|---|
| 882 | : sigma_x: variance in x-direction |
|---|
| 883 | : sigma_y: variance in y-direction |
|---|
| 884 | |
|---|
| 885 | : return: gaussian (value) |
|---|
| 886 | """ |
|---|
| 887 | # phi values at each points (not at the center) |
|---|
| 888 | x_value = x_val - x0_val |
|---|
| 889 | y_value = y_val - y0_val |
|---|
| 890 | phi_i = np.arctan2(y_val, x_val) |
|---|
| 891 | |
|---|
| 892 | # phi correction due to the gravity shift (in phi) |
|---|
| 893 | phi_0 = math.atan2(y0_val, x0_val) |
|---|
| 894 | phi_i = phi_i - phi_0 + self.gravity_phi |
|---|
| 895 | |
|---|
| 896 | sin_phi = np.sin(self.gravity_phi) |
|---|
| 897 | cos_phi = np.cos(self.gravity_phi) |
|---|
| 898 | |
|---|
| 899 | x_p = x_value * cos_phi + y_value * sin_phi |
|---|
| 900 | y_p = -x_value * sin_phi + y_value * cos_phi |
|---|
| 901 | |
|---|
| 902 | new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x) + 1) |
|---|
| 903 | new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y) + 1) |
|---|
| 904 | new_x = x_p * cos_phi / new_sig_x - y_p * sin_phi |
|---|
| 905 | new_x /= sigma_x |
|---|
| 906 | new_y = x_p * sin_phi / new_sig_y + y_p * cos_phi |
|---|
| 907 | new_y /= sigma_y |
|---|
| 908 | |
|---|
| 909 | nu_value = -0.5 * (new_x * new_x + new_y * new_y) |
|---|
| 910 | |
|---|
| 911 | gaussian = np.exp(nu_value) |
|---|
| 912 | # normalizing factor correction |
|---|
| 913 | gaussian /= gaussian.sum() |
|---|
| 914 | |
|---|
| 915 | return gaussian |
|---|
| 916 | |
|---|
| 917 | def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, |
|---|
| 918 | sigma_x, sigma_y, sigma_r): |
|---|
| 919 | """ |
|---|
| 920 | Calculate 2D Gaussian distribution for polar coodinate |
|---|
| 921 | : x_val: x value |
|---|
| 922 | : y_val: y value |
|---|
| 923 | : x0_val: mean value in x-axis |
|---|
| 924 | : y0_val: mean value in y-axis |
|---|
| 925 | : sigma_x: variance in r-direction |
|---|
| 926 | : sigma_y: variance in phi-direction |
|---|
| 927 | : sigma_r: wavelength variance in r-direction |
|---|
| 928 | |
|---|
| 929 | : return: gaussian (value) |
|---|
| 930 | """ |
|---|
| 931 | sigma_x = sqrt(sigma_x * sigma_x + sigma_r * sigma_r) |
|---|
| 932 | # call gaussian1d |
|---|
| 933 | gaussian = self._gaussian1d(x_val, x0_val, sigma_x) |
|---|
| 934 | gaussian *= self._gaussian1d(y_val, y0_val, sigma_y) |
|---|
| 935 | |
|---|
| 936 | # normalizing factor correction |
|---|
| 937 | if sigma_x != 0 and sigma_y != 0: |
|---|
| 938 | gaussian *= sqrt(2 * pi) |
|---|
| 939 | return gaussian |
|---|
| 940 | |
|---|
| 941 | def _gaussian1d(self, value, mean, sigma): |
|---|
| 942 | """ |
|---|
| 943 | Calculate 1D Gaussian distribution |
|---|
| 944 | : value: value |
|---|
| 945 | : mean: mean value |
|---|
| 946 | : sigma: variance |
|---|
| 947 | |
|---|
| 948 | : return: gaussian (value) |
|---|
| 949 | """ |
|---|
| 950 | # default |
|---|
| 951 | gaussian = 1.0 |
|---|
| 952 | if sigma != 0: |
|---|
| 953 | # get exponent |
|---|
| 954 | nu_value = (value - mean) / sigma |
|---|
| 955 | nu_value *= nu_value |
|---|
| 956 | nu_value *= -0.5 |
|---|
| 957 | gaussian *= np.exp(nu_value) |
|---|
| 958 | gaussian /= sigma |
|---|
| 959 | # normalize |
|---|
| 960 | gaussian /= sqrt(2 * pi) |
|---|
| 961 | |
|---|
| 962 | return gaussian |
|---|
| 963 | |
|---|
| 964 | def _atan_phi(self, qy_value, qx_value): |
|---|
| 965 | """ |
|---|
| 966 | Find the angle phi of q on the detector plane for qx_value, qy_value given |
|---|
| 967 | : qx_value: x component of q |
|---|
| 968 | : qy_value: y component of q |
|---|
| 969 | |
|---|
| 970 | : return phi: the azimuthal angle of q on x-y plane |
|---|
| 971 | """ |
|---|
| 972 | phi = math.atan2(qy_value, qx_value) |
|---|
| 973 | return phi |
|---|
| 974 | |
|---|
| 975 | def _get_detector_qxqy_pixels(self): |
|---|
| 976 | """ |
|---|
| 977 | Get the pixel positions of the detector in the qx_value-qy_value space |
|---|
| 978 | """ |
|---|
| 979 | |
|---|
| 980 | # update all param values |
|---|
| 981 | self.get_all_instrument_params() |
|---|
| 982 | |
|---|
| 983 | # wavelength |
|---|
| 984 | wavelength = self.wave.wavelength |
|---|
| 985 | # Gavity correction |
|---|
| 986 | delta_y = self._get_beamcenter_drop() # in cm |
|---|
| 987 | |
|---|
| 988 | # detector_pix size |
|---|
| 989 | detector_pix_size = self.detector_pix_size |
|---|
| 990 | # Square or circular pixel |
|---|
| 991 | if len(detector_pix_size) == 1: |
|---|
| 992 | pix_x_size = detector_pix_size[0] |
|---|
| 993 | pix_y_size = detector_pix_size[0] |
|---|
| 994 | # rectangular pixel pixel |
|---|
| 995 | elif len(detector_pix_size) == 2: |
|---|
| 996 | pix_x_size = detector_pix_size[0] |
|---|
| 997 | pix_y_size = detector_pix_size[1] |
|---|
| 998 | else: |
|---|
| 999 | raise ValueError, " Input value format error..." |
|---|
| 1000 | # Sample to detector distance = sample slit to detector |
|---|
| 1001 | # minus sample offset |
|---|
| 1002 | sample2detector_distance = self.sample2detector_distance[0] - \ |
|---|
| 1003 | self.sample2sample_distance[0] |
|---|
| 1004 | # detector offset in x-direction |
|---|
| 1005 | detector_offset = 0 |
|---|
| 1006 | try: |
|---|
| 1007 | detector_offset = self.sample2detector_distance[1] |
|---|
| 1008 | except: |
|---|
| 1009 | logger.error(sys.exc_value) |
|---|
| 1010 | |
|---|
| 1011 | # detector size in [no of pix_x,no of pix_y] |
|---|
| 1012 | detector_pix_nums_x = self.detector_size[0] |
|---|
| 1013 | |
|---|
| 1014 | # get pix_y if it exists, otherwse take it from [0] |
|---|
| 1015 | try: |
|---|
| 1016 | detector_pix_nums_y = self.detector_size[1] |
|---|
| 1017 | except: |
|---|
| 1018 | detector_pix_nums_y = self.detector_size[0] |
|---|
| 1019 | |
|---|
| 1020 | # detector offset in pix number |
|---|
| 1021 | offset_x = detector_offset / pix_x_size |
|---|
| 1022 | offset_y = delta_y / pix_y_size |
|---|
| 1023 | |
|---|
| 1024 | # beam center position in pix number (start from 0) |
|---|
| 1025 | center_x, center_y = self._get_beamcenter_position(detector_pix_nums_x, |
|---|
| 1026 | detector_pix_nums_y, |
|---|
| 1027 | offset_x, offset_y) |
|---|
| 1028 | # distance [cm] from the beam center on detector plane |
|---|
| 1029 | detector_ind_x = np.arange(detector_pix_nums_x) |
|---|
| 1030 | detector_ind_y = np.arange(detector_pix_nums_y) |
|---|
| 1031 | |
|---|
| 1032 | # shif 0.5 pixel so that pix position is at the center of the pixel |
|---|
| 1033 | detector_ind_x = detector_ind_x + 0.5 |
|---|
| 1034 | detector_ind_y = detector_ind_y + 0.5 |
|---|
| 1035 | |
|---|
| 1036 | # the relative postion from the beam center |
|---|
| 1037 | detector_ind_x = detector_ind_x - center_x |
|---|
| 1038 | detector_ind_y = detector_ind_y - center_y |
|---|
| 1039 | |
|---|
| 1040 | # unit correction in cm |
|---|
| 1041 | detector_ind_x = detector_ind_x * pix_x_size |
|---|
| 1042 | detector_ind_y = detector_ind_y * pix_y_size |
|---|
| 1043 | |
|---|
| 1044 | qx_value = np.zeros(len(detector_ind_x)) |
|---|
| 1045 | qy_value = np.zeros(len(detector_ind_y)) |
|---|
| 1046 | i = 0 |
|---|
| 1047 | |
|---|
| 1048 | for indx in detector_ind_x: |
|---|
| 1049 | qx_value[i] = self._get_qx(indx, sample2detector_distance, wavelength) |
|---|
| 1050 | i += 1 |
|---|
| 1051 | i = 0 |
|---|
| 1052 | for indy in detector_ind_y: |
|---|
| 1053 | qy_value[i] = self._get_qx(indy, sample2detector_distance, wavelength) |
|---|
| 1054 | i += 1 |
|---|
| 1055 | |
|---|
| 1056 | # qx_value and qy_value values in array |
|---|
| 1057 | qx_value = qx_value.repeat(detector_pix_nums_y) |
|---|
| 1058 | qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y) |
|---|
| 1059 | qy_value = qy_value.repeat(detector_pix_nums_x) |
|---|
| 1060 | qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x) |
|---|
| 1061 | qy_value = qy_value.transpose() |
|---|
| 1062 | |
|---|
| 1063 | # p min and max values among the center of pixels |
|---|
| 1064 | self.qx_min = np.min(qx_value) |
|---|
| 1065 | self.qx_max = np.max(qx_value) |
|---|
| 1066 | self.qy_min = np.min(qy_value) |
|---|
| 1067 | self.qy_max = np.max(qy_value) |
|---|
| 1068 | |
|---|
| 1069 | # Appr. min and max values of the detector display limits |
|---|
| 1070 | # i.e., edges of the last pixels. |
|---|
| 1071 | self.qy_min += self._get_qx(-0.5 * pix_y_size, |
|---|
| 1072 | sample2detector_distance, wavelength) |
|---|
| 1073 | self.qy_max += self._get_qx(0.5 * pix_y_size, |
|---|
| 1074 | sample2detector_distance, wavelength) |
|---|
| 1075 | #if self.qx_min == self.qx_max: |
|---|
| 1076 | self.qx_min += self._get_qx(-0.5 * pix_x_size, |
|---|
| 1077 | sample2detector_distance, wavelength) |
|---|
| 1078 | self.qx_max += self._get_qx(0.5 * pix_x_size, |
|---|
| 1079 | sample2detector_distance, wavelength) |
|---|
| 1080 | |
|---|
| 1081 | # min and max values of detecter |
|---|
| 1082 | self.detector_qx_min = self.qx_min |
|---|
| 1083 | self.detector_qx_max = self.qx_max |
|---|
| 1084 | self.detector_qy_min = self.qy_min |
|---|
| 1085 | self.detector_qy_max = self.qy_max |
|---|
| 1086 | |
|---|
| 1087 | # try to set it as a Data2D otherwise pass (not required for now) |
|---|
| 1088 | try: |
|---|
| 1089 | from sas.sascalc.dataloader.data_info import Data2D |
|---|
| 1090 | output = Data2D() |
|---|
| 1091 | inten = np.zeros_like(qx_value) |
|---|
| 1092 | output.data = inten |
|---|
| 1093 | output.qx_data = qx_value |
|---|
| 1094 | output.qy_data = qy_value |
|---|
| 1095 | except: |
|---|
| 1096 | logger.error(sys.exc_value) |
|---|
| 1097 | |
|---|
| 1098 | return output |
|---|
| 1099 | |
|---|
| 1100 | def _get_qx(self, dx_size, det_dist, wavelength): |
|---|
| 1101 | """ |
|---|
| 1102 | :param dx_size: x-distance from beam center [cm] |
|---|
| 1103 | :param det_dist: sample to detector distance [cm] |
|---|
| 1104 | |
|---|
| 1105 | :return: q-value at the given position |
|---|
| 1106 | """ |
|---|
| 1107 | # Distance from beam center in the plane of detector |
|---|
| 1108 | plane_dist = dx_size |
|---|
| 1109 | # full scattering angle on the x-axis |
|---|
| 1110 | theta = np.arctan(plane_dist / det_dist) |
|---|
| 1111 | qx_value = (2.0 * pi / wavelength) * np.sin(theta) |
|---|
| 1112 | return qx_value |
|---|
| 1113 | |
|---|
| 1114 | def _get_polar_value(self, qx_value, qy_value): |
|---|
| 1115 | """ |
|---|
| 1116 | Find qr_value and phi from qx_value and qy_value values |
|---|
| 1117 | |
|---|
| 1118 | : return qr_value, phi |
|---|
| 1119 | """ |
|---|
| 1120 | # find |q| on detector plane |
|---|
| 1121 | qr_value = sqrt(qx_value*qx_value + qy_value*qy_value) |
|---|
| 1122 | # find angle phi |
|---|
| 1123 | phi = self._atan_phi(qy_value, qx_value) |
|---|
| 1124 | |
|---|
| 1125 | return qr_value, phi |
|---|
| 1126 | |
|---|
| 1127 | def _get_beamcenter_position(self, num_x, num_y, offset_x, offset_y): |
|---|
| 1128 | """ |
|---|
| 1129 | :param num_x: number of pixel in x-direction |
|---|
| 1130 | :param num_y: number of pixel in y-direction |
|---|
| 1131 | :param offset: detector offset in x-direction in pix number |
|---|
| 1132 | |
|---|
| 1133 | :return: pix number; pos_x, pos_y in pix index |
|---|
| 1134 | """ |
|---|
| 1135 | # beam center position |
|---|
| 1136 | pos_x = num_x / 2 |
|---|
| 1137 | pos_y = num_y / 2 |
|---|
| 1138 | |
|---|
| 1139 | # correction for offset |
|---|
| 1140 | pos_x += offset_x |
|---|
| 1141 | # correction for gravity that is always negative |
|---|
| 1142 | pos_y -= offset_y |
|---|
| 1143 | |
|---|
| 1144 | return pos_x, pos_y |
|---|
| 1145 | |
|---|
| 1146 | def _get_beamcenter_drop(self): |
|---|
| 1147 | """ |
|---|
| 1148 | Get the beam center drop (delta y) in y diection due to gravity |
|---|
| 1149 | |
|---|
| 1150 | :return delta y: the beam center drop in cm |
|---|
| 1151 | """ |
|---|
| 1152 | # Check if mass == 0 (X-ray). |
|---|
| 1153 | if self.mass == 0: |
|---|
| 1154 | return 0 |
|---|
| 1155 | # Covert unit from A to cm |
|---|
| 1156 | unit_cm = 1e-08 |
|---|
| 1157 | # Velocity of neutron in horizontal direction (~ actual velocity) |
|---|
| 1158 | velocity = _PLANK_H / (self.mass * self.wave.wavelength * unit_cm) |
|---|
| 1159 | # Compute delta y |
|---|
| 1160 | delta_y = 0.5 |
|---|
| 1161 | delta_y *= _GRAVITY |
|---|
| 1162 | sampletodetector = self.sample2detector_distance[0] - \ |
|---|
| 1163 | self.sample2sample_distance[0] |
|---|
| 1164 | delta_y *= sampletodetector |
|---|
| 1165 | delta_y *= (self.source2sample_distance[0] + self.sample2detector_distance[0]) |
|---|
| 1166 | delta_y /= (velocity * velocity) |
|---|
| 1167 | |
|---|
| 1168 | return delta_y |
|---|
