diff --git a/.gitignore b/.gitignore index fa997c4..5205740 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ *~ *.bak -__py_cache__ +*.log +__pycache__ raw - +.idea +test_data diff --git a/libeos/__init__.py b/libeos/__init__.py new file mode 100644 index 0000000..1099066 --- /dev/null +++ b/libeos/__init__.py @@ -0,0 +1,6 @@ +""" +Package to handle data redction at AMOR instrument to be used by eos.py script. +""" + +__version__ = '2.0' +__date__ = '2024-03-04' diff --git a/libeos/command_line.py b/libeos/command_line.py new file mode 100644 index 0000000..7fae9e5 --- /dev/null +++ b/libeos/command_line.py @@ -0,0 +1,193 @@ +import argparse +from datetime import datetime + +from .logconfig import update_loglevel +from .options import ReaderConfig, EOSConfig, ExperimentConfig, OutputConfig, ReductionConfig + + +def commandLineArgs(): + """ + Process command line argument. + The type of the default values is used for conversion and validation. + """ + msg = "eos reads data from (one or several) raw file(s) of the .hdf format, \ + performs various corrections, conversations and projections and exports\ + the resulting reflectivity in an orso-compatible format." + clas = argparse.ArgumentParser(description = msg) + + input_data = clas.add_argument_group('input data') + input_data.add_argument("-n", "--fileIdentifier", + default = ['0'], + nargs = '+', + help = "file number(s) or offset (if negative)") + input_data.add_argument("-r", "--normalisationFileIdentifier", + default = [], + nargs = '+', + help = "file number(s) of normalisation measurement") + input_data.add_argument("-d", "--dataPath", + type = str, + default = '.', + help = "relative path to directory with .hdf files") + input_data.add_argument("-Y", "--year", + default = datetime.year, + type = int, + help = "year the measurement was performed") + input_data.add_argument("-sub", "--subtract", + help = "R(q_z) curve to be subtracted (in .Rqz.ort format)") + + output = clas.add_argument_group('output') + output.add_argument("-o", "--outputName", + default = "fromEOS", + help = "output file name (withot suffix)") + output.add_argument("-of", "--outputFormat", + nargs = '+', + default = ['Rqz.ort']) + output.add_argument("--offSpecular", + type = bool, + default = False, + ) + output.add_argument("-a", "--qResolution", + default = 0.01, + type = float, + help = "q_z resolution") + output.add_argument("-ts", "--timeSlize", + nargs = '+', + type = float, + help = "time slizing ,[ [,stop]]") + output.add_argument("-s", "--scale", + nargs = '+', + default = [1], + type = float, + help = "scaling factor for R(q_z)") + output.add_argument("-S", "--autoscale", + nargs = 2, + type = float, + help = "scale to 1 in the given q_z range") + + masks = clas.add_argument_group('masks') + masks.add_argument("-l", "--lambdaRange", + default = [2., 15.], + nargs = 2, + type = float, + help = "wavelength range") + masks.add_argument("-t", "--thetaRange", + default = [-12., 12.], + nargs = 2, + type = float, + help = "absolute theta range") + masks.add_argument("-T", "--thetaRangeR", + default = [-12., 12.], + nargs = 2, + type = float, + help = "relative theta range") + masks.add_argument("-y", "--yRange", + default = [11, 41], + nargs = 2, + type = int, + help = "detector y range") + masks.add_argument("-q", "--qzRange", + default = [0.005, 0.30], + nargs = 2, + type = float, + help = "q_z range") + + overwrite = clas.add_argument_group('overwrite') + overwrite.add_argument("-cs", "--chopperSpeed", + type = float, + help = "chopper speed in rpm") + overwrite.add_argument("-cp", "--chopperPhase", + default = -13.5, + type = float, + help = "chopper phase") + overwrite.add_argument("-co", "--chopperPhaseOffset", + default = -5, + type = float, + help = "phase offset between chopper opening and trigger pulse") + overwrite.add_argument("-m", "--muOffset", + default = 0., + type = float, + help = "mu offset") + overwrite.add_argument("-mu", "--mu", + default = 0, + type = float, + help ="value of mu") + overwrite.add_argument("-nu", "--nu", + default = 0, + type = float, + help = "value of nu") + overwrite.add_argument("-sm", "--sampleModel", + type = str, + help = "1-line orso sample model description") + + misc = clas.add_argument_group('misc') + misc.add_argument('-v', '--verbose', action='store_true') + misc.add_argument('-vv', '--debug', action='store_true') + + return clas.parse_args() + + +def expand_file_list(short_notation): + """Evaluate string entry for file number lists""" + #log().debug('Executing get_flist') + file_list=[] + for i in short_notation.split(','): + if '-' in i: + if ':' in i: + step = i.split(':', 1)[1] + file_list += range(int(i.split('-', 1)[0]), int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, int(step)) + else: + step = 1 + file_list += range(int(i.split('-', 1)[0]), int(i.split('-', 1)[1])+1, int(step)) + else: + file_list += [int(i)] + + return sorted(file_list) + + +def output_format_list(outputFormat): + format_list = [] + if 'ort' in outputFormat or 'Rqz.ort' in outputFormat or 'Rqz' in outputFormat: + format_list.append('Rqz.ort') + if 'ort' in outputFormat or 'Rlt.ort' in outputFormat or 'Rlt' in outputFormat: + format_list.append('Rlt.ort') + if 'orb' in outputFormat or 'Rqz.orb' in outputFormat or 'Rqz' in outputFormat: + format_list.append('Rqz.orb') + if 'orb' in outputFormat or 'Rlt.orb' in outputFormat or 'Rlt' in outputFormat: + format_list.append('Rlt.orb') + return sorted(format_list, reverse=True) + +def command_line_options(): + clas = commandLineArgs() + update_loglevel(clas.verbose, clas.debug) + + reader_config = ReaderConfig( + year=clas.year, + dataPath=clas.dataPath) + experiment_config = ExperimentConfig( + sampleModel=clas.sampleModel, + chopperPhase=clas.chopperPhase, + chopperPhaseOffset=clas.chopperPhaseOffset, + yRange=clas.yRange, + lambdaRange=clas.lambdaRange, + qzRange=clas.qzRange, + offSpecular=clas.offSpecular, + mu=clas.mu, + nu=clas.nu, + muOffset=clas.muOffset + ) + reduction_config = ReductionConfig( + qResolution=clas.qResolution, + autoscale=clas.autoscale, + thetaRange=clas.thetaRange, + thetaRangeR=clas.thetaRangeR, + fileIdentifier=clas.fileIdentifier, + scale=clas.scale, + subtract=clas.subtract, + normalisationFileIdentifier=clas.normalisationFileIdentifier, + timeSlize=clas.timeSlize + ) + output_config = OutputConfig( + outputFormats=output_format_list(clas.outputFormat), + outputName=clas.outputName + ) + return EOSConfig(reader_config, experiment_config, reduction_config, output_config) \ No newline at end of file diff --git a/libeos/const.py b/libeos/const.py new file mode 100644 index 0000000..27cd24a --- /dev/null +++ b/libeos/const.py @@ -0,0 +1,6 @@ +""" +Constants used in data reduction. +""" + +hdm = 6.626176e-34/1.674928e-27 # h / m +lamdaCut = 2.5 # Aa diff --git a/libeos/file_reader.py b/libeos/file_reader.py new file mode 100644 index 0000000..9824c9d --- /dev/null +++ b/libeos/file_reader.py @@ -0,0 +1,351 @@ +import logging +import os +import subprocess +import sys +from datetime import datetime +from typing import List + +import h5py +import numpy as np +from orsopy import fileio + +from . import const +from .header import Header +from .instrument import Detector +from .options import ExperimentConfig, ReaderConfig + +try: + from . import nb_helpers +except Exception: + nb_helpers = None + + +class AmorData: + """read meta-data and event streams from .hdf file(s), apply filters and conversions""" + chopperDetectorDistance: float + chopperDistance: float + chopperPhase: float + chopperSpeed: float + ctime: float + div: float + data_file_numbers: List[int] + delta_z: np.ndarray + detZ_e: np.ndarray + lamda_e: np.ndarray + wallTime_e: np.ndarray + kad: float + kap: float + lambdaMax: float + lambda_e: np.ndarray + monitor1: float + monitor2: float + mu: float + nu: float + tau: float + tofCut: float + start_date: str + + #------------------------------------------------------------------------------------------------- + def __init__(self, header: Header, reader_config: ReaderConfig, config: ExperimentConfig, + short_notation:str, norm=False): + self.startTime = reader_config.startTime + self.header = header + self.config = config + self.reader_config = reader_config + self.expand_file_list(short_notation) + self.read_data(norm=norm) + + #------------------------------------------------------------------------------------------------- + def read_data(self, norm=False): + self.file_list = [] + for number in self.data_file_numbers: + self.file_list.append(self.path_generator(number)) + ## read specific meta data and measurement from first file + if norm: + self.readHeaderInfo = False + else: + self.readHeaderInfo = True + + _detZ_e = [] + _lamda_e = [] + _wallTime_e = [] + for file in self.file_list: + self.read_individual_data(file, norm) + _detZ_e = np.append(_detZ_e, self.detZ_e) + _lamda_e = np.append(_lamda_e, self.lamda_e) + _wallTime_e = np.append(_wallTime_e, self.wallTime_e) + self.detZ_e = _detZ_e + self.lamda_e = _lamda_e + self.wallTime_e = _wallTime_e + + #------------------------------------------------------------------------------------------------- + def path_generator(self, number): + fileName = f'amor{self.reader_config.year}n{number:06d}.hdf' + if os.path.exists(os.path.join(self.reader_config.dataPath,fileName)): + path = self.reader_config.dataPath + elif os.path.exists(fileName): + path = '.' + elif os.path.exists(os.path.join('.','raw', fileName)): + path = os.path.join('.','raw') + elif os.path.exists(os.path.join('..','raw', fileName)): + path = os.path.join('..','raw') + elif os.path.exists(f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}/{fileName}'): + path = f'/afs/psi.ch/project/sinqdata/{self.reader_config.year}/amor/{int(number/1000)}' + else: + sys.exit(f'# ERROR: the file {fileName} is nowhere to be found!') + return os.path.join(path, fileName) + #------------------------------------------------------------------------------------------------- + def expand_file_list(self, short_notation): + """Evaluate string entry for file number lists""" + #log().debug('Executing get_flist') + file_list=[] + for i in short_notation.split(','): + if '-' in i: + if ':' in i: + step = i.split(':', 1)[1] + file_list += range(int(i.split('-', 1)[0]), int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, int(step)) + else: + step = 1 + file_list += range(int(i.split('-', 1)[0]), int(i.split('-', 1)[1])+1, int(step)) + else: + file_list += [int(i)] + self.data_file_numbers=sorted(file_list) + #------------------------------------------------------------------------------------------------- + def resolve_pixels(self): + """determine spatial coordinats and angles from pixel number""" + nPixel = Detector.nWires * Detector.nStripes * Detector.nBlades + pixelID = np.arange(nPixel) + (bladeNr, bPixel) = np.divmod(pixelID, Detector.nWires * Detector.nStripes) + (bZi, detYi) = np.divmod(bPixel, Detector.nStripes) # z index on blade, y index on detector + detZi = bladeNr * Detector.nWires + bZi # z index on detector + detX = bZi * Detector.dX # x position in detector + # detZ = Detector.zero - bladeNr * Detector.bladeZ - bZi * Detector.dZ # z position on detector + bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*Detector.bladeZ / Detector.distance) ) + delta = (Detector.nBlades/2. - bladeNr) * bladeAngle \ + - np.rad2deg( np.arctan(bZi*Detector.dZ / ( Detector.distance + bZi * Detector.dX) ) ) + self.delta_z = delta[detYi==1] + return np.vstack((detYi.T, detZi.T, detX.T, delta.T)).T + #return matr + #------------------------------------------------------------------------------------------------- + def read_individual_data(self, fileName, norm=False): + self.hdf = h5py.File(fileName, 'r', swmr=True) + + if self.readHeaderInfo: + self.read_header_info() + + logging.info(f'# data from file: {fileName}') + self.read_individual_header() + + # add header content + if self.readHeaderInfo: + self.readHeaderInfo = False + self.header.measurement_instrument_settings = fileio.InstrumentSettings( + incident_angle = fileio.ValueRange(round(self.mu+self.kap+self.kad-0.5*self.div, 3), + round(self.mu+self.kap+self.kad+0.5*self.div, 3), + 'deg'), + wavelength = fileio.ValueRange(const.lamdaCut, self.config.lambdaRange[1], 'angstrom'), + polarization = fileio.Polarization.unpolarized, + ) + self.header.measurement_instrument_settings.mu = fileio.Value(round(self.mu, 3), 'deg', comment='sample angle to horizon') + self.header.measurement_instrument_settings.nu = fileio.Value(round(self.nu, 3), 'deg', comment='detector angle to horizon') + if norm: + self.header.measurement_additional_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=self.fileDate)) + else: + self.header.measurement_data_files.append(fileio.File(file=fileName.split('/')[-1], timestamp=self.fileDate)) + logging.info(f'# mu = {self.mu:6.3f}, nu = {self.nu:6.3f}, kap = {self.kap:6.3f}, kad = {self.kap:6.3f}') + + # TODO: should extract monitor from counts or beam current times time + self.monitor1 = self.ctime + self.monitor2 = self.monitor1 + + self.read_event_stream() + totalNumber = np.shape(self.tof_e)[0] + + self.extract_walltime(norm) + + if True: + self.filter_strange_times() + self.merge_frames() + + self.filter_project_x() + + # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau + # TODO: check for correctness + if not self.config.offSpecular: + self.tof_e -= ( self.delta_e / 180. ) * self.tau + + self.calculate_derived_properties() + + self.filter_qz_range(norm) + + logging.info(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') + + def filter_qz_range(self, norm): + if self.config.qzRange[1]<0.3 and not norm: + self.mask_e = np.logical_and(self.mask_e, + (self.config.qzRange[0]<=self.qz_e) & (self.qz_e<=self.config.qzRange[1])) + self.detZ_e = self.detZ_e[self.mask_e] + self.lamda_e = self.lamda_e[self.mask_e] + self.wallTime_e = self.wallTime_e[self.mask_e] + + def calculate_derived_properties(self): + self.lamdaMax = const.lamdaCut+1.e13*self.tau*const.hdm/(self.chopperDetectorDistance+124.) + if nb_helpers and not self.config.offSpecular: + self.lamda_e, self.qz_e, self.mask_e = nb_helpers.calculate_derived_properties_focussing( + self.tof_e, self.detXdist_e, self.delta_e, self.mask_e, + self.config.lambdaRange[0], self.config.lambdaRange[1], self.nu, self.mu, + self.chopperDetectorDistance, const.hdm + ) + return + # lambda + self.lamda_e = (1.e13*const.hdm)*self.tof_e/(self.chopperDetectorDistance+self.detXdist_e) + self.mask_e = np.logical_and(self.mask_e, (self.config.lambdaRange[0]<=self.lamda_e) & ( + self.lamda_e<=self.config.lambdaRange[1])) + # alpha_f + alphaF_e = self.nu-self.mu+self.delta_e + # q_z + if self.config.offSpecular: + alphaI = self.kap+self.kad+self.mu + self.qz_e = 2*np.pi*((np.sin(np.deg2rad(alphaF_e))+np.sin(np.deg2rad(alphaI)))/self.lamda_e) + self.qx_e = 2*np.pi*((np.cos(np.deg2rad(alphaF_e))-np.cos(np.deg2rad(alphaI)))/self.lamda_e) + self.header.measurement_scheme = 'energy-dispersive', + else: + self.qz_e = 4*np.pi*(np.sin(np.deg2rad(alphaF_e))/self.lamda_e) + # qx_e = 0. + self.header.measurement_scheme = 'angle- and energy-dispersive' + + def filter_project_x(self): + pixelLookUp = self.resolve_pixels() + if nb_helpers: + (self.detZ_e, self.detXdist_e, self.delta_e, self.mask_e) = nb_helpers.filter_project_x( + pixelLookUp, self.pixelID_e.astype(np.int64), self.config.yRange[0], self.config.yRange[1] + ) + else: + # resolve pixel ID into y and z indicees, x position and angle + (detY_e, self.detZ_e, self.detXdist_e, self.delta_e) = pixelLookUp[np.int_(self.pixelID_e)-1, :].T + # define mask and filter y range + self.mask_e = (self.config.yRange[0]<=detY_e) & (detY_e<=self.config.yRange[1]) + + def merge_frames(self): + total_offset = self.tofCut+self.tau*self.config.chopperPhaseOffset/180. + if nb_helpers: + self.tof_e = nb_helpers.merge_frames(self.tof_e, self.tofCut, self.tau, total_offset) + else: + self.tof_e = np.remainder(self.tof_e-(self.tofCut-self.tau), self.tau)+total_offset # tof shifted to 1 frame + + def filter_strange_times(self): + # filter 'strange' tof times > 2 tau + filter_e = (self.tof_e<=2*self.tau) + self.tof_e = self.tof_e[filter_e] + self.pixelID_e = self.pixelID_e[filter_e] + self.wallTime_e = self.wallTime_e[filter_e] + if np.shape(filter_e)[0]-np.shape(self.tof_e)[0]>0.5: + logging.warning(f'# strange times: {np.shape(filter_e)[0]-np.shape(self.tof_e)[0]}') + + def extract_walltime(self, norm): + if nb_helpers: + self.wallTime_e = nb_helpers.extract_walltime(self.tof_e, self.dataPacket_p, self.dataPacketTime_p) + else: + totalNumber = np.shape(self.tof_e)[0] + self.wallTime_e = np.empty(totalNumber) + for i in range(len(self.dataPacket_p)-1): + self.wallTime_e[self.dataPacket_p[i]:self.dataPacket_p[i+1]] = self.dataPacketTime_p[i] + self.wallTime_e[self.dataPacket_p[-1]:] = self.dataPacketTime_p[-1] + if not self.startTime and not norm: + self.startTime = self.wallTime_e[0] + self.wallTime_e -= self.startTime + logging.debug(f'wall time from {self.wallTime_e[0]} to {self.wallTime_e[-1]}') + + def read_event_stream(self): + self.tof_e = np.array(self.hdf['/entry1/Amor/detector/data/event_time_offset'][:])/1.e9 + self.pixelID_e = np.array(self.hdf['/entry1/Amor/detector/data/event_id'][:], dtype=int) + self.dataPacket_p = np.array(self.hdf['/entry1/Amor/detector/data/event_index'][:], dtype=np.uint64) + self.dataPacketTime_p = np.array(self.hdf['/entry1/Amor/detector/data/event_time_zero'][:], dtype=np.uint64)/1e9 + + def read_individual_header(self): + self.chopperDistance = float(np.take(self.hdf['entry1/Amor/chopper/pair_separation'], 0)) + self.detectorDistance = float(np.take(self.hdf['entry1/Amor/detector/transformation/distance'], 0)) + self.chopperDetectorDistance = self.detectorDistance-float(np.take(self.hdf['entry1/Amor/chopper/distance'], 0)) + self.tofCut = const.lamdaCut*self.chopperDetectorDistance/const.hdm*1.e-13 + + try: + self.mu = float(np.take(self.hdf['/entry1/Amor/master_parameters/mu/value'], 0)) + self.nu = float(np.take(self.hdf['/entry1/Amor/master_parameters/nu/value'], 0)) + self.kap = float(np.take(self.hdf['/entry1/Amor/master_parameters/kap/value'], 0)) + self.kad = float(np.take(self.hdf['/entry1/Amor/master_parameters/kad/value'], 0)) + self.div = float(np.take(self.hdf['/entry1/Amor/master_parameters/div/value'], 0)) + self.chopperSpeed = float(np.take(self.hdf['/entry1/Amor/chopper/rotation_speed/value'], 0)) + self.chopperPhase = float(np.take(self.hdf['/entry1/Amor/chopper/phase/value'], 0)) + except(KeyError, IndexError): + logging.warning(" using parameters from nicos cache") + year_date = str(self.start_date).replace('-', '/', 1) + cachePath = '/home/amor/nicosdata/amor/cache/' + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-mu/{year_date}')).split('\t')[-1] + self.mu = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-nu/{year_date}')).split('\t')[-1] + self.nu = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kap/{year_date}')).split('\t')[-1] + self.kap = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-kad/{year_date}')).split('\t')[-1] + self.kad = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-div/{year_date}')).split('\t')[-1] + self.div = float(value) + value = str(subprocess.getoutput(f'/usr/bin/grep "value" {cachePath}nicos-ch1_speed/{year_date}')).split('\t')[-1] + self.chopperSpeed = float(value) + self.chopperPhase = self.config.chopperPhase + self.tau = 30. / self.chopperSpeed + + if self.config.muOffset: + self.mu += self.config.muOffset + if self.config.mu: + self.mu = self.config.mu + if self.config.nu: + self.nu = self.config.nu + + # TODO: figure out real stop time.... + self.ctime=(self.hdf['/entry1/Amor/detector/data/event_time_zero'][-1] + - self.hdf['/entry1/Amor/detector/data/event_time_zero'][0]) / 1.e9 + self.fileDate = datetime.fromisoformat( self.hdf['/entry1/start_time'][0].decode('utf-8') ) + + def read_header_info(self): + # read general information and first data set + logging.info(f'# meta data from: {self.file_list[0]}') + self.hdf = h5py.File(self.file_list[0], 'r', swmr=True) + title = self.hdf['entry1/title'][0].decode('utf-8') + proposal_id = self.hdf['entry1/proposal_id'][0].decode('utf-8') + user_name = self.hdf['entry1/user/name'][0].decode('utf-8') + user_affiliation = 'unknown' + user_email = self.hdf['entry1/user/email'][0].decode('utf-8') + user_orcid = None + sampleName = self.hdf['entry1/sample/name'][0].decode('utf-8') + model = self.hdf['entry1/sample/model'][0].decode('utf-8') + instrumentName = 'Amor' + source = self.hdf['entry1/Amor/source/name'][0].decode('utf-8') + sourceProbe = 'neutron' + start_time = self.hdf['entry1/start_time'][0].decode('utf-8') + self.start_date = start_time.split(' ')[0] + if self.config.sampleModel: + model = self.config.sampleModel + # assembling orso header information + self.header.owner = fileio.Person( + name=user_name, + affiliation=user_affiliation, + contact=user_email, + ) + if user_orcid: + self.header.owner.orcid = user_orcid + self.header.experiment = fileio.Experiment( + title=title, + instrument=instrumentName, + start_date=self.start_date, + probe=sourceProbe, + facility=source, + proposalID=proposal_id + ) + self.header.sample = fileio.Sample( + name=sampleName, + model=model, + sample_parameters=None, + ) + self.header.measurement_scheme = 'angle- and energy-dispersive' + diff --git a/libeos/header.py b/libeos/header.py new file mode 100644 index 0000000..e46fc5e --- /dev/null +++ b/libeos/header.py @@ -0,0 +1,66 @@ +""" +Class to handle Orso header information that changes gradually during the reduction process. +""" + +import platform +import sys +from datetime import datetime + +from orsopy import fileio + +from . import __version__ + + +class Header: + """orso compatible output file header content""" + + def __init__(self): + self.owner = None + self.experiment = None + self.sample = None + self.measurement_instrument_settings = None + self.measurement_scheme = None + self.measurement_data_files = [] + self.measurement_additional_files = [] + + self.reduction = fileio.Reduction( + software = fileio.Software('eos', version=__version__), + call = ' '.join(sys.argv), + computer = platform.node(), + timestamp = datetime.now(), + creator = None, + corrections = ['histogramming in lambda and alpha_f', + 'gravity'], + ) + #------------------------------------------------------------------------------------------------- + def data_source(self): + return fileio.DataSource( + self.owner, + self.experiment, + self.sample, + fileio.Measurement( + instrument_settings = self.measurement_instrument_settings, + scheme = self.measurement_scheme, + data_files = self.measurement_data_files, + additional_files = self.measurement_additional_files, + ), + ) + #------------------------------------------------------------------------------------------------- + def columns(self): + cols = [ + fileio.Column('Qz', '1/angstrom', 'normal momentum transfer'), + fileio.Column('R', '', 'specular reflectivity'), + fileio.ErrorColumn(error_of='R', error_type='uncertainty', distribution='gaussian', value_is='sigma'), + fileio.ErrorColumn(error_of='Qz', error_type='resolution', distribution='gaussian', value_is='sigma'), + ] + return cols + + def orso_header(self, columns=None, extra_columns=[]): + """ + Generate ORSO header from a copy of this class' data. + """ + ds = fileio.DataSource.from_dict(self.data_source().to_dict()) + red = fileio.Reduction.from_dict(self.reduction.to_dict()) + if columns is None: + columns = self.columns() + return fileio.Orso(ds, red, columns+extra_columns) diff --git a/libeos/instrument.py b/libeos/instrument.py new file mode 100644 index 0000000..9a48817 --- /dev/null +++ b/libeos/instrument.py @@ -0,0 +1,68 @@ +""" +Classes describing the AMOR instrument configuration used during reduction. +""" + +import logging +import numpy as np + +from . import const + +class Detector: + nBlades = 14 # number of active blades in the detector + nWires = 32 # number of wires per blade + nStripes = 64 # number of stipes per blade + angle = np.deg2rad(5.1) # deg angle of incidence of the beam on the blades (def: 5.1) + dZ = 4.0*np.sin(angle) # mm height-distance of neighboring pixels on one blade + dX = 4.0*np.cos(angle) # mm depth-distance of neighboring pixels on one blace + bladeZ = 10.455 # mm distance between detector blades + zero = 0.5*nBlades*bladeZ # mm vertical center of the detector + distance = 4000. # mm distance from focal point to leading blade edge + +class Grid: + + def __init__(self, qResolution): + self.lamdaCut = const.lamdaCut + self.dldl = 0.005 # Delta lambda / lambda + self.qResolution = qResolution + + def q(self): + resolutions = [0.005, 0.01, 0.02, 0.025, 0.04, 0.05, 0.1, 1] + a, b = np.histogram([self.qResolution], bins = resolutions) + dqdq = np.matmul(b[:-1],a) + if dqdq != self.qResolution: + logging.info(f'# changed resolution to {dqdq}') + qq = 0.01 + # linear up to qq + q_grid = np.arange(0, qq, qq*dqdq) + # exponential from qq on + q_grid = np.append(q_grid, qq*(1.+dqdq)**np.arange(int(np.log(0.3/qq)/np.log(1+dqdq)))) + return q_grid + + def lamda(self): + lamdaMax = 16 + lamdaMin = self.lamdaCut + lamda_grid = lamdaMin*(1+self.dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+self.dldl)+1)) + return lamda_grid + + def z(self): + return np.arange(Detector.nBlades*Detector.nWires+1) + + def lz(self): + return np.ones(( np.shape(self.lamda()[:-1])[0], np.shape(self.z()[:-1])[0] )) + + def delta(self, detectorDistance): + # unused for now + bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*Detector.bladeZ / detectorDistance) ) + blade_grid = np.arctan( np.arange(33) * Detector.dZ / ( detectorDistance + np.arange(33) * Detector.dX) ) + blade_grid = np.rad2deg(blade_grid) + stepWidth = blade_grid[1] - blade_grid[0] + blade_grid = blade_grid - 0.2 * stepWidth + + delta_grid = [] + for b in np.arange(Detector.nBlades-1): + delta_grid = np.concatenate((delta_grid, blade_grid), axis=None) + blade_grid = blade_grid + bladeAngle + delta_grid = delta_grid[delta_grid0, filter_q, False) + if len(filter_q[filter_q]) > 0: + scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) + else: + logging.warning(f'# automatic scaling not possible') + scale = 1. + else: + filter_q = np.where(np.isnan(pR_q*R_q), False, True) + filter_q = np.where(R_q>0, filter_q, False) + filter_q = np.where(pR_q>0, filter_q, False) + if len(filter_q[filter_q]) > 0: + scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ + / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) + else: + logging.warning(f'# automatic scaling not possible') + scale = 1. + R_q /= scale + dR_q /= scale + logging.debug(f'# scaling factor = {scale}') + + return R_q, dR_q + + def project_on_qz(self, q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz): + q_q = self.grid.q() + mask_lzf = mask_lz.flatten() + q_lzf = q_lz.flatten()[mask_lzf] + R_lzf = R_lz.flatten()[mask_lzf] + dR_lzf = dR_lz.flatten()[mask_lzf] + dq_lzf = dq_lz.flatten()[mask_lzf] + norm_lzf = norm_lz.flatten()[mask_lzf] + + N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + + R_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf * R_lzf )[0] + R_q = R_q / N_q + + dR_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dR_lzf)**2 )[0] + dR_q = np.sqrt( dR_q ) / N_q + + # TODO: different error propagations for dR and dq! + N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf**2 )[0] + N_q = np.where(N_q > 0, N_q, np.nan) + dq_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dq_lzf)**2 )[0] + dq_q = np.sqrt( dq_q / N_q ) + + q_q = 0.5 * (q_q + np.roll(q_q, 1)) + + return q_q[1:], R_q, dR_q, dq_q + + def loadRqz(self, name): + fname = os.path.join(self.reader_config.dataPath, name) + if os.path.exists(fname): + fileName = fname + elif os.path.exists(f'{fname}.Rqz.ort'): + fileName = f'{fname}.Rqz.ort' + else: + sys.exit(f'### the background file \'{fname}\' does not exist! => stopping') + + q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) + + return q_q, Sq_q, dS_q, fileName + + def create_normalisation_map(self, short_notation): + dataPath = self.reader_config.dataPath + normalisation_list = expand_file_list(short_notation) + name = str(normalisation_list[0]) + for i in range(1, len(normalisation_list), 1): + name = f'{name}_{normalisation_list[i]}' + n_path = os.path.join(dataPath, f'{name}.norm') + if os.path.exists(n_path): + logging.info(f'# normalisation matrix: found and using {n_path}') + self.norm_lz = np.loadtxt(f'{dataPath}/{name}.norm') + with open(n_path, 'r') as fh: + fh.readline() + self.normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') + self.normAngle = float(fh.readline().split('= ')[1]) + for i, entry in enumerate(self.normFileList): + self.normFileList[i] = entry.split('/')[-1] + self.header.measurement_additional_files = self.normFileList + else: + logging.info(f'# normalisation matrix: using the files {normalisation_list}') + fromHDF = AmorData(header=self.header, + reader_config=self.reader_config, + config=self.experiment_config, + short_notation=short_notation, norm=True) + self.normAngle = fromHDF.nu - fromHDF.mu + lamda_e = fromHDF.lamda_e + detZ_e = fromHDF.detZ_e + self.norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (self.grid.lamda(), self.grid.z())) + self.norm_lz = np.where(self.norm_lz>0, self.norm_lz, np.nan) + # correct for the SM reflectivity + lamda_l = self.grid.lamda() + theta_z = self.normAngle + fromHDF.delta_z + lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T + theta_lz = self.grid.lz()*theta_z + qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz + Rsm_lz = np.ones(np.shape(qz_lz)) + Rsm_lz = np.where(qz_lz>0.0217, 1-(qz_lz-0.0217)*(0.0625/0.0217), Rsm_lz) + Rsm_lz = np.where(qz_lz>0.0217*5, np.nan, Rsm_lz) + self.norm_lz = self.norm_lz / Rsm_lz + if len(lamda_e) > 1e6: + head = ('normalisation matrix based on the measurements\n' + f'{fromHDF.file_list}\n' + f'nu - mu = {self.normAngle}\n' + f'shape= {np.shape(self.norm_lz)} (lambda, z)\n' + f'measured at mu = {fromHDF.mu:6.3f} deg\n' + f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') + head = head.replace('../', '') + head = head.replace('./', '') + head = head.replace('raw/', '') + np.savetxt(f'{dataPath}/{name}.norm', self.norm_lz, header = head) + self.normFileList = fromHDF.file_list + self.header.reduction.corrections.append('normalisation with \'additional files\'') + + def project_on_lz(self, fromHDF, norm_lz, normAngle, lamda_e, detZ_e): + # projection on lambda-z-grid + lamda_l = self.grid.lamda() + theta_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z + lamda_lz = (self.grid.lz().T*lamda_l[:-1]).T + theta_lz = self.grid.lz()*theta_z + + thetaN_z = fromHDF.delta_z + normAngle + thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z + thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) + + mask_lz = np.where(np.isnan(norm_lz), False, True) + mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) + mask_lz = np.logical_and(mask_lz, np.where(np.absolute(theta_lz)>5e-3, True, False)) + if self.reduction_config.thetaRange[1]<12: + mask_lz = np.logical_and(mask_lz, np.where(theta_lz >= self.reduction_config.thetaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(theta_lz <= self.reduction_config.thetaRange[1], True, False)) + if self.reduction_config.thetaRangeR[1]<12: + t0 = fromHDF.nu - fromHDF.mu + mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= self.reduction_config.thetaRangeR[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= self.reduction_config.thetaRangeR[1], True, False)) + if self.experiment_config.lambdaRange[1]<15: + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= self.experiment_config.lambdaRange[0], True, False)) + mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= self.experiment_config.lambdaRange[1], True, False)) + + # gravity correction + #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) + theta_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) + + z_z = enumerate(theta_z) + qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz + int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, self.grid.z())) + # cut normalisation sample horizon + int_lz = np.where(mask_lz, int_lz, np.nan) + thetaF_lz = np.where(mask_lz, theta_lz, np.nan) + + ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) + err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) + + res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(theta_z)[0])) * 0.022**2 + res_lz = res_lz + (0.008/theta_lz)**2 + res_lz = qz_lz * np.sqrt(res_lz) + + return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz + + + @staticmethod + def histogram2d_lz(lamda_e, detZ_e, bins): + """ + Perform binning operation equivalent to numpy bin2d for the sepcial case + of the second dimension using integer positions (pre-defined pixels). + Based on the devide_bin algorithm below. + """ + dimension = bins[1].shape[0]-1 + if not (np.array(bins[1])==np.arange(dimension+1)).all(): + raise ValueError("histogram2d_lz requires second bin dimension to be contigous integer range") + binning = AmorReduction.devide_bin(lamda_e, detZ_e.astype(np.int64), bins[0], dimension) + return np.array(binning), bins[0], bins[1] + + @staticmethod + def devide_bin(lambda_e, position_e, lamda_edges, dimension): + ''' + Use a divide and conquer strategy to bin the data. For the actual binning the + numpy bincount function is used, as it is much faster than histogram for + counting of integer values. + + :param lambda_e: Array of wavelength for each event + :param position_e: Array of positional indices for each event + :param lamda_edges: The edges of bins to be used for the histogram + :param dimension: position number of buckets in output arrray + + :return: 2D list of dimensions (lambda, x) of counts + ''' + if len(lambda_e)==0: + # no more events in range, return empty bins + return [np.zeros(dimension, dtype=int).tolist()]*(len(lamda_edges)-1) + if len(lamda_edges)==2: + # deepest recursion reached, all items should be within the two ToF edges + return [np.bincount(position_e, minlength=dimension).tolist()] + # split all events into two time of flight regions + split_idx = len(lamda_edges)//2 + left_region = lambda_e 2 tau - if True: - filter_e = (tof_e <= 2*self.tau) - tof_e = tof_e[filter_e] - pixelID_e = pixelID_e[filter_e] - wallTime_e = wallTime_e[filter_e] - if np.shape(filter_e)[0]-np.shape(tof_e)[0] > 0.5 : - print(f'# strange times: {np.shape(filter_e)[0]-np.shape(tof_e)[0]}') - tof_e = np.remainder( tof_e - self.tofCut + self.tau, self.tau) + self.tofCut # tof shifted to 1 frame - tof_e = tof_e + self.tau * clas.chopperPhaseOffset / 180. # correction for time offset between chopper pulse and tof zero - - # resolve pixel ID into y and z indicees, x position and angle - (detY_e, detZ_e, detXdist_e, delta_e) = pixelLookUp[np.int_(pixelID_e)-1,:].T - - # define mask and filter y range - mask_e = (clas.yRange[0] <= detY_e) & (detY_e <= clas.yRange[1]) - - # correct tof for beam size effect at chopper: t_cor = (delta / 180 deg) * tau - # TODO: check for correctness - if not clas.offSpecular: - tof_e -= ( delta_e / 180. ) * self.tau - - # lambda - lamda_e = 1.e13 * tof_e * defs.hdm / (self.chopperDetectorDistance + detXdist_e) - self.lamdaMax = defs.lamdaCut + 1.e13 * self.tau * defs.hdm / (self.chopperDetectorDistance + 124.) - mask_e = np.logical_and(mask_e, (clas.lambdaRange[0] <= lamda_e) & (lamda_e <= clas.lambdaRange[1])) - - # alpha_f - alphaF_e = self.nu - self.mu + delta_e - - # q_z - if clas.offSpecular: - alphaI = self.kap + self.kad + self.mu - qz_e = 2 * np.pi * ( np.sin( np.deg2rad(alphaF_e) ) + np.sin( np.deg2rad( alphaI ) ) ) / lamda_e - qx_e = 2 * np.pi * ( np.cos( np.deg2rad(alphaF_e) ) - np.cos( np.deg2rad( alphaI ) ) ) / lamda_e - header.measurement_scheme = 'energy-dispersive', - else: - qz_e = 4 * np.pi * np.sin( np.deg2rad(alphaF_e) ) / lamda_e - # qx_e = 0. - header.measurement_scheme = 'angle- and energy-dispersive' - - # filter q_z range - if clas.qzRange[1] < 0.3 and not norm: - mask_e = np.logical_and(mask_e, (clas.qzRange[0] <= qz_e) & (qz_e <= clas.qzRange[1])) - - self.detZ_e = detZ_e[mask_e] - self.lamda_e = lamda_e[mask_e] - self.wallTime_e = wallTime_e[mask_e] - - print(f'# number of events: total = {totalNumber:7d}, filtered = {np.shape(self.lamda_e)[0]:7d}') -#===================================================================================================== -class Detector: - def __init__(self): - self.nBlades = 14 # number of active blades in the detector - self.nWires = 32 # number of wires per blade - self.nStripes = 64 # number of stipes per blade - self.angle = np.deg2rad(5.1) # deg angle of incidence of the beam on the blades (def: 5.1) - self.dZ = 4.0 * np.sin(self.angle) # mm height-distance of neighboring pixels on one blade - self.dX = 4.0 * np.cos(self.angle) # mm depth-distance of neighboring pixels on one blace - self.bladeZ = 10.455 # mm distance between detector blades - self.zero = 0.5 * self.nBlades * self.bladeZ # mm vertical center of the detector - self.distance = 4000. # mm distance from focal point to leading blade edge -#===================================================================================================== -def expand_file_list(short_notation): - '''Evaluate string entry for file number lists''' - #log().debug('Executing get_flist') - file_list=[] - for i in short_notation.split(','): - if '-' in i: - if ':' in i: - step = i.split(':', 1)[1] - file_list += range(int(i.split('-', 1)[0]), int((i.rsplit('-', 1)[1]).split(':', 1)[0])+1, int(step)) - else: - step = 1 - file_list += range(int(i.split('-', 1)[0]), int(i.split('-', 1)[1])+1, int(step)) - else: - file_list += [int(i)] - - return sorted(file_list) -#===================================================================================================== -def normalisation_map(short_notation): - fromHDF = AmorData() - normalisation_list = expand_file_list(short_notation) - name = str(normalisation_list[0]) - for i in range(1, len(normalisation_list), 1): - name = f'{name}_{normalisation_list[i]}' - if os.path.exists(f'{clas.dataPath}/{name}.norm'): - print(f'# normalisation matrix: found and using {clas.dataPath}/{name}.norm') - norm_lz = np.loadtxt(f'{clas.dataPath}/{name}.norm') - fh = open(f'{clas.dataPath}/{name}.norm', 'r') - fh.readline() - normFileList = fh.readline().split('[')[1].split(']')[0].replace('\'', '').split(', ') - normAngle = float(fh.readline().split('= ')[1]) - fh.close() - for i, entry in enumerate(normFileList): - normFileList[i] = entry.split('/')[-1] - header.measurement_additional_files = normFileList - else: - print(f'# normalisation matrix: using the files {normalisation_list}') - fromHDF.read_data(short_notation, norm=True) - normAngle = fromHDF.nu - fromHDF.mu - lamda_e = fromHDF.lamda_e - detZ_e = fromHDF.detZ_e - norm_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (grid.lamda(), grid.z())) - norm_lz = np.where(norm_lz>0, norm_lz, np.nan) - # correct for the SM reflectivity - lamda_l = grid.lamda() - theta_z = normAngle + fromHDF.delta_z - lamda_lz = (grid.lz().T*lamda_l[:-1]).T - theta_lz = grid.lz()*theta_z - qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz - Rsm_lz = np.ones(np.shape(qz_lz)) - Rsm_lz = np.where(qz_lz>0.0217, 1-(qz_lz-0.0217)*(0.0625/0.0217), Rsm_lz) - Rsm_lz = np.where(qz_lz>0.0217*5, np.nan, Rsm_lz) - norm_lz = norm_lz / Rsm_lz - if len(lamda_e) > 1e6: - head = ('normalisation matrix based on the measurements\n' - f'{fromHDF.file_list}\n' - f'nu - mu = {normAngle}\n' - f'shape= {np.shape(norm_lz)} (lambda, z)\n' - f'measured at mu = {fromHDF.mu:6.3f} deg\n' - f'N(l_lambda, z) = theta(z) / sum_i=-1..1 I(l_lambda+i, z)') - head = head.replace('../', '') - head = head.replace('./', '') - head = head.replace('raw/', '') - np.savetxt(f'{clas.dataPath}/{name}.norm', norm_lz, header = head) - normFileList = fromHDF.file_list - return norm_lz, normAngle, normFileList -#===================================================================================================== -def output_format_list(outputFormat): - format_list = [] - if 'ort' in outputFormat or 'Rqz.ort' in outputFormat or 'Rqz' in outputFormat: - format_list.append('Rqz.ort') - if 'ort' in outputFormat or 'Rlt.ort' in outputFormat or 'Rlt' in outputFormat: - format_list.append('Rlt.ort') - if 'orb' in outputFormat or 'Rqz.orb' in outputFormat or 'Rqz' in outputFormat: - format_list.append('Rqz.orb') - if 'orb' in outputFormat or 'Rlt.orb' in outputFormat or 'Rlt' in outputFormat: - format_list.append('Rlt.orb') - - return sorted(format_list, reverse=True) -#===================================================================================================== -def project_on_lz(fromHDF, norm_lz, normAngle, lamda_e, detZ_e): - # projection on lambda-z-grid - lamda_l = grid.lamda() - theta_z = fromHDF.nu - fromHDF.mu + fromHDF.delta_z - lamda_lz = (grid.lz().T*lamda_l[:-1]).T - theta_lz = grid.lz()*theta_z - - thetaN_z = fromHDF.delta_z + normAngle - thetaN_lz = np.ones(np.shape(norm_lz))*thetaN_z - thetaN_lz = np.where(np.absolute(thetaN_lz)>5e-3, thetaN_lz, np.nan) - - mask_lz = np.where(np.isnan(norm_lz), False, True) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(thetaN_lz)>5e-3, True, False)) - mask_lz = np.logical_and(mask_lz, np.where(np.absolute(theta_lz)>5e-3, True, False)) - if clas.thetaRange[1]<12: - mask_lz = np.logical_and(mask_lz, np.where(theta_lz >= clas.thetaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(theta_lz <= clas.thetaRange[1], True, False)) - if clas.thetaRangeR[1]<12: - t0 = fromHDF.nu - fromHDF.mu - mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 >= clas.thetaRangeR[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(theta_lz-t0 <= clas.thetaRangeR[1], True, False)) - if clas.lambdaRange[1]<15: - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz >= clas.lambdaRange[0], True, False)) - mask_lz = np.logical_and(mask_lz, np.where(lamda_lz <= clas.lambdaRange[1], True, False)) - - # gravity correction - #theta_lz += np.rad2deg( np.arctan( 3.07e-10 * (fromHDF.detectorDistance + detXdist_e) * lamda_lz**2 ) ) - theta_lz += np.rad2deg( np.arctan( 3.07e-10 * fromHDF.detectorDistance * lamda_lz**2 ) ) - - z_z = enumerate(theta_z) - qz_lz = 4.0*np.pi * np.sin(np.deg2rad(theta_lz)) / lamda_lz - int_lz, bins_l, bins_z = np.histogram2d(lamda_e, detZ_e, bins = (lamda_l, grid.z())) - # cut normalisation sample horizon - int_lz = np.where(mask_lz, int_lz, np.nan) - thetaF_lz = np.where(mask_lz, theta_lz, np.nan) - - ref_lz = (int_lz * np.absolute(thetaN_lz)) / (norm_lz * np.absolute(thetaF_lz)) - err_lz = ref_lz * np.sqrt( 1/(int_lz+.1) + 1/norm_lz ) - - res_lz = np.ones((np.shape(lamda_l[:-1])[0], np.shape(theta_z)[0])) * 0.022**2 - res_lz = res_lz + (0.008/theta_lz)**2 - res_lz = qz_lz * np.sqrt(res_lz) - - return qz_lz, ref_lz, err_lz, res_lz, lamda_lz, theta_lz, int_lz, mask_lz -#===================================================================================================== -def project_on_qz(q_lz, R_lz, dR_lz, dq_lz, norm_lz, mask_lz): - q_q = grid.q() - mask_lzf = mask_lz.flatten() - q_lzf = q_lz.flatten()[mask_lzf] - R_lzf = R_lz.flatten()[mask_lzf] - dR_lzf = dR_lz.flatten()[mask_lzf] - dq_lzf = dq_lz.flatten()[mask_lzf] - norm_lzf = norm_lz.flatten()[mask_lzf] - - N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - - R_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf * R_lzf )[0] - R_q = R_q / N_q - - dR_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dR_lzf)**2 )[0] - dR_q = np.sqrt( dR_q ) / N_q - - # TODO: different error propagations for dR and dq! - N_q = np.histogram(q_lzf, bins = q_q, weights = norm_lzf**2 )[0] - N_q = np.where(N_q > 0, N_q, np.nan) - dq_q = np.histogram(q_lzf, bins = q_q, weights = (norm_lzf * dq_lzf)**2 )[0] - dq_q = np.sqrt( dq_q / N_q ) - - q_q = 0.5 * (q_q + np.roll(q_q, 1)) - - return q_q[1:], R_q, dR_q, dq_q -#===================================================================================================== -def autoscale(q_q, R_q, dR_q, pR_q=[], pdR_q=[]): - if len(pR_q) == 0: - filter_q = np.where((clas.autoscale[0]<=q_q)&(q_q<=clas.autoscale[1]), True, False) - filter_q = np.where(dR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**2/dR_q[filter_q]) / np.sum(R_q[filter_q]/dR_q[filter_q]) - else: - print(f'# automatic scaling not possible') - scale = 1. - else: - filter_q = np.where(np.isnan(pR_q*R_q), False, True) - filter_q = np.where(R_q>0, filter_q, False) - filter_q = np.where(pR_q>0, filter_q, False) - if len(filter_q[filter_q]) > 0: - scale = np.sum(R_q[filter_q]**3 * pR_q[filter_q] / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) \ - / np.sum(R_q[filter_q]**2 * pR_q[filter_q]**2 / (dR_q[filter_q]**2 * pdR_q[filter_q]**2)) - else: - print(f'# automatic scaling not possible') - scale = 1. - R_q /= scale - dR_q /= scale - #print(f'# scaling factor = {scale}') - - return R_q, dR_q -#===================================================================================================== -def loadRqz(name): - - if os.path.exists(f'{clas.dataPath}/{name}'): - fileName = f'{clas.dataPath}/{name}' - elif os.path.exists(f'{clas.dataPath}/{name}.Rqz.ort'): - fileName = f'{clas.dataPath}/{name}.Rqz.ort' - else: - sys.exit(f'### the background file \'{clas.dataPath}/{name}\' does not exist! => stopping') - - q_q, Sq_q, dS_q = np.loadtxt(fileName, usecols=(0, 1, 2), comments='#', unpack=True) - - return q_q, Sq_q, dS_q, fileName -#===================================================================================================== -class Grid: - - def __init__(self): - self.det = Detector() - self.lamdaCut = Defs.lamdaCut - self.dldl = 0.005 # Delta lambda / lambda - - def q(self): - resolutions = [0.005, 0.01, 0.02, 0.025, 0.04, 0.05, 0.1, 1] - a, b = np.histogram([clas.qResolution], bins = resolutions) - dqdq = np.matmul(b[:-1],a) - if dqdq != clas.qResolution: - print(f'# changed resolution to {dqdq}') - qq = 0.01 - # linear up to qq - q_grid = np.arange(0, qq, qq*dqdq) - # exponential from qq on - q_grid = np.append(q_grid, qq*(1.+dqdq)**np.arange(int(np.log(0.3/qq)/np.log(1+dqdq)))) - return q_grid - - def lamda(self): - lamdaMax = 16 - lamdaMin = self.lamdaCut - lamda_grid = lamdaMin*(1+self.dldl)**np.arange(int(np.log(lamdaMax/lamdaMin)/np.log(1+self.dldl)+1)) - return lamda_grid - - def z(self): - return np.arange(self.det.nBlades*self.det.nWires+1) - - def lz(self): - return np.ones(( np.shape(self.lamda()[:-1])[0], np.shape(self.z()[:-1])[0] )) - - def delta(self, detectorDistance): - # unused for now - bladeAngle = np.rad2deg( 2. * np.arcsin(0.5*self.det.bladeZ / detectorDistance) ) - blade_grid = np.arctan( np.arange(33) * self.det.dZ / ( detectorDistance + np.arange(33) * self.det.dX) ) - blade_grid = np.rad2deg(blade_grid) - stepWidth = blade_grid[1] - blade_grid[0] - blade_grid = blade_grid - 0.2 * stepWidth - - delta_grid = [] - for b in np.arange(self.det.nBlades-1): - delta_grid = np.concatenate((delta_grid, blade_grid), axis=None) - blade_grid = blade_grid + bladeAngle - delta_grid = delta_grid[delta_grid