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