Source code for xconf.spectrometer

#!/usr/bin/python
# -*- coding: utf-8 -*-

"""
Handle **von Hamos** type and **Johann** type X-ray spectrometers.
"""

__copyright__ = "Copyright (C) 2019 Florian Otte"

__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""

import math

import xconf.misc.units as units
from xconf.crystal import diffraction_allowed


[docs]class Spectrometer(object): """Spectrometer base class. This spectrometer consists of *analyzers*.""" def __init__(self, source, analyzers, detector): self.source = source self.analyzers = analyzers self.detector = detector
[docs] def compute_configurations(self, energy, max_order, min_angle, detector_pos, vhs_tilt, filter = []): """ Return all feasible spectrometer configuration, considering each analyzer. This might be used to test different analyzer crystals, or to evaluate the same analyzer crystal at different positions. Giving a list of analyzer, i.e. [['Ge', 1, 1, 0],...] as *filter* will then only return configurations with those analyzers. """ configurations = { 'id':[], 'analyzer':[], 'energy':[], 'brangle':[], 'vhs_tilt':[], 'hkl':[], 'order':[], 'c_distance':[], 'c_horizontal':[], 'c_vertical':[], 'd_distance':[], 'd_horizontal':[], 'd_vertical':[], 'xtal_e_min':[], 'xtal_e_max':[], 'total_e_min':[], 'total_e_max':[], 'total_e_range':[], 'geometric_bandwidth':[], 'xtal_bandwidth':[], 'total_bandwidth':[], 'det_e_min':[], 'det_e_max':[], 'solid_angle':[] } for analyzer in self.analyzers: allowed = diffraction_allowed(analyzer.hkl, max_order) for hkl in allowed: h,k,l = hkl label = list([analyzer.material_string, h, k, l]) if len(filter) > 0 and not label in filter: continue angle = analyzer.get_bragg_angle(energy, hkl) if angle is not None and angle >= min_angle: geom = self.calculate_geometry( energy, hkl, self.source, analyzer, self.detector, vhs_tilt, detector_pos) # Spectrometer energy range xtal_e_min = analyzer.get_energy(geom['xtal_th_max'], hkl) xtal_e_max = analyzer.get_energy(geom['xtal_th_min'], hkl) det_e_min =analyzer.get_energy(geom['det_th_max'], hkl) det_e_max =analyzer.get_energy(geom['det_th_min'], hkl) total_e_min = max(det_e_min, xtal_e_min) total_e_max = min(det_e_max, xtal_e_max) total_range = total_e_max - total_e_min # Spectrometer resolution try: # If there is no compatible arrangement of analyzer and # detector, no geometric resolution contribution can be # calculated geom_bw = self.get_geometric_resolution( self.source, analyzer, self.detector, geom, hkl) except: geom_bw = [units.Energy(0.0, 'eV')] xtal_bw = analyzer.get_spectral_bandwidth( energy, hkl, angle) total_bw = xtal_bw + max(geom_bw) # Covered solid angle solid_angle = analyzer.get_solid_angle(angle) # Build dict label = '{}({},{},{})'.format( analyzer.material_string, *list(analyzer.hkl)) d = { 'id':None, 'analyzer':label, 'energy':energy.convert('eV'), 'brangle':angle.convert('deg'), 'vhs_tilt':geom['vhs_tilt'], 'hkl':hkl, 'order':int(hkl[0] / analyzer.hkl[0]), 'c_distance':geom['analyzer_dis'], 'c_horizontal':geom['analyzer_horizontal'], 'c_vertical':geom['analyzer_vertical'], 'd_distance':geom['detector_pos'], 'd_horizontal':geom['detector_horizontal'], 'd_vertical':geom['detector_vertical'], 'xtal_e_min':xtal_e_min.convert('eV'), 'xtal_e_max':xtal_e_max.convert('eV'), 'det_e_min':det_e_min.convert('eV'), 'det_e_max':det_e_max.convert('eV'), 'total_e_min':total_e_min.convert('eV'), 'total_e_max':total_e_max.convert('eV'), 'total_e_range':total_range.convert('eV'), 'geometric_bandwidth':max(geom_bw).convert('eV'), 'xtal_bandwidth':xtal_bw.convert('eV'), 'total_bandwidth':total_bw.convert('eV'), 'solid_angle':solid_angle.convert('msr') } for key in d: configurations[key].append(d[key]) return configurations
[docs]class vonHamosSpectrometer(Spectrometer): """Describe the **von Hamos** type spectrometer.""" def __init__(self, *args, **kwargs): Spectrometer.__init__(self, *args, **kwargs)
[docs] def calculate_geometry(self, energy, hkl, source, analyzer, detector, vhs_tilt, detector_pos = None): """ Calculate the spectrometer geometry, i.e. the distances and angles between *source*, *analyzer* and *detector* for a given *energy*. The detector position *detector_pos* can be forced to a fixed position. The return value is a dict (see below): ============ ======================================================= key value ============ ======================================================= detector_pos Distance from sample to detector. analyzer_dis Distance from sample to the analyzer crystal. analyzer_pos Offset of the analyzer crystal from sample along a line between sample and detector. xtal_th_min Minimal angle on xtal surface xtal_th_max Maximal angle on xtal surface det_th_min Minimal angle on detector surface det_th_max Maximal angle on detector surface ============ ======================================================= """ bragg_angle = analyzer.get_bragg_angle(energy, hkl) if detector_pos is not None: units.check_parameter("Detector position", detector_pos) roc = analyzer.roc d_crystal = roc / math.sin(bragg_angle.convert('rad').free()) h_crystal = roc / math.tan(bragg_angle.convert('rad').free()) d_detector = 2* h_crystal if not detector_pos else detector_pos h_analyzer_size = 0.5*analyzer.height h_detector_size = 0.5*detector.get_chip_size() ang = vhs_tilt.convert('rad').free() + bragg_angle.convert('rad').free() analyzer_vertical = math.cos(ang) * d_crystal analyzer_horizontal = math.sin(ang) * d_crystal detector_vertical = math.cos(vhs_tilt.convert('rad').free()) * d_detector detector_horizontal = math.sin(vhs_tilt.convert('rad').free()) * d_detector xtal_th_min = math.atan((roc/((h_crystal + h_analyzer_size))).free()) xtal_th_max = math.atan((roc/((h_crystal - h_analyzer_size))).free()) det_th_min = math.atan((2*roc/(d_detector + h_detector_size)).free()) det_th_max = math.atan((2*roc/(d_detector - h_detector_size)).free()) geometry = dict( detector_pos = d_detector, detector_horizontal = detector_horizontal, detector_vertical = detector_vertical, analyzer_dis = d_crystal, analyzer_pos = h_crystal, analyzer_horizontal = analyzer_horizontal, analyzer_vertical = analyzer_vertical, xtal_th_min = units.Angle(xtal_th_min, "rad"), xtal_th_max = units.Angle(xtal_th_max, "rad"), det_th_min = units.Angle(det_th_min, "rad"), det_th_max = units.Angle(det_th_max, "rad"), vhs_tilt = vhs_tilt ) return geometry
[docs] def get_geometric_resolution(self, source, analyzer, detector, geometry, hkl): """ Compute the minimal and maximal angle under which X-ray photons are allowed to bragg-scatter from the analyzer surface to still reach a certain pixel on the detector. Those angles determine the geometric energy bandwidth of the pixel. Specifiy a *source*, *analyzer* and *detector*, set up in a *geometry*. The Miller indizes *hkl* are needed to determine the energy for a given angle on the analyzer. """ pos_a = geometry['analyzer_pos'] pos_d = geometry['detector_pos'] a1 = pos_a - 0.5*analyzer.height a2 = pos_a + 0.5*analyzer.height d1 = pos_d - 0.5*detector.get_chip_size() d2 = pos_d + 0.5*detector.get_chip_size() half_bs = 0.5*source.beam_size roc = analyzer.roc t1 = (a1 + half_bs) / roc * 2*roc - half_bs t2 = (a2 - half_bs) /roc * 2*roc + half_bs w1 = max(t1, d1) w2 = min(t2, d2) if w1 >= w2: raise Exception("Invalid geometry.") resolutions = [.0] * detector.pixels for p in range(detector.pixels): dpix1 = max(d1 + p*detector.ps, w1) dpix2 = min(d1 + (p+1)*detector.ps, w2) if dpix1 >= dpix2: resolutions[p] = units.Energy(0.0, 'eV') else: theta_w1 = math.atan(((2 *roc) / (dpix1 - half_bs)).free()) theta_w2 = math.atan(((2 *roc) / (dpix2 + half_bs)).free()) e1 = analyzer.get_energy(units.Angle(theta_w1, 'rad'), hkl) e2 = analyzer.get_energy(units.Angle(theta_w2, 'rad'), hkl) resolutions[p] = e2-e1 return resolutions
[docs]class JohannSpectrometer(Spectrometer): """Describe the **Johann** type spectrometer.""" def __init__(self, *args, **kwargs): Spectrometer.__init__(self, *args, **kwargs) raise NotImplementedError('Johann geometry has not yet been added.')