Source code for xrayutilities.analysis.sample_align

# This file is part of xrayutilities.
#
# xrayutilities is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>.
#
# Copyright (c) 2011-2023 Dominik Kriegner <dominik.kriegner@gmail.com>

"""
functions to help with experimental alignment during experiments, especially
for experiments with linear and area detectors
"""

import glob
import math
import numbers
import re
import time

import numpy
from numpy import cos, degrees, radians, sin, tan
from scipy import odr, optimize
from scipy.ndimage import center_of_mass

from .. import config, cxrayutilities
from .. import math as xumath
from .. import utilities
from ..exception import InputError
from ..math import fwhm_exp

# regular expression to check goniometer circle syntax
circleSyntax = re.compile("[xyz][+-]")

#################################################
# channel per degree calculation
#################################################


[docs] def psd_chdeg(angles, channels, stdev=None, usetilt=True, plot=True, datap="xk", modelline="--r", modeltilt="-b", fignum=None, mlabel="fit", mtiltlabel="fit w/tilt", dlabel="data", figtitle=True): """ function to determine the channels per degree using a linear fit of the function nchannel = center_ch+chdeg*tan(angles) or the equivalent including a detector tilt Parameters ---------- angles : array-like detector angles for which the position of the beam was measured channels : array-like detector channels where the beam was found stdev : array-like, optional standard deviation of the beam position plot : bool, optional flag to specify if a visualization of the fit should be done usetilt : bool, optional whether to use model considering a detector tilt, i.e. deviation angle of the pixel direction from orthogonal to the primary beam (default: True) Other Parameters ---------------- datap : str, optional plot format of data points modelline : str, optional plot format of modelline modeltilt : str, optional plot format of modeltilt fignum : int or str, optional figure number to use for the plot mlabel : str label of the model w/o tilt to be used in the plot mtiltlabel : str label of the model with tilt to be used in the plot dlabel : str label of the data line to be used in the plot figtitle : bool flag to tell if the figure title should show the fit parameters Returns ------- pixelwidth : float the width of one detector channel @ 1m distance, which is negative in case the hit channel number decreases upon an increase of the detector angle. centerch : float center channel of the detector tilt : float tilt of the detector from perpendicular to the beam (will be zero in case of usetilt=False) Note: L/pixelwidth*pi/180 = channel/degree for large detector distance with the sample detector disctance L """ if stdev is None: stdevu = numpy.ones(len(channels)) else: stdevu = stdev # define detector model and other functions needed for the tilt def straight_tilt(p, x): """ model for straight-linear detectors including tilt Parameters ---------- p : list [L/w_pix*pi/180 ~= channel/degree, center_channel, detector_tilt] with L sample detector disctance, and w_pix the width of one detector channel x : array-like independent variable of the model: detector angle (degree) """ rad = radians(x) r = math.degrees(p[0]) * sin(rad) / \ cos(rad - math.radians(p[2])) + p[1] return r def straight_tilt_der_x(p, x): """ derivative of straight-linear detector model with respect to the angle for parameter description see straigt_tilt """ rad = radians(x) p2 = math.radians(p[2]) r = math.degrees(p[0]) * \ (cos(rad) / cos(rad - p2) + sin(rad) / cos(rad - p2) ** 2 * sin(rad - p2)) return r def straight_tilt_der_p(p, x): """ derivative of straight-linear detector model with respect to the paramters for parameter description see straigt_tilt """ rad = radians(x) p2 = math.radians(p[2]) r = numpy.concatenate([degrees(sin(rad)) / cos(rad - p2), numpy.ones(x.shape, dtype=float), - math.degrees(p[0]) * sin(rad) / cos(rad - p2) ** 2 * sin(rad - p2)]) r.shape = (3,) + x.shape return r # fit linear model = odr.unilinear data = odr.RealData(angles, channels, sy=stdevu) my_odr = odr.ODR(data, model) # fit type 2 for least squares my_odr.set_job(fit_type=2) fitlin = my_odr.run() # fit linear with tangens angle model = odr.unilinear data = odr.RealData(degrees(tan(radians(angles))), channels, sy=stdevu) my_odr = odr.ODR(data, model) # fit type 2 for least squares my_odr.set_job(fit_type=2) fittan = my_odr.run() if usetilt: # fit tilted straight detector model model = odr.Model(straight_tilt, fjacd=straight_tilt_der_x, fjacb=straight_tilt_der_p) data = odr.RealData(angles, channels, sy=stdevu) my_odr = odr.ODR(data, model, beta0=[fittan.beta[0], fittan.beta[1], 0]) # fit type 2 for least squares my_odr.set_job(fit_type=2) fittilt = my_odr.run() if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis.psd_chdeg') if plot: markersize = 6.0 markeredgewidth = 1.5 linewidth = 2.0 if fignum is None: plt.figure() else: plt.figure(fignum) # first plot to show linear model ax1 = plt.subplot(211) angr = angles.max() - angles.min() angp = numpy.linspace(angles.min() - angr * 0.1, angles.max() + angr * .1, 1000) if modelline: plt.plot(angp, odr.unilinear.fcn(fittan.beta, degrees(tan(radians(angp)))), modelline, label=mlabel, lw=linewidth) plt.plot(angp, odr.unilinear.fcn(fitlin.beta, angp), '-k', label='') if usetilt: plt.plot(angp, straight_tilt(fittilt.beta, angp), modeltilt, label=mtiltlabel, lw=linewidth) if stdev is None: plt.plot(angles, channels, datap, ms=markersize, mew=markeredgewidth, mec=datap[-1], mfc='none', label=dlabel) else: plt.errorbar(angles, channels, fmt=datap, yerr=stdevu, ms=markersize, mew=markeredgewidth, mec=datap[-1], mfc='none', label=dlabel, ecolor='0.5') plt.grid(True) leg = plt.legend(numpoints=1) leg.get_frame().set_alpha(0.8) plt.ylabel("channel number") # lower plot to show deviations from linear model plt.subplot(212, sharex=ax1) if modelline: plt.plot(angp, odr.unilinear.fcn(fittan.beta, degrees(tan(radians(angp)))) - odr.unilinear.fcn(fitlin.beta, angp), modelline, label=mlabel, lw=linewidth) if usetilt: plt.plot(angp, straight_tilt(fittilt.beta, angp) - odr.unilinear.fcn(fitlin.beta, angp), modeltilt, label=mtiltlabel, lw=linewidth) if stdev is None: plt.plot(angles, channels - odr.unilinear.fcn(fitlin.beta, angles), datap, ms=markersize, mew=markeredgewidth, mec=datap[-1], mfc='none', label=dlabel) else: plt.errorbar(angles, channels - odr.unilinear.fcn(fitlin.beta, angles), fmt=datap, yerr=stdevu, ms=markersize, mew=markeredgewidth, mec=datap[-1], mfc='none', label=dlabel, ecolor='0.5') plt.xlabel("detector angle (deg)") plt.ylabel("ch. num. - linear trend") plt.grid(True) plt.hlines(0, angp.min(), angp.max()) plt.tight_layout() if figtitle: if usetilt: plt.suptitle( "L/w*pi/180: %8.2f; center channel: %8.2f; tilt: %5.2fdeg" % (fittilt.beta[0], fittilt.beta[1], fittilt.beta[2])) else: plt.suptitle("L/w*pi/180: %8.2f; center channel: %8.2f" % (fittan.beta[0], fittan.beta[1])) if usetilt: fit = fittilt else: fit = fittan if config.VERBOSITY >= config.INFO_LOW: if usetilt: print("XU.analysis.psd_chdeg: channelwidth@1m / center channel /" " tilt: %8.4e / %8.2f / %6.3fdeg" % (abs(1 / math.degrees(fit.beta[0])), fit.beta[1], fit.beta[2])) print("XU.analysis.psd_chdeg: error of channelwidth / " "center channel / tilt: %8.4e / %8.3f / %6.3fdeg" % (math.radians(fit.sd_beta[0] / fit.beta[0] ** 2), fit.sd_beta[1], fit.sd_beta[2])) else: print("XU.analysis.psd_chdeg: channelwidth@1m / center channel: " "%8.4e / %8.2f" % (1 / math.degrees(fit.beta[0]), fit.beta[1])) print("XU.analysis.psd_chdeg: error of channelwidth / " "center channel: %8.4e / %8.3f" % (math.radians(fit.sd_beta[0] / fit.beta[0] ** 2), fit.sd_beta[1])) if usetilt: return (1. / math.degrees(fit.beta[0]), fit.beta[1], fit.beta[2]) return (1. / math.degrees(fit.beta[0]), fit.beta[1], 0.)
################################################# # channel per degree calculation from scan with # linear detector (determined maximum by peak_fit) #################################################
[docs] def linear_detector_calib(angle, mca_spectra, **keyargs): """ function to calibrate the detector distance/channel per degrees for a straight linear detector mounted on a detector arm Parameters ---------- angle : array-like array of angles in degree of measured detector spectra mca_spectra : array-like corresponding detector spectra (shape: (len(angle), Nchannels) r_i : str, optional primary beam direction as vector [xyz][+-]; default: 'y+' detaxis : str, optional detector arm rotation axis [xyz][+-]; default: 'x+' Other parameters ---------------- plot : bool flag to specify if a visualization of the fit should be done usetilt : bool whether to use model considering a detector tilt, i.e. deviation angle of the pixel direction from orthogonal to the primary beam (default: True) Returns ------- pixelwidth : float width of the pixel at one meter distance, pixelwidth is negative in case the hit channel number decreases upon an increase of the detector angle center_channel : float central channel of the detector detector_tilt : float, optional if usetilt=True the fitted tilt of the detector is also returned Note: L/pixelwidth*pi/180 ~= channel/degree, with the sample detector distance L The function also prints out how a linear detector can be initialized using the results obtained from this calibration. Carefully check the results See Also -------- psd_chdeg : low level function with more configurable options """ if "detaxis" in keyargs: detrotaxis = keyargs["detaxis"] keyargs.pop("detaxis") else: # use default detrotaxis = 'x+' if "r_i" in keyargs: r_i = keyargs["r_i"] keyargs.pop("r_i") else: # use default r_i = 'y+' # max intensity per spectrum mca_int = mca_spectra.sum(axis=1) mca_avg = numpy.average(mca_int) mca_rowmax = numpy.max(mca_int) mca_std = numpy.std(mca_int) # determine positions pos = [] posstd = [] ang = [] nignored = 0 for i in range(len(mca_spectra)): row = mca_spectra[i, :] row_int = row.sum() if ((abs(row_int - mca_avg) > 3 * mca_std) or (row_int - mca_rowmax * 0.7 < 0)): if config.VERBOSITY >= config.DEBUG: print("XU.analysis.linear_detector_calib: spectrum #%d " "out of intensity range -> ignored" % i) nignored += 1 continue maxp = numpy.argmax(row) fwhm = fwhm_exp(numpy.arange(row.size), row) N = int(7 * numpy.ceil(fwhm)) // 2 * 2 # fit beam position # determine maximal usable length of array around peak position Nuse = min(maxp + N // 2, len(row) - 1) - max(maxp - N // 2, 0) param, perr, _ = xumath.peak_fit( numpy.arange(Nuse), row[max(maxp - N // 2, 0):min(maxp + N // 2, len(row) - 1)], peaktype='PseudoVoigt') if param[0] > 0 and param[0] < Nuse and perr[0] < Nuse/2.: param[0] += max(maxp - N // 2, 0) pos.append(param[0]) posstd.append(perr[0]) ang.append(angle[i]) ang = numpy.array(ang) pos = numpy.array(pos) posstd = numpy.array(posstd) if config.VERBOSITY >= config.INFO_ALL: print("XU.analysis.linear_detector_calib: using %d out of %d given " "spectra." % (len(ang), len(angle))) if config.VERBOSITY >= config.DEBUG: print("XU.analysis.linear_detector_calib: determined peak positions:") print(zip(pos, posstd)) detparam = psd_chdeg(ang, pos, stdev=posstd, **keyargs) if numpy.sign(detparam[0]) > 0: sign = '-' else: sign = '+' detaxis = ' ' detd = numpy.cross(xumath.getVector(detrotaxis), xumath.getVector(r_i)) argm = numpy.abs(detd).argmax() def flipsign(char, val): if numpy.sign(val) < 0: if char == '+': return '-' return '+' return char if argm == 0: detaxis = 'x' + flipsign(sign, detd[argm]) elif argm == 1: detaxis = 'y' + flipsign(sign, detd[argm]) elif argm == 2: detaxis = 'z' + flipsign(sign, detd[argm]) if config.VERBOSITY >= config.INFO_LOW: print("XU.analysis.linear_detector_calib:\n\tused/total spectra: " f"{mca_spectra.shape[0] - nignored:d}/{mca_spectra.shape[0]:d}") print("\tdetector rotation axis (given by user/default input): " f"{detrotaxis}") if len(detparam) == 3: tilt = detparam[2] else: tilt = 0 print(f"\tdetector initialization with: init_linear('{detaxis}', " f"{abs(detparam[1]):.2f}, {mca_spectra.shape[1]:d}, " f"pixelwidth={abs(detparam[0]):.4e}, distance=1., " f"tilt={tilt:.2f})") return detparam
###################################################### # detector parameter calculation from scan with # area detector (determine maximum by center of mass) ######################################################
[docs] def area_detector_calib(angle1, angle2, ccdimages, detaxis, r_i, plot=True, cut_off=0.7, start=(None, None, 1, 0, 0, 0, 0), fix=(False, False, True, False, False, False, False), fig=None, wl=None, plotlog=False, nwindow=50, debug=False): """ function to calibrate the detector parameters of an area detector it determines the detector tilt possible rotations and offsets in the detector arm angles Parameters ---------- angle1 : array-like outer detector arm angle angle2 : array-like inner detector arm angle ccdimages : array-like images of the ccd taken at the angles given above detaxis : list of str detector arm rotation axis; default: ['z+', 'y-'] r_i : str primary beam direction [xyz][+-]; default 'x+' plot : bool, optional flag to determine if results and intermediate results should be plotted; default: True cut_off : float, optional cut off intensity to decide if image is used for the determination or not; default: 0.7 = 70% start : tuple, optional sequence of start values of the fit for parameters, which can not be estimated automatically or might want to be fixed. These are: pwidth1, pwidth2, distance, tiltazimuth, tilt, detector_rotation, outerangle_offset. By default (None, None, 1, 0, 0, 0, 0) is used. fix : tuple of bool fix parameters of start (default: (False, False, True, False, False, False, False)) It is strongly recommended to either fix the distance or the pwidth1, 2 values. fig : Figure, optional matplotlib figure used for plotting the error default: None (creates own figure) wl : float or str wavelength of the experiment in angstrom (default: config.WAVELENGTH) value does not really matter here but does affect the scaling of the error plotlog : bool flag to specify if the created error plot should be on log-scale nwindow : int window size for determination of the center of mass position after the center of mass of every full image is determined, the center of mass is determined again using a window of size nwindow in order to reduce the effect of hot pixels. debug : bool flag to specify that you want to see verbose output and saving of images to show if the CEN determination works """ if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis.area_' 'detector_calib') if wl is None: wl = config.WAVELENGTH else: wl = utilities.wavelength(wl) for i in [0, 1]: if fix[i] and not numpy.isscalar(start[i]): raise ValueError("XU.analysis.area_detector_calib: start value for" " pwidth%d must be given if it should be fixed " "during the fit" % (i+1)) t0 = time.time() Npoints = len(angle1) if debug: print("number of given images: %d" % Npoints) # determine center of mass position from detector images # also use only images with an intensity larger than 70% of the average # intensity n1 = numpy.zeros(0, dtype=numpy.double) n2 = n1 ang1 = n1 ang2 = n1 avg = 0 for i in range(Npoints): avg += numpy.sum(ccdimages[i]) avg /= float(Npoints) (N1, N2) = ccdimages[0].shape if debug: print(f"average intensity per image: {avg:.1f}") for i in range(Npoints): if debug and i == 0: print("angle1, angle2, cen1, cen2") img = ccdimages[i] if numpy.sum(img) > cut_off * avg: cen1, cen2 = _peak_position(img, nwindow, plot=debug and plot) n1 = numpy.append(n1, cen1) n2 = numpy.append(n2, cen2) ang1 = numpy.append(ang1, angle1[i]) ang2 = numpy.append(ang2, angle2[i]) if debug: print(f"{angle1[i]:8.3f} {angle2[i]:8.3f} \t" f"{cen1:.2f} {cen2:.2f}") Nused = len(ang1) if debug: print("Nused / Npoints: %d / %d" % (Nused, Npoints)) # determine detector directions detdir1, detdir2 = _determine_detdir( ang1 - start[6], ang2, n1, n2, detaxis, r_i) if debug: print(f"determined detector directions:[{detdir1}, {detdir2}]") epslist = [] paramlist = [] epsmin = numpy.inf print("tiltaz tilt detrot offset: error (relative) (fittime)") print("------------------------------------------------------------") # find optimal detector rotation (however keep other parameters free) detrot = start[5] if not fix[5]: for detrotstart in numpy.linspace(start[5] - 1, start[5] + 1, 40): start = start[:5] + (detrotstart,) + (start[6],) eps, param, fit = _area_detector_calib_fit( ang1, ang2, n1, n2, detaxis, r_i, detdir1, detdir2, start=start, fix=fix, full_output=True, wl=wl, debug=debug) epslist.append(eps) paramlist.append(param) if epslist[-1] < epsmin: epsmin = epslist[-1] parammin = param fitmin = fit detrot = param[7] if debug: print(eps, param) Ntiltaz = 1 if fix[3] else 5 Ntilt = 1 if fix[4] else 6 Noffset = 1 if fix[6] else 100 if fix[6]: Ntilt = Ntilt * 8 if not fix[4] else Ntilt Ntiltaz = Ntiltaz * 7 if not fix[3] else Ntiltaz startparam = start[:5] + (detrot,) + (start[6],) if debug: print(f"start params: {str(startparam)}") Ntot = Ntiltaz * Ntilt * Noffset ict = 0 for tiltazimuth in numpy.linspace(startparam[3] if fix[3] else 0, 360, Ntiltaz, endpoint=False): for tilt in numpy.linspace( startparam[4] if fix[4] else max(0, startparam[4]-2), max(0, startparam[4]-2)+4, Ntilt): for offset in numpy.linspace( startparam[6] if fix[6] else - 3 + startparam[6], 3 + startparam[6], Noffset): t1 = time.time() start = (startparam[0], startparam[1], startparam[2], tiltazimuth, tilt, startparam[5], offset) eps, param, fit = _area_detector_calib_fit( ang1, ang2, n1, n2, detaxis, r_i, detdir1, detdir2, start=start, fix=fix, full_output=True, wl=wl) epslist.append(eps) paramlist.append(param) t2 = time.time() print("%d/%d\t%6.1f %6.2f %8.3f %8.3f: %10.4e (%4.2f) " "(%5.2fsec)" % ((ict, Ntot) + start[3:] + (epslist[-1], epslist[-1] / epsmin, t2 - t1))) ict += 1 if epslist[-1] < epsmin: print("************************") print("new global minimum found") epsmin = epslist[-1] parammin = param fitmin = fit print("new best parameters: %.2f %.2f %10.4e %10.4e %8.4f " "%.1f %.2f %.3f %.3f" % parammin) print("************************\n") (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset) = parammin if plot: if fig: plt.figure(fig.number) else: plt.figure("CCD Calib fit") plt.clf() nparams = numpy.array(paramlist) neps = numpy.array(epslist) labels = ( 'cch1 (1)', 'cch2 (1)', r'pwidth1 ($\mu$m)', r'pwidth2 ($\mu$m)', 'distance (m)', 'tiltazimuth (deg)', 'tilt (deg)', 'detrot (deg)', 'outerangle offset (deg)') xscale = (1., 1., 1.e6, 1.e6, 1., 1., 1., 1., 1.) for p in range(9): ax = plt.subplot(3, 3, p + 1) if plotlog: plt.semilogy(nparams[:, p] * xscale[p], neps, '.k') else: plt.scatter(nparams[:, p] * xscale[p], neps, c=nparams[:, -1], s=10, marker='o', cmap="gnuplot", edgecolor='none') plt.xlabel(labels[p]) if plotlog: plt.semilogy(parammin[p] * xscale[p], epsmin, 'ok', ms=8, mew=2.5, mec='k', mfc='w') else: plt.plot(parammin[p] * xscale[p], epsmin, 'ok', ms=8, mew=2.5, mec='k', mfc='w') plt.ylim(epsmin * 0.7, epsmin * 2.) plt.locator_params(nbins=4, axis='x') if p > 1: if fix[p-2]: ax.set_facecolor('0.85') plt.tight_layout() if config.VERBOSITY >= config.INFO_LOW: print(f"total time needed for fit: {time.time() - t0:.2f}sec") print("fitted parameters: epsilon: %10.4e (%d,%s) " # pylint: disable-next=no-member % (epsmin, fitmin.info, repr(fitmin.stopreason))) print("param: (cch1, cch2, pwidth1, pwidth2, tiltazimuth, tilt, " "detrot, outerangle_offset)") print("param: %.2f %.2f %10.4e %10.4e %.4f %.1f %.2f %.3f %.3f" % (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset)) if config.VERBOSITY > 0: print("please check the resulting data (consider setting plot=True)") print("detector rotation axis / primary beam direction " "(given by user): %s / %s" % (repr(detaxis), r_i)) print(f"detector pixel directions / distance: {detdir1} {detdir2} / " f"{1.0:g}") print("\tdetector initialization with: init_area('%s', '%s', " "cch1=%.2f, cch2=%.2f, Nch1=%d, Nch2=%d, pwidth1=%.4e, " "pwidth2=%.4e, distance=%.5f, detrot=%.3f, tiltazimuth=%.1f, " "tilt=%.3f)" % (detdir1, detdir2, cch1, cch2, N1, N2, pwidth1, pwidth2, distance, detrot, tiltazimuth, tilt)) print("AND ALWAYS USE an (additional) OFFSET of " f"{outerangle_offset:.4f}deg in the OUTER DETECTOR ANGLE!") return (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset), epsmin
def _peak_position(img, nwindow, plot=False): """ function to determine the peak position on the detector using the center of mass (COM) Parameters ---------- img : array-like detector image data as 2D array nwindow : int to avoid influence of hot pixels far away from the peak position the center of mass approach is repeated with a window around the COM of the full image. COM of the size (nwindow, nwindow) is returned plot : bool, optional the result of the of the determination can be saved as a plot """ nw = nwindow // 2 [cen1r, cen2r] = center_of_mass(img) for i in range(11): # refine center of mass multiple times [cen1, cen2] = center_of_mass( img[max(int(cen1r) - nw, 0): min(int(cen1r) + nw, img.shape[0]), max(int(cen2r) - nw, 0): min(int(cen2r) + nw, img.shape[1])]) cen1 += max(int(cen1r) - nw, 0) cen2 += max(int(cen2r) - nw, 0) if numpy.linalg.norm((cen1 - cen1r, cen2 - cen2r)) > 3: cen1r, cen2r = (cen1, cen2) else: break if i == 10 and config.VERBOSITY >= config.INFO_LOW: print("XU.analysis._peak_position: Warning: peak position " "determination not converged, consider debug mode!") if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis._peak_' 'position') if plot: plt.figure("_ccd") plt.imshow(utilities.maplog(img), origin='lower') plt.plot(cen2, cen1, 'ow', mfc='none') plt.axis([cen2 - nw, cen2 + nw, cen1 - nw, cen1 + nw]) plt.colorbar() fnr = len(glob.glob('xu_calib_ccd_img*.png')) plt.savefig("xu_calib_ccd_img%d.png" % (fnr + 1)) plt.close("_ccd") return cen1, cen2 def _determine_detdir(ang1, ang2, n1, n2, detaxis, r_i): """ determines detector pixel direction from correlation analysis of linear fits to the observed pixel numbers of the primary beam. """ # center channel and detector pixel direction and pixel size (s1, _), r1 = xumath.linregress(ang1, n1) (s2, _), r2 = xumath.linregress(ang1, n2) (s3, _), r3 = xumath.linregress(ang2, n1) (s4, _), r4 = xumath.linregress(ang2, n2) # determine detector directions s = ord('x') + ord('y') + ord('z') c1 = ord(detaxis[0][0]) + ord(r_i[0]) c2 = ord(detaxis[1][0]) + ord(r_i[0]) sign1 = numpy.sign(numpy.sum(numpy.cross(xumath.getVector(detaxis[0]), xumath.getVector(r_i)))) sign2 = numpy.sign(numpy.sum(numpy.cross(xumath.getVector(detaxis[1]), xumath.getVector(r_i)))) if ((r1 + r2 > r3 + r4) and r1 > r2) or ((r1 + r2 < r3 + r4) and r3 < r4): detdir1 = chr(s - c1) detdir2 = chr(s - c2) if numpy.sign(s1) > 0: if sign1 > 0: detdir1 += '-' else: detdir1 += '+' else: if sign1 > 0: detdir1 += '+' else: detdir1 += '-' if numpy.sign(s4) > 0: if sign2 > 0: detdir2 += '-' else: detdir2 += '+' else: if sign2 > 0: detdir2 += '+' else: detdir2 += '-' else: detdir1 = chr(s - c2) detdir2 = chr(s - c1) if numpy.sign(s3) > 0: if sign2 > 0: detdir1 += '-' else: detdir1 += '+' else: if sign2 > 0: detdir1 += '+' else: detdir1 += '-' if numpy.sign(s2) > 0: if sign1 > 0: detdir2 += '-' else: detdir2 += '+' else: if sign1 > 0: detdir2 += '+' else: detdir2 += '-' return detdir1, detdir2 def _area_detector_calib_fit(ang1, ang2, n1, n2, detaxis, r_i, detdir1, detdir2, start=(None, None, 1, 0, 0, 0, 0), fix=(False, False, True, False, False, False, False), full_output=False, wl=1., debug=False): """ INTERNAL FUNCTION function to calibrate the detector parameters of an area detector it determines the detector tilt possible rotations and offsets in the detector arm angles parameters ---------- angle1 : array-like outer detector arm angle angle2 : array-like inner detector arm angle n1, n2 : int pixel number at which the primary beam was observed detaxis : list of str detector arm rotation axis; default: ['z+', 'y-'] r_i : str primary beam direction [xyz][+-]; default 'x+' detdir1, detdir2 : str detector pixel directions of first and second pixel coordinates; e.g. 'y+' start : tuple, optional sequence of start values of the fit for parameters, which can not be estimated automatically or might want to be fixed. These are: pwidth1, pwidth2, distance, tiltazimuth, tilt, detector_rotation, outerangle_offset. By default (None, None, 1, 0, 0, 0, 0) is used. fix : tuple of bool fix parameters of start (default: (False, False, True, False, False, False, False)) It is strongly recommended to either fix the distance or the pwidth1, 2 values. fig : Figure, optional matplotlib figure used for plotting the error default: None (creates own figure) full_output : bool flag to tell if to return fit object with final parameters and detector directions wl : float or str wavelength of the experiment in angstrom (default: 1) value does not really matter here but does affect the scaling of the error debug : bool flag to tell if you want to see debug output of the script (switch this to true only if you can handle it :)) Returns ------- float final epsilon of the fit param : list, optional if full_output: fit parameters fit : object, optional if full_output: fit object """ def areapixel(params, detectorDir1, detectorDir2, r_i, detectorAxis, *args, **kwargs): """ angular to momentum space conversion for pixels of an area detector the center pixel is in direction of self.r_i when detector angles are zero the detector geometry must be given to the routine Parameters ---------- params : list parameters of the detector calibration model (cch1, cch2, pwidth1, pwidth2, tiltazimuth, tilt, detrot) - cch1, cch2: center pixel, in direction of self.r_i at zero detectorAngles; - pwidth1, pwidth2: width of one pixel (same unit as distance); - distance: distance of center pixel from center of rotation; - tiltazimuth: direction of the tilt vector in the detector plane (in degree); - tilt: tilt of the detector plane around an axis normal to the direction given by the tiltazimuth; - detrot: detector rotation around the primary beam direction as given by r_i; detectorDir1 : str direction of the detector (along the pixel direction 1); e.g. 'z+' means higher pixel numbers at larger z positions detectorDir2 : str direction of the detector (along the pixel direction 2); e.g. 'x+' r_i : str primary beam direction e.g. 'x+' detectorAxis : list or tuple detector circles e.g. ['z+', 'y-'] would mean a detector arm with a two rotations *args : array-like detector angles and channel numbers; *dAngles* as numpy array, lists or Scalars in total len(detectorAxis) must be given starting with the outer most circle. All arguments must have the same shape or length. *channel numbers* n1 and n2 where the primary beam hits the detector with same length as the detector values delta : list, optional giving delta angles to correct the given ones for misalignment delta must be an numpy array or list of len(*dAngles) used angles are than *args - delta wl : float or str, optional x-ray wavelength in angstrom (default: 1 (since it does not matter here)) deg : bool, optional flag to tell if angles are passed as degree (default: True) Returns ------- ndarray reciprocal space position of detector pixels n1, n2 in a numpy.ndarray of shape ( len(args) , 3 ) """ # check detector circle argument if isinstance(detectorAxis, (str, list, tuple)): if isinstance(detectorAxis, str): dAxis = list([detectorAxis]) else: dAxis = list(detectorAxis) for circ in dAxis: if not isinstance(circ, str) or len(circ) != 2: raise InputError("QConversionPixel: incorrect detector " "circle type or syntax (%s)" % repr(circ)) if not circleSyntax.search(circ): raise InputError("QConversionPixel: incorrect detector " "circle syntax (%s)" % circ) else: raise TypeError("Qconversion error: invalid type for detectorAxis," " must be str, list or tuple") # add detector rotation around primary beam dAxis += [r_i] _detectorAxis_str = '' for circ in dAxis: _detectorAxis_str += circ Nd = len(dAxis) Nargs = Nd + 2 - 1 # check detectorDir if not isinstance(detectorDir1, str) or len(detectorDir1) != 2: raise InputError("QConversionPixel: incorrect detector direction1 " "type or syntax (%s)" % repr(detectorDir1)) if not circleSyntax.search(detectorDir1): raise InputError("QConversionPixel: incorrect detector direction1 " "syntax (%s)" % detectorDir1) _area_detdir1 = detectorDir1 if not isinstance(detectorDir2, str) or len(detectorDir2) != 2: raise InputError("QConversionPixel: incorrect detector direction2 " "type or syntax (%s)" % repr(detectorDir2)) if not circleSyntax.search(detectorDir2): raise InputError("QConversionPixel: incorrect detector direction2 " "syntax (%s)" % detectorDir2) _area_detdir2 = detectorDir2 # parse parameter arguments _area_cch1 = float(params[0]) _area_cch2 = float(params[1]) _area_pwidth1 = float(params[2]) _area_pwidth2 = float(params[3]) _area_distance = float(params[4]) _area_tiltazimuth = math.radians(params[5]) _area_tilt = math.radians(params[6]) _area_rot = float(params[7]) _area_ri = xumath.getVector(r_i) * _area_distance # kwargs wl = utilities.wavelength(kwargs.get('wl', 1.)) deg = kwargs.get('deg', True) delta = numpy.asarray(kwargs.get('delta', numpy.zeros(Nd)), dtype=numpy.double) if delta.size != Nd - 1: raise InputError("QConversionPixel: keyword argument delta " "does not have an appropriate shape") # prepare angular arrays from *args # need one sample angle and one detector angle array if len(args) != Nargs: raise InputError("QConversionPixel: wrong amount (%d) of arguments" " given, number of arguments should be %d" % (len(args), Nargs)) try: Npoints = len(args[0]) except (TypeError, IndexError): Npoints = 1 dAngles = numpy.array((), dtype=numpy.double) for i in range(Nd - 1): arg = args[i] if not isinstance(arg, (numbers.Number, tuple, list, numpy.ndarray)): raise TypeError("QConversionPixel: invalid type for one of " "the detector coordinates, must be scalar, " "list or array") if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) arg = arg - delta[i] dAngles = numpy.concatenate((dAngles, arg)) # add detector rotation around primary beam dAngles = numpy.concatenate((dAngles, numpy.ones(arg.shape, dtype=numpy.double) * _area_rot)) # read channel numbers n1 = numpy.array((), dtype=numpy.double) n2 = numpy.array((), dtype=numpy.double) arg = args[Nd - 1] if not isinstance(arg, (numbers.Number, tuple, list, numpy.ndarray)): raise TypeError("QConversionPixel: invalid type for one of the " "detector coordinates, must be scalar, list or " "array") if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) n1 = arg arg = args[Nd] if not isinstance(arg, (numbers.Number, tuple, list, numpy.ndarray)): raise TypeError("QConversionPixel: invalid type for one of the " "detector coordinates, must be scalar, list or " "array") if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) n2 = arg dAngles.shape = (Nd, Npoints) dAngles = dAngles.transpose() if deg: dAngles = radians(dAngles) qpos = cxrayutilities.ang2q_conversion_area_pixel( dAngles, n1, n2, _area_ri, _detectorAxis_str, _area_cch1, _area_cch2, _area_pwidth1, _area_pwidth2, _area_detdir1, _area_detdir2, _area_tiltazimuth, _area_tilt, wl, config.NTHREADS) return qpos[:, 0], qpos[:, 1], qpos[:, 2] def afunc(param, x, detectorDir1, detectorDir2, r_i, detectorAxis, wl): """ function for fitting the detector parameters basically this is a wrapper for the areapixel function parameters ---------- param : list fit parameters (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset) x : array-like independent variables (angle1, angle2, n1, n2) with shape (4, Npoints) detectorDir1 : str direction of the detector (along the pixel direction 1); e.g. 'z+' means higher pixel numbers at larger z positions detectorDir2 : str direction of the detector (along the pixel direction 2); e.g. 'x+' r_i : str primary beam direction e.g. 'x+' detectorAxis : list or tuple detector circles e.g. ['z+', 'y-'] would mean a detector arm with a two rotations wl : float or str wavelength of the experiment in angstrom Returns ------- ndarray reciprocal space position of detector pixels n1, n2 in a numpy.ndarray of shape (3, x.shape[1]) """ angle1 = x[0, :] angle2 = x[1, :] n1 = x[2, :] n2 = x[3, :] # use only positive tilt param[6] = abs(param[6]) (qx, qy, qz) = areapixel( param[:-1], detectorDir1, detectorDir2, r_i, detectorAxis, angle1, angle2, n1, n2, delta=[param[-1], 0.], wl=wl) return qx ** 2 + qy ** 2 + qz ** 2 Npoints = len(ang1) # guess initial parameters # center channel and detector pixel direction and pixel size (s1, i1), r1 = xumath.linregress(ang1 - start[6], n1) (s2, i2), r2 = xumath.linregress(ang1 - start[6], n2) (s3, i3), r3 = xumath.linregress(ang2, n1) (s4, i4), r4 = xumath.linregress(ang2, n2) distance = start[2] if ((r1 + r2 > r3 + r4) and r1 > r2) or ((r1 + r2 < r3 + r4) and r3 < r4): cch1 = i1 cch2 = i4 pwidth1 = 2 * distance / numpy.abs(s1) * math.tan(math.radians(0.5)) pwidth2 = 2 * distance / numpy.abs(s4) * math.tan(math.radians(0.5)) else: cch1 = i3 cch2 = i2 pwidth1 = 2 * distance / numpy.abs(s3) * math.tan(math.radians(0.5)) pwidth2 = 2 * distance / numpy.abs(s2) * math.tan(math.radians(0.5)) if numpy.isscalar(start[0]): pwidth1 = start[0] if numpy.isscalar(start[1]): pwidth2 = start[1] if numpy.isscalar(start[0]) or numpy.isscalar(start[1]): # find biggest correlation and recalculate distance idxmax = numpy.argmax((r1, r2, r3, r4)) s = (s1, s2, s3, s4)[idxmax] distance = abs(s) / math.tan(math.radians(0.5)) / 2 distance *= pwidth1 if idxmax < 2 else pwidth2 tilt = abs(start[4]) tiltazimuth = start[3] detrot = start[5] outerangle_offset = start[6] # parameters for the fitting param = (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset) if debug: print("initial parameters: ") print("primary beam / detector pixel directions / distance: " "%s / %s %s / %e" % (r_i, detdir1, detdir2, distance)) print("param: (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, " "tilt, detrot, outerangle_offset)") print("param: %.2f %.2f %10.4e %10.4e %.3f %.1f %.2f %.3f %.3f" % param) # set data x = numpy.empty((4, Npoints), dtype=numpy.double) x[0, :] = ang1 x[1, :] = ang2 x[2, :] = n1 x[3, :] = n2 data = odr.Data(x, y=1) # define model for fitting model = odr.Model(afunc, extra_args=(detdir1, detdir2, r_i, detaxis, wl), implicit=True) # check if parameters need to be fixed ifixb = () for i in range(len(fix)): ifixb += (int(not fix[i]),) my_odr = odr.ODR(data, model, beta0=param, ifixb=(1, 1) + ifixb, ifixx=(0, 0, 0, 0), stpb=(0.4, 0.4, pwidth1/50., pwidth2/50., distance/1000, 2, 0.125, 0.01, 0.01), sclb=(1/abs(cch1), 1/abs(cch2), 1/pwidth1, 1/pwidth2, 1/distance, 1/90., 1/0.2, 1/0.2, 1/0.2), maxit=200, ndigit=12, sstol=1e-11, partol=1e-11) if debug: my_odr.set_iprint(final=1) my_odr.set_iprint(iter=2) fit = my_odr.run() (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset) = fit.beta # fix things in parameters tiltazimuth = tiltazimuth % 360. tilt = abs(tilt) final_q = afunc([cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset], x, detdir1, detdir2, r_i, detaxis, wl) final_error = numpy.mean(final_q) if debug: # pylint: disable-next=no-member print("fitted parameters: (%e, %d, %s) " % (final_error, fit.info, repr(fit.stopreason))) print("primary beam / detector pixel directions / distance: " "%s / %s %s / %e" % (r_i, detdir1, detdir2, distance)) print("param: (cch1, cch2, pwidth1, pwidth2, disance, tiltazimuth, " "tilt, detrot, outerangle_offset)") print("param: %.2f %.2f %10.4e %10.4e %.3f %.1f %.2f %.3f %.3f" % (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset)) if full_output: return final_error, (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset), fit return final_error # ##################################################### # detector parameter calculation from scan with # area detector (determine maximum by center of mass) # #####################################################
[docs] def area_detector_calib_hkl(sampleang, angle1, angle2, ccdimages, hkls, experiment, material, detaxis, r_i, plot=True, cut_off=0.7, start=(None, None, 1, 0, 0, 0, 0, 0, 0, 'config'), fix=(False, False, True, False, False, False, False, False, False, False), fig=None, plotlog=False, nwindow=50, debug=False): """ function to calibrate the detector parameters of an area detector it determines the detector tilt possible rotations and offsets in the detector arm angles in this variant not only scans through the primary beam but also scans at a set of symmetric reflections can be used for the detector parameter determination. for this not only the detector parameters but in addition the sample orientation and wavelength need to be fit. Both images from the primary beam hkl = (0, 0, 0) and from a symmetric reflection hkl = (h, k, l) need to be given for a successful run. Parameters ---------- sampleang : array-like sample rocking angle (needed to align the reflections (same rotation direction as inner detector rotation)) other sample angle are not allowed to be changed during the scans angle1 : array-like outer detector arm angle angle2 : array-like inner detector arm angle ccdimages : array-like images of the ccd taken at the angles given above hkls : list or array-like hkl values for every image experiment : Experiment Experiment class object needed to get the UB matrix for the hkl peak treatment material : Crystal material used as reference crystal detaxis : list of str detector arm rotation axis; default: ['z+', 'y-'] r_i : str primary beam direction [xyz][+-]; default 'x+' plot : bool, optional flag to determine if results and intermediate results should be plotted; default: True cut_off : float, optional cut off intensity to decide if image is used for the determination or not; default: 0.7 = 70% start : tuple, optional sequence of start values of the fit for parameters, which can not be estimated automatically or might want to be fixed. These are: pwidth1, pwidth2, distance, tiltazimuth, tilt, detector_rotation, outerangle_offset, sampletilt, sampletiltazimuth, wavelength. By default (None, None, 1, 0, 0, 0, 0, 0, 0, 'config'). fix : tuple of bool fix parameters of start (default: (False, False, True, False, False, False, False, False, False, False)) It is strongly recommended to either fix the distance or the pwidth1, 2 values. fig : Figure, optional matplotlib figure used for plotting the error default: None (creates own figure) plotlog : bool flag to specify if the created error plot should be on log-scale nwindow : int window size for determination of the center of mass position after the center of mass of every full image is determined, the center of mass is determined again using a window of size nwindow in order to reduce the effect of hot pixels. debug : bool flag to specify that you want to see verbose output and saving of images to show if the CEN determination works """ if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis.area_' 'detector_calib_hkl') if start[-1] == 'config': start[-1] = config.WAVELENGTH elif isinstance(start[-1], str): start[-1] = utilities.wavelength(start[-1]) t0 = time.time() Npoints = len(angle1) if debug: print("number of given images: %d" % Npoints) # determine center of mass position from detector images also use only # images with an intensity larger than cut_off of the average intensity. # the image selection is only performed for images in the primary beam n1 = numpy.zeros(0, dtype=numpy.double) n2 = n1 ang1 = n1 ang2 = n1 sang = n1 usedhkls = [] avg = 0 imgpbcnt = 0 for i in range(Npoints): if numpy.all(hkls[i] == (0, 0, 0)): avg += numpy.sum(ccdimages[i]) imgpbcnt += 1 if imgpbcnt > 0: avg /= float(imgpbcnt) else: raise ValueError("XU.analyis.area_detector_calib_hkl: no required " "images in the primary beam given!") (N1, N2) = ccdimages[0].shape if debug: print(f"average intensity per image in the primary beam: {avg:.1f}") for i in range(Npoints): if debug and i == 0: print("angle1, angle2, cen1, cen2") img = ccdimages[i] if ((numpy.sum(img) > cut_off * avg) or (numpy.all(hkls[i] != (0, 0, 0)))): cen1, cen2 = _peak_position(img, nwindow, plot=debug and plot) n1 = numpy.append(n1, cen1) n2 = numpy.append(n2, cen2) ang1 = numpy.append(ang1, angle1[i]) ang2 = numpy.append(ang2, angle2[i]) sang = numpy.append(sang, sampleang[i]) usedhkls.append(hkls[i]) if debug: print(f"{angle1[i]:8.3f} {angle2[i]:8.3f} \t" f"{cen1:.2f} {cen2:.2f}") Nused = len(ang1) usedhkls = numpy.array(usedhkls) if debug: print("Nused / Npoints: %d / %d" % (Nused, Npoints)) # guess initial parameters n10 = numpy.zeros(0, dtype=numpy.double) n20 = n10 ang10 = n10 ang20 = n10 for i in range(Nused): if numpy.all(usedhkls[i] == (0, 0, 0)): n10 = numpy.append(n10, n1[i]) n20 = numpy.append(n20, n2[i]) ang10 = numpy.append(ang10, ang1[i]) ang20 = numpy.append(ang20, ang2[i]) detdir1, detdir2 = _determine_detdir(ang10 - start[3], ang20, n10, n20, detaxis, r_i) epslist = [] paramlist = [] epsmin = numpy.inf fitmin = None print("tiltaz tilt detrot offset sampletilt+azimuth wavelength: " "error (relative) (fittime)") print("------------------------------------------------------------") # find optimal detector rotation (however keep other parameters free) detrot = start[5] if not fix[5]: for detrotstart in numpy.linspace(start[5] - 1, start[5] + 1, 20): start = start[:5] + (detrotstart,) + start[6:] eps, param, fit = _area_detector_calib_fit2( sang, ang1, ang2, n1, n2, usedhkls, experiment, material, detaxis, r_i, detdir1, detdir2, start=start, fix=fix, full_output=True) epslist.append(eps) paramlist.append(param) if epslist[-1] < epsmin: epsmin = epslist[-1] parammin = param fitmin = fit detrot = param[7] if debug: print("single fit") print(param) Ntiltaz = 1 if fix[3] else 5 Ntilt = 1 if fix[4] else 6 Noffset = 1 if fix[6] else 100 if fix[6]: Ntilt = Ntilt * 8 if not fix[4] else Ntilt Ntiltaz = Ntiltaz * 7 if not fix[3] else Ntiltaz startparam = start[:5] + (detrot,) + start[6:] Ntot = Ntiltaz * Ntilt * Noffset ict = 0 for tiltazimuth in numpy.linspace(startparam[3] if fix[3] else 0, 360, Ntiltaz, endpoint=False): for tilt in numpy.linspace(startparam[4] if fix[4] else 0, 4, Ntilt): for offset in numpy.linspace( startparam[6] if fix[6] else - 3 + startparam[6], 3 + startparam[6], Noffset): t1 = time.time() start = (start[:3] + (tiltazimuth, tilt, detrot, offset) + start[7:]) eps, param, fit = _area_detector_calib_fit2( sang, ang1, ang2, n1, n2, usedhkls, experiment, material, detaxis, r_i, detdir1, detdir2, start=start, fix=fix, full_output=True) epslist.append(eps) paramlist.append(param) t2 = time.time() print("%d/%d\t%6.1f %6.2f %8.3f %8.3f %8.3f %7.2f %8.4f: " "%10.4e (%4.2f) (%5.2fsec)" % ((ict, Ntot) + start[3:] + (epslist[-1], epslist[-1] / epsmin, t2 - t1))) ict += 1 if epslist[-1] < epsmin: print("************************") print("new global minimum found") epsmin = epslist[-1] parammin = param fitmin = fit print("new best parameters: %.2f %.2f %10.4e %10.4e %8.4f " "%.1f %.2f %.3f %.3f %.3f %.2f %.4f" % parammin) print("************************\n") (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, stilt, stazimuth, wavelength) = parammin if plot: if fig: plt.figure(fig.number) else: figlabel = "CCD Calib fit %d" i = 1 while figlabel % i in plt.get_figlabels(): i += 1 plt.figure(figlabel % i) nparams = numpy.array(paramlist) neps = numpy.array(epslist) labels = ( 'cch1 (1)', 'cch2 (1)', r'pwidth1 ($\mu$m@1m)', r'pwidth2 ($\mu$m@1m)', 'distance (m)', 'tiltazimuth (deg)', 'tilt (deg)', 'detrot (deg)', 'outerangle offset (deg)', 'sample tilt (deg)', 'st azimuth (deg)', 'wavelength (AA)') xscale = (1., 1., 1.e6, 1.e6, 1., 1., 1., 1., 1., 1., 1., 1.) for p in range(12): ax = plt.subplot(3, 4, p + 1) if plotlog: plt.semilogy(nparams[:, p] * xscale[p], neps, '.k') else: plt.scatter(nparams[:, p] * xscale[p], neps, c=nparams[:, -1], s=10, marker='o', cmap="gnuplot", edgecolor='none') plt.xlabel(labels[p]) if plotlog: plt.semilogy(parammin[p] * xscale[p], epsmin, 'ok', ms=8, mew=2.5, mec='k', mfc='w') else: plt.plot(parammin[p] * xscale[p], epsmin, 'ok', ms=8, mew=2.5, mec='k', mfc='w') plt.ylim(epsmin * 0.7, epsmin * 2.) plt.locator_params(nbins=4, axis='x') if p > 1: if fix[p-2]: ax.set_facecolor('0.85') plt.tight_layout() if config.VERBOSITY >= config.INFO_LOW: print(f"total time needed for fit: {time.time() - t0:.2f}sec") print("fitted parameters: epsilon: %10.4e (%d,%s) " # pylint: disable-next=no-member % (epsmin, fitmin.info, repr(fitmin.stopreason))) print("param: (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, " "tilt, detrot, outerangle_offset, sampletilt, stazimuth, " "wavelength)") print("param: %.2f %.2f %10.4e %10.4e %.4f %.1f %.2f %.3f %.3f %.3f " "%.2f %.4f" % (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, stilt, stazimuth, wavelength)) if config.VERBOSITY > 0: print("please check the resulting data (consider setting plot=True)") print("detector rotation axis / primary beam direction (given by user)" ": %s / %s" % (repr(detaxis), r_i)) print("detector pixel directions / distance: %s %s / %e" % (detdir1, detdir2, distance)) print("\tdetector initialization with: init_area('%s', '%s', " "cch1=%.2f, cch2=%.2f, Nch1=%d, Nch2=%d, pwidth1=%.4e, " "pwidth2=%.4e, distance=%.5f, detrot=%.3f, tiltazimuth=%.1f, " "tilt=%.3f)" % (detdir1, detdir2, cch1, cch2, N1, N2, pwidth1, pwidth2, distance, detrot, tiltazimuth, tilt)) print("AND ALWAYS USE an (additional) OFFSET of %.4fdeg in the " "OUTER DETECTOR ANGLE!" % (outerangle_offset)) return (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, stilt, stazimuth, wavelength), epsmin
def _area_detector_calib_fit2(sang, ang1, ang2, n1, n2, hkls, experiment, material, detaxis, r_i, detdir1, detdir2, start=(None, None, 1, 0, 0, 0, 0, 0, 0, 1.0), fix=(False, False, True, False, False, False, False, False, False, False), full_output=False, debug=False): """ INTERNAL FUNCTION function to calibrate the detector parameters of an area detector it determines the detector tilt possible rotations and offsets in the detector arm angles Parameters ---------- sang : array-like sample rocking angle (rotation direction of inner detector rotation) angle1 : array-like outer detector arm angle angle2 : array-like inner detector arm angle n1, n2 : array-like pixel number at which the beam was observed hkls : tuple or list Miller indices of the reflection were the images were taken (use (0, 0, 0)) for primary beam experiment : Experiment Experiment class object needed to get the UB matrix needed for the hkl peak treatment material : Crystal material used as reference crystal detaxis : str detector arm rotation axis default: ['z+', 'y-'] detdir1, detdir2 : str detector pixel directions of first and second pixel coordinates; e.g. 'y+' r_i : str primary beam direction [xyz][+-] default 'x+' start : tuple or list, optional sequence of start values of the fit for parameters, which can not be estimated automatically. these are: pwidth1, pwidth2, distance, tiltazimuth, tilt, detector_rotation, outerangle_offset, sampletilt, sampletiltazimuth, wavelength. By default: (None, None, 1, 0, 0, 0, 0, 0, 0, 1.0) fix : tuple or list, optional fix parameters of start full_output : bool, optional flag to tell if to return fit object with final parameters and detector directions debug : bool, optional flag to tell if you want to see debug output of the script (switch this to true only if you can handle it :)) Returns ------- float final epsilon of the fit param : list, optional if full_output: fit parameters fit : object, optional if full_output: fit object """ def areapixel2(params, detectorDir1, detectorDir2, r_i, detectorAxis, *args, **kwargs): """ angular to momentum space conversion for pixels of an area detector the center pixel is in direction of self.r_i when detector angles are zero the detector geometry must be given to the routine Parameters ---------- params : list or tuple parameters of the detector calibration model (cch1, cch2, pwidth1, pwidth2, tiltazimuth, tilt, detrot): cch1, 2: center pixel, in direction of r_i at zero detectorAngles; pwidth1, 2: width of one pixel (same unit as distance); distance: distance of center pixel from center of rotation; tiltazimuth: direction of the tilt vector in the detector plane (in degree); tilt: tilt of the detector plane around an axis normal to the direction given by the tiltazimuth; detrot: detector rotation around the primary beam direction as given by r_i detectorDir1 : str direction of the detector (along the pixel direction 1); e.g. 'z+' means higher pixel numbers at larger z positions detectorDir2 : str direction of the detector (along the pixel direction 2); e.g. 'x+' r_i : str primary beam direction e.g. 'x+' detectorAxis : list or tuple detector circles; e.g. ['z+', 'y-'] would mean a detector arm with a two rotations *args : array-like sample, detector angles and channel numbers; *sAngle* sample rocking angle as numpy array, lists or Scalars. *dAngles* as numpy array, lists or Scalars. In total len(detectorAxis) values must be given starting with the outer most circle. All arguments must have the same shape or length. *channel numbers* `n1` and `n2` where the primary beam hits the detector with the same length as the detector values delta : list of array-like, optional giving delta angles to correct the given ones for misalignment. delta must be an numpy array or list of len(*dAngles). used angles are than *args - delta UB : array-like, optional orientation matrix of the sample wl : float or str, optional x-ray wavelength in angstrom deg : bool, optional flag to tell if angles are passed as degree (default: True) Returns ------- ndarray reciprocal space position of detector pixels n1, n2 in a numpy.ndarray of shape (len(args) , 3) """ # check detector circle argument if isinstance(detectorAxis, str): dAxis = list([detectorAxis]) else: dAxis = list(detectorAxis) _sampleAxis_str = dAxis[-1] # add detector rotation around primary beam dAxis += [r_i] _detectorAxis_str = '' for circ in dAxis: _detectorAxis_str += circ Nd = len(dAxis) Nargs = Nd + 2 - 1 + 1 # check detectorDir _area_detdir1 = detectorDir1 _area_detdir2 = detectorDir2 # parse parameter arguments _area_cch1 = float(params[0]) _area_cch2 = float(params[1]) _area_pwidth1 = float(params[2]) _area_pwidth2 = float(params[3]) _area_distance = float(params[4]) _area_tiltazimuth = math.radians(params[5]) _area_tilt = math.radians(params[6]) _area_rot = float(params[7]) _area_ri = xumath.getVector(r_i) * _area_distance # kwargs wl = utilities.wavelength(kwargs.get('wl', 1.)) deg = kwargs.get('deg', True) UB = kwargs.get('UB', numpy.identity(3)) delta = numpy.asarray(kwargs.get('delta', numpy.zeros(Nd)), dtype=numpy.double) if delta.size != Nd - 1 + 1: raise InputError("QConversionPixel: keyword argument delta " "does not have an appropriate shape") # prepare angular arrays from *args # need one sample angle and one detector angle array if len(args) != Nargs: raise InputError("QConversionPixel: wrong amount (%d) of " "arguments given, number of arguments should be " "%d" % (len(args), Nargs)) try: Npoints = len(args[0]) except (TypeError, IndexError): Npoints = 1 sAngles = numpy.array((), dtype=numpy.double) for i in range(1): arg = args[i] if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) arg = arg - delta[i] sAngles = numpy.concatenate((sAngles, arg)) dAngles = numpy.array((), dtype=numpy.double) for i in range(1, Nd): arg = args[i] if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) arg = arg - delta[i] dAngles = numpy.concatenate((dAngles, arg)) # add detector rotation around primary beam dAngles = numpy.concatenate(( dAngles, numpy.ones(arg.shape, dtype=numpy.double) * _area_rot)) # read channel numbers n1 = numpy.array((), dtype=numpy.double) n2 = numpy.array((), dtype=numpy.double) arg = args[Nd] if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) n1 = arg arg = args[Nd + 1] if isinstance(arg, numbers.Number): arg = numpy.array([arg], dtype=numpy.double) elif isinstance(arg, list): arg = numpy.array(arg, dtype=numpy.double) n2 = arg sAngles.shape = (1, Npoints) sAngles = sAngles.transpose() dAngles.shape = (Nd, Npoints) dAngles = dAngles.transpose() if deg: sAngles = radians(sAngles) dAngles = radians(dAngles) qpos = cxrayutilities.ang2q_conversion_area_pixel2( sAngles, dAngles, n1, n2, _area_ri, _sampleAxis_str, _detectorAxis_str, _area_cch1, _area_cch2, _area_pwidth1, _area_pwidth2, _area_detdir1, _area_detdir2, _area_tiltazimuth, _area_tilt, UB, wl, config.NTHREADS) return qpos[:, 0], qpos[:, 1], qpos[:, 2] def afunc(param, x, detectorDir1, detectorDir2, r_i, detectorAxis): """ function for fitting the detector parameters basically this is a wrapper for the areapixel function parameters ---------- param : list or tuple fit parameters (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, sampletiltazimuth, wavelength) x : array-like independent variables; contains (sang, angle1, angle2, n1, n2, hkls) with shape (8, Npoints) detectorDir1 : str direction of the detector (along the pixel direction 1); e.g. 'z+' means higher pixel numbers at larger z positions detectorDir2 : str direction of the detector (along the pixel direction 2); e.g. 'x+' r_i : str primary beam direction e.g. 'x+' detectorAxis : list or tuple detector circles; e.g. ['z+', 'y-'] would mean a detector arm with a two rotations Returns ------- ndarray reciprocal space position of detector pixels n1, n2 in a numpy.ndarray of shape (3, x.shape[1]) """ sang = x[0, :] angle1 = x[1, :] angle2 = x[2, :] n1 = x[3, :] n2 = x[4, :] H = x[5, :] K = x[6, :] L = x[7, :] # use only positive tilt and sample tilt param[6] = abs(param[6]) param[9] = abs(param[9]) wl = param[11] cp = math.cos(math.radians(param[10])) sp = math.sin(math.radians(param[10])) cc = math.cos(math.radians(param[9])) sc = math.sin(math.radians(param[9])) # UB matrix due to tilt at symmetric peak U1 = numpy.array(((cp * cc, -sp, cp * sc), (sp * cc, cp, sp * sc), (-sc, 0, cc)), dtype=numpy.double) U2 = experiment._transform.matrix U = numpy.dot(U1, U2) B = material.B ubmat = numpy.dot(U, B) (qx, qy, qz) = areapixel2( param[:-4], detectorDir1, detectorDir2, r_i, detectorAxis, sang, angle1, angle2, n1, n2, delta=[0, param[8], 0.], distance=1., UB=ubmat, wl=wl) return (qx - H) ** 2 + (qy - K) ** 2 + (qz - L) ** 2 Npoints = len(ang1) # guess initial parameters n10 = numpy.zeros(0, dtype=numpy.double) n20 = n10 ang10 = n10 ang20 = n10 n1s = numpy.zeros(0, dtype=numpy.double) n2s = n1s ang1s = n1s ang2s = n1s sangs = n1s for i in range(Npoints): if numpy.all(hkls[i] == (0, 0, 0)): n10 = numpy.append(n10, n1[i]) n20 = numpy.append(n20, n2[i]) ang10 = numpy.append(ang10, ang1[i]) ang20 = numpy.append(ang20, ang2[i]) else: n1s = numpy.append(n1s, n1[i]) n2s = numpy.append(n2s, n2[i]) ang1s = numpy.append(ang1s, ang1[i]) ang2s = numpy.append(ang2s, ang2[i]) sangs = numpy.append(sangs, sang[i]) # center channel and detector pixel direction and pixel size (s1, i1), r1 = xumath.linregress(ang10 - start[3], n10) (s2, i2), r2 = xumath.linregress(ang10 - start[3], n20) (s3, i3), r3 = xumath.linregress(ang20, n10) (s4, i4), r4 = xumath.linregress(ang20, n20) distance = start[2] if ((r1 + r2 > r3 + r4) and r1 > r2) or ((r1 + r2 < r3 + r4) and r3 < r4): cch1 = i1 cch2 = i4 pwidth1 = 2 * distance / numpy.abs(s1) * math.tan(math.radians(0.5)) pwidth2 = 2 * distance / numpy.abs(s4) * math.tan(math.radians(0.5)) else: cch1 = i3 cch2 = i2 pwidth1 = 2 * distance / numpy.abs(s3) * math.tan(math.radians(0.5)) pwidth2 = 2 * distance / numpy.abs(s2) * math.tan(math.radians(0.5)) if numpy.isscalar(start[0]): pwidth1 = start[0] if numpy.isscalar(start[1]): pwidth2 = start[1] if numpy.isscalar(start[0]) or numpy.isscalar(start[1]): # find biggest correlation and recalculate distance idxmax = numpy.argmax((r1, r2, r3, r4)) s = (s1, s2, s3, s4)[idxmax] distance = abs(s) / math.tan(math.radians(0.5)) / 2 distance *= pwidth1 if idxmax < 2 else pwidth2 tilt = abs(start[4]) tiltazimuth = start[3] detrot = start[5] outerangle_offset = start[6] sampletilt = start[7] stazimuth = start[8] wavelength = start[9] # parameters for the fitting param = (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, stazimuth, wavelength) # determine better start values for sample tilt and azimuth (qx, qy, qz) = areapixel2( param[:-4], detdir1, detdir2, r_i, detaxis, sangs, ang1s, ang2s, n1s, n2s, delta=[0, param[8], 0.], wl=wavelength) if debug: print(f"average qx: {numpy.average(qx):.3f}({numpy.std(qx):.3f})") print(f"average qy: {numpy.average(qy):.3f}({numpy.std(qy):.3f})") print(f"average qz: {numpy.average(qz):.3f}({numpy.std(qz):.3f})") qvecav = (numpy.average(qx), numpy.average(qy), numpy.average(qz)) sampletilt = xumath.VecAngle(experiment.Transform(experiment.ndir), qvecav, deg=True) stazimuth = -xumath.VecAngle( experiment.Transform(experiment.idir), qvecav - xumath.VecDot(experiment.Transform(experiment.ndir), qvecav) * experiment.Transform(experiment.ndir), deg=True) param = (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, stazimuth, wavelength) if debug: print("initial parameters: ") print("primary beam / detector pixel directions / distance: %s / " "%s %s / %e" % (r_i, detdir1, detdir2, distance)) print("param: (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, " "tilt, detrot, outerangle_offset, sampletilt, stazimuth, " "wavelength)") print("param: %.2f %.2f %10.4e %10.4e %.3f %.1f %.2f %.3f %.3f %.3f " "%.2f %.4f" % param) # set data x = numpy.empty((8, Npoints), dtype=numpy.double) x[0, :] = sang x[1, :] = ang1 x[2, :] = ang2 x[3, :] = n1 x[4, :] = n2 x[5, :] = hkls[:, 0] x[6, :] = hkls[:, 1] x[7, :] = hkls[:, 2] data = odr.Data(x, y=1) # define model for fitting model = odr.Model(afunc, extra_args=(detdir1, detdir2, r_i, detaxis), implicit=True) # check if parameters need to be fixed ifixb = () for i in range(len(fix)): ifixb += (int(not fix[i]),) my_odr = odr.ODR(data, model, beta0=param, ifixb=(1, 1) + ifixb, ifixx=(0, 0, 0, 0, 0, 0, 0, 0), stpb=(0.4, 0.4, pwidth1/50., pwidth2/50., distance/1000, 2, 0.125, 0.01, 0.01, 0.01, 1., 0.0001), sclb=(1/abs(cch1), 1/abs(cch2), 1/pwidth1, 1/pwidth2, 1/distance, 1/90., 1/0.2, 1/0.2, 1/0.2, 1/0.1, 1/90., 1.), maxit=200, ndigit=12, sstol=1e-11, partol=1e-11) if debug: my_odr.set_iprint(final=1) my_odr.set_iprint(iter=2) fit = my_odr.run() (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, stazimuth, wavelength) = fit.beta # fix things in parameters tiltazimuth = tiltazimuth % 360. stazimuth = stazimuth % 360. tilt = abs(tilt) sampletilt = abs(sampletilt) final_q = afunc([cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, stazimuth, wavelength], x, detdir1, detdir2, r_i, detaxis) final_error = numpy.mean(final_q) if debug: # pylint: disable-next=no-member print("fitted parameters: (%e, %d, %s) " % (final_error, fit.info, repr(fit.stopreason))) print("primary beam / detector pixel directions / distance: %s / %s " "%s / %e" % (r_i, detdir1, detdir2, distance)) print("param: (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, " "tilt, detrot, outerangle_offset, sampletilt, sampletiltazimuth," " wavelength)") print("param: %.2f %.2f %10.4e %10.4e %.3f %.1f %.2f %.3f %.3f %.3f " "%.2f %.4f" % (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, stazimuth, wavelength)) if full_output: return final_error, (cch1, cch2, pwidth1, pwidth2, distance, tiltazimuth, tilt, detrot, outerangle_offset, sampletilt, stazimuth, wavelength), fit return final_error #################################################
[docs] def psd_refl_align(primarybeam, angles, channels, plot=True): """ function which calculates the angle at which the sample is parallel to the beam from various angles and detector channels from the reflected beam. The function can be used during the half beam alignment with a linear detector. Parameters ---------- primarybeam : int primary beam channel number angles : list or array-like incidence angles channels : list or array-like corresponding detector channels plot : bool, optional flag to specify if a visualization of the fit is wanted default : True Returns ------- float angle at which the sample is parallel to the beam Examples -------- >>> zeroangle = psd_refl_align(500, [0, 0.1, 0.2, 0.3], ... [550, 600, 640, 700]) XU.analysis.psd_refl_align: sample is parallel to beam at goniometer angle\ -0.0986 (R^2=0.9942) """ if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis.psd_refl_' 'align') p, rsq = xumath.linregress(numpy.asarray(channels), numpy.asarray(angles)) zeropos = numpy.polyval(p, primarybeam) if plot: xmin = min(min(channels), primarybeam) xmax = max(max(channels), primarybeam) ymin = min(min(angles), zeropos) ymax = max(max(angles), zeropos) # open new figure for the plot plt.figure() plt.plot(channels, angles, 'xk', ms=8., mew=2.) plt.plot([xmin - (xmax - xmin) * 0.1, xmax + (xmax - xmin) * 0.1], numpy.polyval(p, [xmin - (xmax - xmin) * 0.1, xmax + (xmax - xmin) * 0.1]), '-g', linewidth=1.5) ax = plt.gca() plt.grid() ax.set_xlim(xmin - (xmax - xmin) * 0.15, xmax + (xmax - xmin) * 0.15) ax.set_ylim(ymin - (ymax - ymin) * 0.15, ymax + (ymax - ymin) * 0.15) plt.vlines(primarybeam, ymin - (ymax - ymin) * 0.1, ymax + (ymax - ymin) * 0.1, linewidth=1.5) plt.xlabel("PSD Channel") plt.ylabel("sample angle") plt.tight_layout() if config.VERBOSITY >= config.INFO_LOW: print("XU.analysis.psd_refl_align: sample is parallel to beam at " "goniometer angle %7.4f (R^2=%6.4f)" % (zeropos, rsq)) return zeropos
################################################# # miscut calculation from alignment in 2 and # more azimuths #################################################
[docs] def miscut_calc(phi, aomega, zeros=None, omega0=None, plot=True): """ function to calculate the miscut direction and miscut angle of a sample by fitting a sinusoidal function to the variation of the aligned omega values of more than two reflections. The function can also be used to fit reflectivity alignment values in various azimuths. Parameters ---------- phi : list, tuple or array-like azimuths in which the reflection was aligned (deg) aomega : list, tuple or array-like aligned omega values (deg) zeros : list, tuple or array-like, optional angles at which surface is parallel to the beam (deg). For the analysis the angles (aomega - zeros) are used. omega0 : float, optional if specified the nominal value of the reflection is not included as fit parameter, but is fixed to the specified value. This value is MANDATORY if ONLY TWO AZIMUTHs are given. plot : bool, optional flag to specify if a visualization of the fit is wanted. default: True Returns ------- omega0 : float the omega value of the reflection should be close to the nominal one phi0 : float the azimuth in which the primary beam looks upstairs miscut : float amplitude of the sinusoidal variation == miscut angle """ if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis.miscut_' 'calc') if zeros is not None: om = (numpy.array(aomega) - numpy.array(zeros)) else: om = numpy.array(aomega) a = numpy.array(phi) if omega0 is None: # first guess for the parameters # omega0, phi0, miscut p0 = (om.mean(), a[om.argmax()], om.max() - om.min()) def fitfunc(p, phi): return abs(p[2]) * \ cos(radians(phi - (p[1] % 360.))) + p[0] else: # first guess for the parameters p0 = (a[om.argmax()], om.max() - om.min()) # omega0, phi0, miscut def fitfunc(p, phi): return abs(p[1]) * \ cos(radians(phi - (p[0] % 360.))) + omega0 def errfunc(p, phi, om): return fitfunc(p, phi) - om p1, success = optimize.leastsq(errfunc, p0, args=(a, om), maxfev=10000) if config.VERBOSITY >= config.INFO_ALL: print("xu.analysis.miscut_calc: leastsq optimization return value: " "%d" % success) if plot: plt.figure() plt.plot(a, om, 'xk', mew=2, ms=8) linx = numpy.linspace(a.min() - 45, a.min() + 360 - 45, num=1000) plt.plot(linx, fitfunc(p1, linx), '-g', linewidth=1.5) plt.grid() plt.xlabel("azimuth") plt.ylabel("aligned sample angle") if omega0 is None: ret = [p1[0], p1[1] % 360., abs(p1[2])] else: ret = [omega0] + [p1[0] % 360., abs(p1[1])] if config.VERBOSITY >= config.INFO_LOW: print("xu.analysis.miscut_calc: \n" "\t fitted reflection angle: %8.3f \n" "\t looking upstairs at phi: %8.2f \n" "\t miscut angle: %8.3f \n" % (ret[0], ret[1], ret[2])) return ret
################################################# # correct substrate Bragg peak position in # reciprocal space maps #################################################
[docs] def fit_bragg_peak(om, tt, psd, omalign, ttalign, exphxrd, frange=(0.03, 0.03), peaktype='Gauss', plot=True): r""" helper function to determine the Bragg peak position in a reciprocal space map used to obtain the position needed for correction of the data. the determination is done by fitting a two dimensional Gaussian (xrayutilities.math.Gauss2d) or Lorentzian (xrayutilities.math.Lorentz2d) PLEASE ALWAYS CHECK THE RESULT CAREFULLY! Parameters ---------- om, tt : array-like angular coordinates of the measurement either with size of psd or of psd.shape[0] psd : array-like intensity values needed for fitting omalign : float aligned omega value, used as first guess in the fit ttalign : float aligned two theta values used as first guess in the fit these values are also used to set the range for the fit: the peak should be within +/-frange\AA^{-1} of those values exphxrd : Experiment experiment class used for the conversion between angular and reciprocal space. frange : tuple of float, optional data range used for the fit in both directions (see above for details default:(0.03, 0.03) unit: \AA^{-1}) peaktype : {'Gauss', 'Lorentz'} peak type to fit plot : bool, optional if True (default) function will plot the result of the fit in comparison with the measurement. Returns ------- omfit, ttfit : float fitted angular values params : list fit parameters (of the Gaussian/Lorentzian) covariance : ndarray covariance matrix of the fit parameters """ if plot: plot, plt = utilities.import_matplotlib_pyplot('XU.analysis.fit_bragg' '_peak') if peaktype == 'Gauss': func = xumath.Gauss2d elif peaktype == 'Lorentz': func = xumath.Lorentz2d else: raise InputError("peaktype must be either 'Gauss' or 'Lorentz'") if om.size != psd.size: [_, qy, qz] = exphxrd.Ang2Q.linear(om, tt) else: [_, qy, qz] = exphxrd.Ang2Q(om, tt) [_, qysub, qzsub] = exphxrd.Ang2Q(omalign, ttalign) params = [qysub, qzsub, 0.001, 0.001, psd.max(), 0, 0.] drange = [qysub - frange[0], qysub + frange[0], qzsub - frange[1], qzsub + frange[1]] params, covariance = xumath.fit_peak2d( qy.flatten(), qz.flatten(), psd.flatten(), params, drange, func, maxfev=10000) # correct params params[6] = params[6] % (numpy.pi) if params[5] < 0: params[5] = 0 [omfit, _, _, ttfit] = exphxrd.Q2Ang((0, params[0], params[1]), trans=False, geometry="real") if config.VERBOSITY >= config.INFO_LOW: print("XU.analysis.fit_bragg_peak:fitted peak angles: \n\tom =%8.4f\n" "\ttt =%8.4f" % (omfit, ttfit)) if plot: plt.figure() plt.clf() from ..gridder2d import Gridder2D gridder = Gridder2D(50, 50) mask = (qy.flatten() > drange[0]) * (qy.flatten() < drange[1]) * \ (qz.flatten() > drange[2]) * (qz.flatten() < drange[3]) gridder(qy.flatten()[mask], qz.flatten()[mask], psd.flatten()[mask]) # calculate intensity which should be plotted INT = utilities.maplog(gridder.data.transpose(), 4, 0) QXm = gridder.xmatrix QZm = gridder.ymatrix cl = plt.contour(gridder.xaxis, gridder.yaxis, utilities.maplog(func(QXm, QZm, *params), 4, 0).T, 8, colors='k', linestyles='solid') cf = plt.contourf(gridder.xaxis, gridder.yaxis, INT, 35) cf.collections[0].set_label('data') cl.collections[0].set_label('fit') # plt.legend() # for some reason not working? plt.colorbar(extend='min') plt.title("plot shows only coarse data! fit used raw data!") return omfit, ttfit, params, covariance