""" This object is a small tool to allow user to quickly determine the variance in q from the instrumental parameters. """ from instrument import Sample from instrument import Detector from instrument import TOF as Neutron from instrument import Aperture # import math stuffs from math import pi from math import sqrt import math import numpy as np import sys import logging #Plank's constant in cgs unit _PLANK_H = 6.62606896E-27 #Gravitational acc. in cgs unit _GRAVITY = 981.0 class ResolutionCalculator(object): """ compute resolution in 2D """ def __init__(self): # wavelength self.wave = Neutron() # sample self.sample = Sample() # aperture self.aperture = Aperture() # detector self.detector = Detector() # 2d image of the resolution self.image = [] self.image_lam = [] # resolutions # lamda in r-direction self.sigma_lamd = 0 # x-dir (no lamda) self.sigma_1 = 0 #y-dir (no lamda) self.sigma_2 = 0 # 1D total self.sigma_1d = 0 self.gravity_phi = None # q min and max self.qx_min = -0.3 self.qx_max = 0.3 self.qy_min = -0.3 self.qy_max = 0.3 # q min and max of the detector self.detector_qx_min = -0.3 self.detector_qx_max = 0.3 self.detector_qy_min = -0.3 self.detector_qy_max = 0.3 # possible max qrange self.qxmin_limit = 0 self.qxmax_limit = 0 self.qymin_limit = 0 self.qymax_limit = 0 # plots self.plot = None # instrumental params defaults self.mass = 0 self.intensity = 0 self.wavelength = 0 self.wavelength_spread = 0 self.source_aperture_size = [] self.source2sample_distance = [] self.sample2sample_distance = [] self.sample_aperture_size = [] self.sample2detector_distance = [] self.detector_pix_size = [] self.detector_size = [] self.get_all_instrument_params() # max q range for all lambdas self.qxrange = [] self.qyrange = [] def compute_and_plot(self, qx_value, qy_value, qx_min, qx_max, qy_min, qy_max, coord='cartesian'): """ Compute the resolution : qx_value: x component of q : qy_value: y component of q """ # make sure to update all the variables need. # except lambda, dlambda, and intensity self.get_all_instrument_params() # wavelength etc. lamda_list, dlamb_list = self.get_wave_list() intens_list = [] sig1_list = [] sig2_list = [] sigr_list = [] sigma1d_list = [] num_lamda = len(lamda_list) for num in range(num_lamda): lam = lamda_list[num] # wavelength spread dlam = dlamb_list[num] intens = self.setup_tof(lam, dlam) intens_list.append(intens) # cehck if tof if num_lamda > 1: tof = True else: tof = False # compute 2d resolution _, _, sigma_1, sigma_2, sigma_r, sigma1d = \ self.compute(lam, dlam, qx_value, qy_value, coord, tof) # make image image = self.get_image(qx_value, qy_value, sigma_1, sigma_2, sigma_r, qx_min, qx_max, qy_min, qy_max, coord, False) if qx_min > self.qx_min: qx_min = self.qx_min if qx_max < self.qx_max: qx_max = self.qx_max if qy_min > self.qy_min: qy_min = self.qy_min if qy_max < self.qy_max: qy_max = self.qy_max # set max qranges self.qxrange = [qx_min, qx_max] self.qyrange = [qy_min, qy_max] sig1_list.append(sigma_1) sig2_list.append(sigma_2) sigr_list.append(sigma_r) sigma1d_list.append(sigma1d) # redraw image in global 2d q-space. self.image_lam = [] total_intensity = 0 sigma_1 = 0 sigma_r = 0 sigma_2 = 0 sigma1d = 0 for ind in range(num_lamda): lam = lamda_list[ind] dlam = dlamb_list[ind] intens = self.setup_tof(lam, dlam) out = self.get_image(qx_value, qy_value, sig1_list[ind], sig2_list[ind], sigr_list[ind], qx_min, qx_max, qy_min, qy_max, coord) # this is the case of q being outside the detector #if numpy.all(out==0.0): # continue image = out # set variance as sigmas sigma_1 += sig1_list[ind] * sig1_list[ind] * self.intensity sigma_r += sigr_list[ind] * sigr_list[ind] * self.intensity sigma_2 += sig2_list[ind] * sig2_list[ind] * self.intensity sigma1d += sigma1d_list[ind] * sigma1d_list[ind] * self.intensity total_intensity += self.intensity if total_intensity != 0: # average variance image_out = image / total_intensity sigma_1 = sigma_1 / total_intensity sigma_r = sigma_r / total_intensity sigma_2 = sigma_2 / total_intensity sigma1d = sigma1d / total_intensity # set sigmas self.sigma_1 = sqrt(sigma_1) self.sigma_lamd = sqrt(sigma_r) self.sigma_2 = sqrt(sigma_2) self.sigma_1d = sqrt(sigma1d) # rescale max_im_val = 1 if max_im_val > 0: image_out /= max_im_val else: image_out = image * 0.0 # Don't calculate sigmas nor set self.sigmas! sigma_1 = 0 sigma_r = 0 sigma_2 = 0 sigma1d = 0 if len(self.image) > 0: self.image += image_out else: self.image = image_out # plot image return self.plot_image(self.image) def setup_tof(self, wavelength, wavelength_spread): """ Setup all parameters in instrument : param ind: index of lambda, etc """ # set wave.wavelength self.set_wavelength(wavelength) self.set_wavelength_spread(wavelength_spread) self.intensity = self.wave.get_intensity() if wavelength == 0: msg = "Can't compute the resolution: the wavelength is zero..." raise RuntimeError, msg return self.intensity def compute(self, wavelength, wavelength_spread, qx_value, qy_value, coord='cartesian', tof=False): """ Compute the Q resoltuion in || and + direction of 2D : qx_value: x component of q : qy_value: y component of q """ coord = 'cartesian' lamb = wavelength lamb_spread = wavelength_spread # the shape of wavelength distribution if tof: # rectangular tof_factor = 2 else: # triangular tof_factor = 1 # Find polar values qr_value, phi = self._get_polar_value(qx_value, qy_value) # vacuum wave transfer knot = 2*pi/lamb # scattering angle theta; always true for plane detector # aligned vertically to the ko direction if qr_value > knot: theta = pi/2 else: theta = math.asin(qr_value/knot) # source aperture size rone = self.source_aperture_size # sample aperture size rtwo = self.sample_aperture_size # detector pixel size rthree = self.detector_pix_size # source to sample(aperture) distance l_ssa = self.source2sample_distance[0] # sample(aperture) to detector distance l_sad = self.sample2detector_distance[0] # sample (aperture) to sample distance l_sas = self.sample2sample_distance[0] # source to sample distance l_one = l_ssa + l_sas # sample to detector distance l_two = l_sad - l_sas # Sample offset correction for l_one and Lp on variance calculation l1_cor = (l_ssa * l_two) / (l_sas + l_two) lp_cor = (l_ssa * l_two) / (l_one + l_two) # the radial distance to the pixel from the center of the detector radius = math.tan(theta) * l_two #Lp = l_one*l_two/(l_one+l_two) # default polar coordinate comp1 = 'radial' comp2 = 'phi' # in the case of the cartesian coordinate if coord == 'cartesian': comp1 = 'x' comp2 = 'y' # sigma in the radial/x direction # for source aperture sigma_1 = self.get_variance(rone, l1_cor, phi, comp1) # for sample apperture sigma_1 += self.get_variance(rtwo, lp_cor, phi, comp1) # for detector pix sigma_1 += self.get_variance(rthree, l_two, phi, comp1) # for gravity term for 1d sigma_1grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, phi, comp1, 'on') / tof_factor # for wavelength spread # reserve for 1d calculation A_value = self._cal_A_value(lamb, l_ssa, l_sad) sigma_wave_1, sigma_wave_1_1d = self.get_variance_wave(A_value, radius, l_two, lamb_spread, phi, 'radial', 'on') sigma_wave_1 /= tof_factor sigma_wave_1_1d /= tof_factor # for 1d variance_1d_1 = (sigma_1 + sigma_1grav1d) / 2 + sigma_wave_1_1d # normalize variance_1d_1 = knot * knot * variance_1d_1 / 12 # for 2d #sigma_1 += sigma_wave_1 # normalize sigma_1 = knot * sqrt(sigma_1 / 12) sigma_r = knot * sqrt(sigma_wave_1 / (tof_factor *12)) # sigma in the phi/y direction # for source apperture sigma_2 = self.get_variance(rone, l1_cor, phi, comp2) # for sample apperture sigma_2 += self.get_variance(rtwo, lp_cor, phi, comp2) # for detector pix sigma_2 += self.get_variance(rthree, l_two, phi, comp2) # for gravity term for 1d sigma_2grav1d = self.get_variance_gravity(l_ssa, l_sad, lamb, lamb_spread, phi, comp2, 'on') / tof_factor # for wavelength spread # reserve for 1d calculation sigma_wave_2, sigma_wave_2_1d = self.get_variance_wave(A_value, radius, l_two, lamb_spread, phi, 'phi', 'on') sigma_wave_2 /= tof_factor sigma_wave_2_1d /= tof_factor # for 1d variance_1d_2 = (sigma_2 + sigma_2grav1d) / 2 + sigma_wave_2_1d # normalize variance_1d_2 = knot * knot * variance_1d_2 / 12 # for 2d #sigma_2 = knot*sqrt(sigma_2/12) #sigma_2 += sigma_wave_2 # normalize sigma_2 = knot * sqrt(sigma_2 / 12) sigma1d = sqrt(variance_1d_1 + variance_1d_2) # set sigmas self.sigma_1 = sigma_1 self.sigma_lamd = sigma_r self.sigma_2 = sigma_2 self.sigma_1d = sigma1d return qr_value, phi, sigma_1, sigma_2, sigma_r, sigma1d def _within_detector_range(self, qx_value, qy_value): """ check if qvalues are within detector range """ # detector range detector_qx_min = self.detector_qx_min detector_qx_max = self.detector_qx_max detector_qy_min = self.detector_qy_min detector_qy_max = self.detector_qy_max if self.qxmin_limit > detector_qx_min: self.qxmin_limit = detector_qx_min if self.qxmax_limit < detector_qx_max: self.qxmax_limit = detector_qx_max if self.qymin_limit > detector_qy_min: self.qymin_limit = detector_qy_min if self.qymax_limit < detector_qy_max: self.qymax_limit = detector_qy_max if qx_value < detector_qx_min or qx_value > detector_qx_max: return False if qy_value < detector_qy_min or qy_value > detector_qy_max: return False return True def get_image(self, qx_value, qy_value, sigma_1, sigma_2, sigma_r, qx_min, qx_max, qy_min, qy_max, coord='cartesian', full_cal=True): """ Get the resolution in polar coordinate ready to plot : qx_value: qx_value value : qy_value: qy_value value : sigma_1: variance in r direction : sigma_2: variance in phi direction : coord: coordinate system of image, 'polar' or 'cartesian' """ # Get qx_max and qy_max... self._get_detector_qxqy_pixels() qr_value, phi = self._get_polar_value(qx_value, qy_value) # Check whether the q value is within the detector range if qx_min < self.qx_min: self.qx_min = qx_min #raise ValueError, msg if qx_max > self.qx_max: self.qx_max = qx_max #raise ValueError, msg if qy_min < self.qy_min: self.qy_min = qy_min #raise ValueError, msg if qy_max > self.qy_max: self.qy_max = qy_max #raise ValueError, msg if not full_cal: return None # Make an empty graph in the detector scale dx_size = (self.qx_max - self.qx_min) / (1000 - 1) dy_size = (self.qy_max - self.qy_min) / (1000 - 1) x_val = np.arange(self.qx_min, self.qx_max, dx_size) y_val = np.arange(self.qy_max, self.qy_min, -dy_size) q_1, q_2 = np.meshgrid(x_val, y_val) #q_phi = numpy.arctan(q_1,q_2) # check whether polar or cartesian if coord == 'polar': # Find polar values qr_value, phi = self._get_polar_value(qx_value, qy_value) q_1, q_2 = self._rotate_z(q_1, q_2, phi) qc_1 = qr_value qc_2 = 0.0 # Calculate the 2D Gaussian distribution image image = self._gaussian2d_polar(q_1, q_2, qc_1, qc_2, sigma_1, sigma_2, sigma_r) else: # catesian coordinate # qx_center qc_1 = qx_value # qy_center qc_2 = qy_value # Calculate the 2D Gaussian distribution image image = self._gaussian2d(q_1, q_2, qc_1, qc_2, sigma_1, sigma_2, sigma_r) # out side of detector if not self._within_detector_range(qx_value, qy_value): image *= 0.0 self.intensity = 0.0 #return self.image # Add it if there are more than one inputs. if len(self.image_lam) > 0: self.image_lam += image * self.intensity else: self.image_lam = image * self.intensity return self.image_lam def plot_image(self, image): """ Plot image using pyplot : image: 2d resolution image : return plt: pylab object """ import matplotlib.pyplot as plt self.plot = plt plt.xlabel('$\\rm{Q}_{x} [A^{-1}]$') plt.ylabel('$\\rm{Q}_{y} [A^{-1}]$') # Max value of the image # max = numpy.max(image) qx_min, qx_max, qy_min, qy_max = self.get_detector_qrange() # Image im = plt.imshow(image, extent=[qx_min, qx_max, qy_min, qy_max]) # bilinear interpolation to make it smoother im.set_interpolation('bilinear') return plt def reset_image(self): """ Reset image to default (=[]) """ self.image = [] def get_variance(self, size=[], distance=0, phi=0, comp='radial'): """ Get the variance when the slit/pinhole size is given : size: list that can be one(diameter for circular) or two components(lengths for rectangular) : distance: [z, x] where z along the incident beam, x // qx_value : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial' : return variance: sigma^2 """ # check the length of size (list) len_size = len(size) # define sigma component direction if comp == 'radial': phi_x = math.cos(phi) phi_y = math.sin(phi) elif comp == 'phi': phi_x = math.sin(phi) phi_y = math.cos(phi) elif comp == 'x': phi_x = 1 phi_y = 0 elif comp == 'y': phi_x = 0 phi_y = 1 else: phi_x = 0 phi_y = 0 # calculate each component # for pinhole w/ radius = size[0]/2 if len_size == 1: x_comp = (0.5 * size[0]) * sqrt(3) y_comp = 0 # for rectangular slit elif len_size == 2: x_comp = size[0] * phi_x y_comp = size[1] * phi_y # otherwise else: raise ValueError, " Improper input..." # get them squared sigma = x_comp * x_comp sigma += y_comp * y_comp # normalize by distance sigma /= (distance * distance) return sigma def get_variance_wave(self, A_value, radius, distance, spread, phi, comp='radial', switch='on'): """ Get the variance when the wavelength spread is given : radius: the radial distance from the beam center to the pix of q : distance: sample to detector distance : spread: wavelength spread (ratio) : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial' : return variance: sigma^2 for 2d, sigma^2 for 1d [tuple] """ if switch.lower() == 'off': return 0, 0 # check the singular point if distance == 0 or comp == 'phi': return 0, 0 else: # calculate sigma^2 for 1d sigma1d = 2 * math.pow(radius/distance*spread, 2) if comp == 'x': sigma1d *= (math.cos(phi)*math.cos(phi)) elif comp == 'y': sigma1d *= (math.sin(phi)*math.sin(phi)) else: sigma1d *= 1 # sigma^2 for 2d # shift the coordinate due to the gravitational shift rad_x = radius * math.cos(phi) rad_y = A_value - radius * math.sin(phi) radius = math.sqrt(rad_x * rad_x + rad_y * rad_y) # new phi phi = math.atan2(-rad_y, rad_x) self.gravity_phi = phi # calculate sigma^2 sigma = 2 * math.pow(radius/distance*spread, 2) if comp == 'x': sigma *= (math.cos(phi)*math.cos(phi)) elif comp == 'y': sigma *= (math.sin(phi)*math.sin(phi)) else: sigma *= 1 return sigma, sigma1d def get_variance_gravity(self, s_distance, d_distance, wavelength, spread, phi, comp='radial', switch='on'): """ Get the variance from gravity when the wavelength spread is given : s_distance: source to sample distance : d_distance: sample to detector distance : wavelength: wavelength : spread: wavelength spread (ratio) : comp: direction of the sigma; can be 'phi', 'y', 'x', and 'radial' : return variance: sigma^2 """ if switch.lower() == 'off': return 0 if self.mass == 0.0: return 0 # check the singular point if d_distance == 0 or comp == 'x': return 0 else: a_value = self._cal_A_value(None, s_distance, d_distance) # calculate sigma^2 sigma = math.pow(a_value / d_distance, 2) sigma *= math.pow(wavelength, 4) sigma *= math.pow(spread, 2) sigma *= 8 return sigma def _cal_A_value(self, lamda, s_distance, d_distance): """ Calculate A value for gravity : s_distance: source to sample distance : d_distance: sample to detector distance """ # neutron mass in cgs unit self.mass = self.get_neutron_mass() # plank constant in cgs unit h_constant = _PLANK_H # gravity in cgs unit gravy = _GRAVITY # m/h m_over_h = self.mass / h_constant # A value a_value = d_distance * (s_distance + d_distance) a_value *= math.pow(m_over_h / 2, 2) a_value *= gravy # unit correction (1/cm to 1/A) for A and d_distance below a_value *= 1.0E-16 # if lamda is give (broad meanning of A) return 2* lamda^2 * A if lamda != None: a_value *= (4 * lamda * lamda) return a_value def get_intensity(self): """ Get intensity """ return self.wave.intensity def get_wavelength(self): """ Get wavelength """ return self.wave.wavelength def get_default_spectrum(self): """ Get default_spectrum """ return self.wave.get_default_spectrum() def get_spectrum(self): """ Get _spectrum """ return self.wave.get_spectrum() def get_wavelength_spread(self): """ Get wavelength spread """ return self.wave.wavelength_spread def get_neutron_mass(self): """ Get Neutron mass """ return self.wave.mass def get_source_aperture_size(self): """ Get source aperture size """ return self.aperture.source_size def get_sample_aperture_size(self): """ Get sample aperture size """ return self.aperture.sample_size def get_detector_pix_size(self): """ Get detector pixel size """ return self.detector.pix_size def get_detector_size(self): """ Get detector size """ return self.detector.size def get_source2sample_distance(self): """ Get detector source2sample_distance """ return self.aperture.sample_distance def get_sample2sample_distance(self): """ Get detector sampleslitsample_distance """ return self.sample.distance def get_sample2detector_distance(self): """ Get detector sample2detector_distance """ return self.detector.distance def set_intensity(self, intensity): """ Set intensity """ self.wave.set_intensity(intensity) def set_wave(self, wavelength): """ Set wavelength list or wavelength """ if wavelength.__class__.__name__ == 'list': self.wave.set_wave_list(wavelength) elif wavelength.__class__.__name__ == 'float': self.wave.set_wave_list([wavelength]) #self.set_wavelength(wavelength) else: raise def set_wave_spread(self, wavelength_spread): """ Set wavelength spread or wavelength spread """ if wavelength_spread.__class__.__name__ == 'list': self.wave.set_wave_spread_list(wavelength_spread) elif wavelength_spread.__class__.__name__ == 'float': self.wave.set_wave_spread_list([wavelength_spread]) else: raise def set_wavelength(self, wavelength): """ Set wavelength """ self.wavelength = wavelength self.wave.set_wavelength(wavelength) def set_spectrum(self, spectrum): """ Set spectrum """ self.spectrum = spectrum self.wave.set_spectrum(spectrum) def set_wavelength_spread(self, wavelength_spread): """ Set wavelength spread """ self.wavelength_spread = wavelength_spread self.wave.set_wavelength_spread(wavelength_spread) def set_wave_list(self, wavelength_list, wavelengthspread_list): """ Set wavelength and its spread list """ self.wave.set_wave_list(wavelength_list) self.wave.set_wave_spread_list(wavelengthspread_list) def get_wave_list(self): """ Set wavelength spread """ return self.wave.get_wave_list() def get_intensity_list(self): """ Set wavelength spread """ return self.wave.get_intensity_list() def set_source_aperture_size(self, size): """ Set source aperture size : param size: [dia_value] or [x_value, y_value] """ if len(size) < 1 or len(size) > 2: raise RuntimeError, "The length of the size must be one or two." self.aperture.set_source_size(size) def set_neutron_mass(self, mass): """ Set Neutron mass """ self.wave.set_mass(mass) self.mass = mass def set_sample_aperture_size(self, size): """ Set sample aperture size : param size: [dia_value] or [xheight_value, yheight_value] """ if len(size) < 1 or len(size) > 2: raise RuntimeError, "The length of the size must be one or two." self.aperture.set_sample_size(size) def set_detector_pix_size(self, size): """ Set detector pixel size """ self.detector.set_pix_size(size) def set_detector_size(self, size): """ Set detector size in number of pixels : param size: [pixel_nums] or [x_pix_num, yx_pix_num] """ self.detector.set_size(size) def set_source2sample_distance(self, distance): """ Set detector source2sample_distance : param distance: [distance, x_offset] """ if len(distance) < 1 or len(distance) > 2: raise RuntimeError, "The length of the size must be one or two." self.aperture.set_sample_distance(distance) def set_sample2sample_distance(self, distance): """ Set detector sample_slit2sample_distance : param distance: [distance, x_offset] """ if len(distance) < 1 or len(distance) > 2: raise RuntimeError, "The length of the size must be one or two." self.sample.set_distance(distance) def set_sample2detector_distance(self, distance): """ Set detector sample2detector_distance : param distance: [distance, x_offset] """ if len(distance) < 1 or len(distance) > 2: raise RuntimeError, "The length of the size must be one or two." self.detector.set_distance(distance) def get_all_instrument_params(self): """ Get all instrumental parameters """ self.mass = self.get_neutron_mass() self.spectrum = self.get_spectrum() self.source_aperture_size = self.get_source_aperture_size() self.sample_aperture_size = self.get_sample_aperture_size() self.detector_pix_size = self.get_detector_pix_size() self.detector_size = self.get_detector_size() self.source2sample_distance = self.get_source2sample_distance() self.sample2sample_distance = self.get_sample2sample_distance() self.sample2detector_distance = self.get_sample2detector_distance() def get_detector_qrange(self): """ get max detector q ranges : return: qx_min, qx_max, qy_min, qy_max tuple """ if len(self.qxrange) != 2 or len(self.qyrange) != 2: return None qx_min = self.qxrange[0] qx_max = self.qxrange[1] qy_min = self.qyrange[0] qy_max = self.qyrange[1] return qx_min, qx_max, qy_min, qy_max def _rotate_z(self, x_value, y_value, theta=0.0): """ Rotate x-y cordinate around z-axis by theta : x_value: numpy array of x values : y_value: numpy array of y values : theta: angle to rotate by in rad :return: x_prime, y-prime """ # rotate by theta x_prime = x_value * math.cos(theta) + y_value * math.sin(theta) y_prime = -x_value * math.sin(theta) + y_value * math.cos(theta) return x_prime, y_prime def _gaussian2d(self, x_val, y_val, x0_val, y0_val, sigma_x, sigma_y, sigma_r): """ Calculate 2D Gaussian distribution : x_val: x value : y_val: y value : x0_val: mean value in x-axis : y0_val: mean value in y-axis : sigma_x: variance in x-direction : sigma_y: variance in y-direction : return: gaussian (value) """ # phi values at each points (not at the center) x_value = x_val - x0_val y_value = y_val - y0_val phi_i = np.arctan2(y_val, x_val) # phi correction due to the gravity shift (in phi) phi_0 = math.atan2(y0_val, x0_val) phi_i = phi_i - phi_0 + self.gravity_phi sin_phi = np.sin(self.gravity_phi) cos_phi = np.cos(self.gravity_phi) x_p = x_value * cos_phi + y_value * sin_phi y_p = -x_value * sin_phi + y_value * cos_phi new_sig_x = sqrt(sigma_r * sigma_r / (sigma_x * sigma_x) + 1) new_sig_y = sqrt(sigma_r * sigma_r / (sigma_y * sigma_y) + 1) new_x = x_p * cos_phi / new_sig_x - y_p * sin_phi new_x /= sigma_x new_y = x_p * sin_phi / new_sig_y + y_p * cos_phi new_y /= sigma_y nu_value = -0.5 * (new_x * new_x + new_y * new_y) gaussian = np.exp(nu_value) # normalizing factor correction gaussian /= gaussian.sum() return gaussian def _gaussian2d_polar(self, x_val, y_val, x0_val, y0_val, sigma_x, sigma_y, sigma_r): """ Calculate 2D Gaussian distribution for polar coodinate : x_val: x value : y_val: y value : x0_val: mean value in x-axis : y0_val: mean value in y-axis : sigma_x: variance in r-direction : sigma_y: variance in phi-direction : sigma_r: wavelength variance in r-direction : return: gaussian (value) """ sigma_x = sqrt(sigma_x * sigma_x + sigma_r * sigma_r) # call gaussian1d gaussian = self._gaussian1d(x_val, x0_val, sigma_x) gaussian *= self._gaussian1d(y_val, y0_val, sigma_y) # normalizing factor correction if sigma_x != 0 and sigma_y != 0: gaussian *= sqrt(2 * pi) return gaussian def _gaussian1d(self, value, mean, sigma): """ Calculate 1D Gaussian distribution : value: value : mean: mean value : sigma: variance : return: gaussian (value) """ # default gaussian = 1.0 if sigma != 0: # get exponent nu_value = (value - mean) / sigma nu_value *= nu_value nu_value *= -0.5 gaussian *= np.exp(nu_value) gaussian /= sigma # normalize gaussian /= sqrt(2 * pi) return gaussian def _atan_phi(self, qy_value, qx_value): """ Find the angle phi of q on the detector plane for qx_value, qy_value given : qx_value: x component of q : qy_value: y component of q : return phi: the azimuthal angle of q on x-y plane """ phi = math.atan2(qy_value, qx_value) return phi def _get_detector_qxqy_pixels(self): """ Get the pixel positions of the detector in the qx_value-qy_value space """ # update all param values self.get_all_instrument_params() # wavelength wavelength = self.wave.wavelength # Gavity correction delta_y = self._get_beamcenter_drop() # in cm # detector_pix size detector_pix_size = self.detector_pix_size # Square or circular pixel if len(detector_pix_size) == 1: pix_x_size = detector_pix_size[0] pix_y_size = detector_pix_size[0] # rectangular pixel pixel elif len(detector_pix_size) == 2: pix_x_size = detector_pix_size[0] pix_y_size = detector_pix_size[1] else: raise ValueError, " Input value format error..." # Sample to detector distance = sample slit to detector # minus sample offset sample2detector_distance = self.sample2detector_distance[0] - \ self.sample2sample_distance[0] # detector offset in x-direction detector_offset = 0 try: detector_offset = self.sample2detector_distance[1] except: logging.error(sys.exc_value) # detector size in [no of pix_x,no of pix_y] detector_pix_nums_x = self.detector_size[0] # get pix_y if it exists, otherwse take it from [0] try: detector_pix_nums_y = self.detector_size[1] except: detector_pix_nums_y = self.detector_size[0] # detector offset in pix number offset_x = detector_offset / pix_x_size offset_y = delta_y / pix_y_size # beam center position in pix number (start from 0) center_x, center_y = self._get_beamcenter_position(detector_pix_nums_x, detector_pix_nums_y, offset_x, offset_y) # distance [cm] from the beam center on detector plane detector_ind_x = np.arange(detector_pix_nums_x) detector_ind_y = np.arange(detector_pix_nums_y) # shif 0.5 pixel so that pix position is at the center of the pixel detector_ind_x = detector_ind_x + 0.5 detector_ind_y = detector_ind_y + 0.5 # the relative postion from the beam center detector_ind_x = detector_ind_x - center_x detector_ind_y = detector_ind_y - center_y # unit correction in cm detector_ind_x = detector_ind_x * pix_x_size detector_ind_y = detector_ind_y * pix_y_size qx_value = np.zeros(len(detector_ind_x)) qy_value = np.zeros(len(detector_ind_y)) i = 0 for indx in detector_ind_x: qx_value[i] = self._get_qx(indx, sample2detector_distance, wavelength) i += 1 i = 0 for indy in detector_ind_y: qy_value[i] = self._get_qx(indy, sample2detector_distance, wavelength) i += 1 # qx_value and qy_value values in array qx_value = qx_value.repeat(detector_pix_nums_y) qx_value = qx_value.reshape(detector_pix_nums_x, detector_pix_nums_y) qy_value = qy_value.repeat(detector_pix_nums_x) qy_value = qy_value.reshape(detector_pix_nums_y, detector_pix_nums_x) qy_value = qy_value.transpose() # p min and max values among the center of pixels self.qx_min = np.min(qx_value) self.qx_max = np.max(qx_value) self.qy_min = np.min(qy_value) self.qy_max = np.max(qy_value) # Appr. min and max values of the detector display limits # i.e., edges of the last pixels. self.qy_min += self._get_qx(-0.5 * pix_y_size, sample2detector_distance, wavelength) self.qy_max += self._get_qx(0.5 * pix_y_size, sample2detector_distance, wavelength) #if self.qx_min == self.qx_max: self.qx_min += self._get_qx(-0.5 * pix_x_size, sample2detector_distance, wavelength) self.qx_max += self._get_qx(0.5 * pix_x_size, sample2detector_distance, wavelength) # min and max values of detecter self.detector_qx_min = self.qx_min self.detector_qx_max = self.qx_max self.detector_qy_min = self.qy_min self.detector_qy_max = self.qy_max # try to set it as a Data2D otherwise pass (not required for now) try: from sas.sascalc.dataloader.data_info import Data2D output = Data2D() inten = np.zeros_like(qx_value) output.data = inten output.qx_data = qx_value output.qy_data = qy_value except: logging.error(sys.exc_value) return output def _get_qx(self, dx_size, det_dist, wavelength): """ :param dx_size: x-distance from beam center [cm] :param det_dist: sample to detector distance [cm] :return: q-value at the given position """ # Distance from beam center in the plane of detector plane_dist = dx_size # full scattering angle on the x-axis theta = np.arctan(plane_dist / det_dist) qx_value = (2.0 * pi / wavelength) * np.sin(theta) return qx_value def _get_polar_value(self, qx_value, qy_value): """ Find qr_value and phi from qx_value and qy_value values : return qr_value, phi """ # find |q| on detector plane qr_value = sqrt(qx_value*qx_value + qy_value*qy_value) # find angle phi phi = self._atan_phi(qy_value, qx_value) return qr_value, phi def _get_beamcenter_position(self, num_x, num_y, offset_x, offset_y): """ :param num_x: number of pixel in x-direction :param num_y: number of pixel in y-direction :param offset: detector offset in x-direction in pix number :return: pix number; pos_x, pos_y in pix index """ # beam center position pos_x = num_x / 2 pos_y = num_y / 2 # correction for offset pos_x += offset_x # correction for gravity that is always negative pos_y -= offset_y return pos_x, pos_y def _get_beamcenter_drop(self): """ Get the beam center drop (delta y) in y diection due to gravity :return delta y: the beam center drop in cm """ # Check if mass == 0 (X-ray). if self.mass == 0: return 0 # Covert unit from A to cm unit_cm = 1e-08 # Velocity of neutron in horizontal direction (~ actual velocity) velocity = _PLANK_H / (self.mass * self.wave.wavelength * unit_cm) # Compute delta y delta_y = 0.5 delta_y *= _GRAVITY sampletodetector = self.sample2detector_distance[0] - \ self.sample2sample_distance[0] delta_y *= sampletodetector delta_y *= (self.source2sample_distance[0] + self.sample2detector_distance[0]) delta_y /= (velocity * velocity) return delta_y