diff --git a/direct_inelastic/ISIS/DG_monovan.py b/direct_inelastic/ISIS/DG_monovan.py index 7545ee9..9fcd63b 100644 --- a/direct_inelastic/ISIS/DG_monovan.py +++ b/direct_inelastic/ISIS/DG_monovan.py @@ -20,6 +20,7 @@ t = time.time() #start the clock +#!begin_params #============================User Input=========================== # Assumes that the thickness of the vanadium and the sample are the same monovan = 59138 #mono van run @@ -37,22 +38,35 @@ # MAPS monitors: m2spec=41475, m3spec=41476 - fixei = False # LET monitors: m2spec=98310, m3spec= None - fixei = True -config['default.instrument'] = 'INSTRUMENT_NAME' +inst = 'INSTRUMENT_NAME' cycle = 'CYCLE_ID' -m2spec = 69636 #specID of monitor2 (post monochromator) -m3spec = 69640 #specID of monitor2 (post monochromator) fixei = False #true for LET since no m3 powdermap = 'RINGS_MAP_XML' #rings mapping file monovan_range = 0.1 #integration range (+/- fraction of Ei) min_theta = 5. #minimum theta to avoid low Q detectors (SANS) idebug = False #keep workspaces #================================================================= -inst = config['default.instrument'] +#!end_params +config['default.instrument'] = inst if inst == 'MARI': source = 'Moderator' -else: + m2spec = 2 # specID of monitor2 (pre-sample) + m3spec = 3 # specID of monitor3 (post-sample) +elif inst == 'MERLIN': + source = 'undulator' + m2spec = 69636 # specID of monitor2 (pre-sample) + m3spec = 69640 # specID of monitor3 (post-sample) +elif inst == 'MAPS': source = 'undulator' + m2spec = 36867 # specID of monitor2 (pre-sample) + m3spec = 36868 # specID of monitor3 (post-sample) +elif inst == 'LET': + source = 'undulator' + m2spec = 98310 # specID of monitor2 (pre-sample) + m3spec = None # specID of monitor3 (post-sample) +else: + raise RuntimeError(f'Unrecognised instrument: {inst}') if cycle.startswith('20'): cycle = cycle[2:] diff --git a/direct_inelastic/ISIS/DG_reduction.py b/direct_inelastic/ISIS/DG_reduction.py index f8aebab..4a46b67 100644 --- a/direct_inelastic/ISIS/DG_reduction.py +++ b/direct_inelastic/ISIS/DG_reduction.py @@ -17,26 +17,38 @@ # import mantid algorithms, numpy and matplotlib from mantid.simpleapi import * from mantid.api import AnalysisDataService as ADS +from mantid.kernel.funcinspect import lhs_info import numpy as np import os, sys import time from importlib import reload +#!begin_params #=======================User Inputs======================= powder = True # powder or 1to1 map sumruns = False # set to True to sum sample runs sample = [93338] # sample runs (list) sample_bg = 93329 # single background run +sample_cd = None # Cadmium run for instrument background subtraction wv_file = 'WV_91329.txt' # white vanadium integral file (mandatory) -Ei_list = [3.71,1.035,1.775] # incident energies (Ei) - from PyChop +Ei_list = ['auto'] # incident energies (Ei) - from PyChop Erange = [-0.8,0.0025,0.8] # energy transfer range to output in fractions of Ei trans = [0.95,0.95,0.95] # elastic line transmission factors for each Ei mask = 'MASK_FILE_XML' # hard mask # what to do if see run with same angle as previous run (and sumruns==False) # - 'replace': sum and replace first nxspe file created # - 'ignore': ignore and create an nxspe for each run -# - 'accumulate': sum all runs with this angle when creating new nxspe same_angle_action = 'ignore' +# to enable creating multiple reduced data files from a single +# "continuous scan" set the following variables to a valid log "block" +# name and the bin size (width) to something larger than zero. +# An optional unit can be given. +# If sumruns=True, the sample list is summed and then divided +# If sumruns=False and multiple sample runs are listed, an error is raised. +# This is also only compatible with same_angle_action='ignore' +cs_block = None +cs_block_unit = '' +cs_bin_size = 0 #======================================================== #==================Absolute Units Inputs================= @@ -52,7 +64,7 @@ fixei = True # True for LET since no monitor 3 powdermap = 'RINGS_MAP_XML' # rings mapping file - must be .xml format file_wait = 30 # wait for data file to appear (seconds) -keepworkspaces = True # should be false for Horace scans +keepworkspaces = False # should be false for Horace scans saveformat = '.nxspe' # format of output, '.nxspe', '.nxs' QENS = False # output Q-binned data for QENS data analysis '_red.nxs' Qbins = 20 # approximate number of Q-bins for QENS @@ -61,7 +73,11 @@ save_dir = f'/instrument/{inst}/RBNumber/USER_RB_FOLDER' # Set to None to avoid reseting psi_motor_name = 'rot' # name of the rotation motor in the logs angles_workspace = 'angles_ws' # name of workspace to store previously seen angles +sumruns_savemem = False # Compresses event in summed ws to save memory + # (causes some loss of data so cannot use event filtering) +cs_smidge = 0.001 # Tolerance on continuous scan bins size #======================================================== +#!end_params # ==============================setup directroties================================ config['default.instrument'] = inst @@ -100,7 +116,7 @@ ws_list = ADS.getObjectNames() # Try to load subsidiary utilities -sys.path.append(instfiles) +sys.path.append(os.path.dirname(__file__)) try: from reduction_utils import * utils_loaded = True @@ -108,6 +124,9 @@ utils_loaded = False def rename_existing_ws(*a, **k): raise RuntimeError('Unable to load utils') def remove_extra_spectra_if_mari(*a, **k): raise RuntimeError('Unable to load utils') + def autoei(*a, **k): + raise RuntimeError('Cannot use Auto-Ei: ' \ + 'You need to copy the reduction_utils.py file somewhere on the Python path') #======================================================== # Helper functions @@ -125,21 +144,32 @@ def tryload(irun): # loops till data file appears Pause(file_wait) continue break - if inst == 'MARI': - remove_extra_spectra_if_mari('ws') return mtd['ws'] -def load_sum(run_list): +def load_sum(run_list, block_name=None): for ii, irun in enumerate(run_list): tryload(irun) if ii == 0: w_buf = CloneWorkspace('ws') w_buf_monitors = CloneWorkspace('ws_monitors') + if block_name: + bval = mtd['ws'].getRun().getLogData(block_name).filtered_value print(f'run #{irun} loaded') else: w_buf = Plus('w_buf', 'ws') w_buf_monitors = Plus('w_buf_monitors', 'ws_monitors') + if block_name: + bv2 = mtd['ws'].getRun().getLogData(block_name).filtered_value + bval = np.concatenate((bval, bv2)) print(f'run #{irun} added') + ADS.remove('ws') + ADS.remove('ws_monitors') + wsout_name = lhs_info('names')[0] + if wsout_name != 'w_buf': + RenameWorkspace('w_buf_monitors', wsout_name+'_monitors') + RenameWorkspace('w_buf', wsout_name) + return (mtd[wsout_name], bval) if block_name else mtd[wsout_name] + #======================================================== @@ -151,6 +181,8 @@ def load_sum(run_list): sample = [sample] if sample_bg is not None and (isinstance(sample_bg, str) or not hasattr(sample_bg, '__iter__')): sample_bg = [sample_bg] +if sample_cd is not None and (isinstance(sample_cd, str) or not hasattr(sample_cd, '__iter__')): + sample_cd = [sample_cd] #==================================load hard mask================================ if mask is None: @@ -171,9 +203,71 @@ def load_sum(run_list): else: print(f'{inst}: Using previously loaded white vanadium - {wv_file}') -# ===============================load monovan file=============================== +# =======================load background runs and sum========================= +if sample_cd is not None: + ws_cd = load_sum(sample_cd) +if sample_bg is not None: + ws_bg = load_sum(sample_bg) + if sample_cd is not None: + ws_bg = ws_bg - ws_cd + ws_bg = NormaliseByCurrent(ws_bg) + +# ==========================continuous scan stuff============================= +if cs_block and cs_bin_size > 0: + if len(sample) > 1 and not sumruns: + raise RuntimeError('Continous scans for multiple (non-summed) files not supported') + if same_angle_action.lower() != 'ignore': + raise RuntimeError(f'Continous scans not compatible with same_angle_action="{same_angle_action}"') + # Sets sumruns false so don't redo the sum below. Instead sum here to get filtered block values. + # This is due to a bug in the Mantid "Plus" algorithm, which gives incorrect filtered values + # https://github.com/mantidproject/mantid/issues/36194 + sumruns = False + ws_full, bval = load_sum(sample, cs_block) + ws_monitors = CloneWorkspace('ws_full_monitors') + bval_range = max(bval) - min(bval) + bval_nbins = int(bval_range / cs_bin_size) + bval_remainder = bval_range - cs_bin_size * bval_nbins + # Overrides the sample list in order to hijack the loop over run numbers below + irun_orig = sample[0] + sample = [x*cs_bin_size + min(bval) for x in range(bval_nbins)] + unit = cs_block_unit + if not unit: + unit = ws_full.getRun().getLogData(cs_block).units + print(f'{cs_block} = {min(bval):.1f} {unit} to {max(bval):.1f} {unit}') + print(f'... filtering in {cs_bin_size:.1f} {unit} steps') + print(f'... N={bval_nbins} bins with {bval_remainder:.2f} {unit} remainder') + +# =======================sum sample runs if required========================= +sumsuf = sumruns and len(sample) > 1 +if sumruns: + ws = load_sum(sample) + sample = [sample[0]] + +# =============================auto-Ei stuff================================= +is_auto = lambda x: isinstance(x, str) and 'auto' in x.lower() +if is_auto(Ei_list) or hasattr(Ei_list, '__iter__') and is_auto(Ei_list[0]): + use_auto_ei = True + try: + Ei_list = autoei(ws) + except NameError: + fn = str(sample[0]) + if not fn.startswith(inst[:3]): fn = f'{inst[:3]}{fn}' + if fn.endswith('.raw'): fn = fn[:-4] + if fn[-4:-2] == '.s': fn = fn[:-4] + '.n0' + fn[-2:] + if not fn.endswith('.nxs') and fn[-5:-2] != '.n0': fn += '.nxs' + Ei_list = autoei(LoadNexusMonitors(fn, OutputWorkspace='ws_tmp_mons')) + print(f"Automatically determined Ei's: {Ei_list}") + if len(trans) < len(Ei_list): + print(f'{inst}: WARNING - not enough transmision values for auto-Ei. ' \ + 'Extending list with last (end) value') + trans += [trans[-1]]*(len(Ei_list) - len(trans)) +else: + use_auto_ei = False + +# ============================load monovan file============================== mv_fac = [] if mv_file is not None and monovan_mass is not None: + use_auto_ei = False if mv_file not in ws_list: print(f'{inst}: Loading monovan calibration factors - {mv_file}') LoadAscii(mv_file,OutputWorkspace=mv_file) @@ -183,7 +277,7 @@ def load_sum(run_list): # check that monovan is compatible with Ei_list) Ei_diff = sum([x-y for x,y in zip(mv_eis,Ei_list)]) if (abs(Ei_diff) > 0): - print('----ERROR: Monovan file Eis not compatible with Ei_list') + raise RuntimeError('----ERROR: Monovan file Eis not compatible with Ei_list') for Ei in Ei_list: ii = np.where(mv_eis == Ei) mvf = mv_cal[ii][0] @@ -195,19 +289,6 @@ def load_sum(run_list): print(f'{inst}: Skipping absolute calibration') mv_fac = [x/x for x in Ei_list] # monovan factors = 1 by default -# =======================load background runs and sum========================= -if sample_bg is not None: - load_sum(sample_bg) - ws_bg = CloneWorkspace('w_buf') - ws_bg = NormaliseByCurrent('ws_bg') - -# =======================sum sample runs if required========================= -if sumruns: - load_sum(sample) - ws = CloneWorkspace('w_buf') - ws_monitors = CloneWorkspace('w_buf_monitors') - sample = [sample[0]] - # =====================angles cache stuff==================================== if utils_loaded: build_angles_ws(sample, angles_workspace, psi_motor_name) @@ -226,24 +307,44 @@ def load_sum(run_list): for ss in nx_list: ADS.remove(ss) + # For continuous scans generate the workspace based on binned log values + # NB. we replaced the sample (run) list with list of block min bin boundaries + if cs_block and cs_bin_size > 0: + val = round(irun + cs_bin_size/2,1) + print(f"Filtering: {cs_block}={val:.1f} ± {round(cs_bin_size/2,2):.2f} {unit}") + ws = FilterByLogValue(ws_full, cs_block, irun, irun + cs_bin_size - cs_smidge) + # if this run is at an angle already seen and we are doing 'replace', skip if same_angle_action.lower() != 'ignore' and irun in runs_with_angles_already_seen: continue print('============') - if not sumruns: + if not sumruns and not cs_block: # Checks rotation angle if same_angle_action.lower() != 'ignore': runs_with_same_angles = get_angle(irun, angles_workspace, psi_motor_name, tryload) if len(runs_with_same_angles) > 1: - load_sum(runs_with_same_angles) - if same_angle_action.lower() == 'replace': - irun = runs_with_same_angles[0] - runs_with_angles_already_seen += runs_with_same_angles + ws = load_sum(runs_with_same_angles) + irun = runs_with_same_angles[0] + runs_with_angles_already_seen += runs_with_same_angles else: tryload(irun) print(f'Loading run# {irun}') + ws = NormaliseByCurrent('ws') + if sumruns and sumruns_savemem: + ws = CompressEvents(ws, Tolerance=1e-5) # Tolerance in microseconds + + # instrument geometry to work out ToF ranges + sampos = ws.getInstrument().getSample().getPos() + l1 = (sampos - ws.getInstrument().getSource().getPos()).norm() + l2 = (ws.getDetector(0).getPos() - sampos).norm() + + # Updates ei_loop if auto _and_ not using monovan + if use_auto_ei: + Ei_list = autoei(ws) + if len(trans) < len(Ei_list): trans += [trans[-1]]*(len(Ei_list) - len(trans)) + if len(mv_fac) < len(Ei_list): mv_fac += [1]*(len(Ei_list) - len(mv_fac)) # ============================= Ei loop ===================================== for ienergy in range(len(Ei_list)): @@ -253,15 +354,25 @@ def load_sum(run_list): mvf = mv_fac[ienergy] print(f'\n{inst}: Reducing data for Ei={Ei:.2f} meV') - ws_corrected = Scale('ws', 1 / tr, 'Multiply') + tofs = ws.readX(0) + tof_min = np.sqrt(l1**2 * 5.227e6 / Ei) - t_shift + tof_max = tof_min + np.sqrt(l2**2 * 5.226e6 / (Ei*(1-Erange[-1]))) + ws_rep = CropWorkspace(ws, max(min(tofs), tof_min), min(max(tofs), tof_max)) + + if sample_cd is not None: + print(f'... subtracting Cd background') + ws_rep = ws_rep - ws_cd + if sample_bg is not None: print(f'... subtracting background - transmission factor = {tr:.2f}') - ws_corrected = ws/tr - ws_bg + ws_rep = ws_rep/tr - ws_bg + else: + ws_rep = Scale('ws_rep', 1 / tr, 'Multiply') # normalise to WB vanadium and apply fixed mask print('... normalising/masking data') - ws_norm = Divide('ws_corrected', wv_file) # white beam normalisation - MaskDetectors(ws_norm,MaskedWorkspace=mask,ForceInstrumentMasking=True) + ws_rep = Divide('ws_rep', wv_file) # white beam normalisation + MaskDetectors(ws_rep, MaskedWorkspace=mask, ForceInstrumentMasking=True) # t2e section print('... t2e section') @@ -270,48 +381,52 @@ def load_sum(run_list): index = spectra.index(m2spec) m2pos = ws.detectorInfo().position(index)[2] - if inst == 'MARI' and utils_loaded and origEi < 4.01: - # Shifts data / monitors into second frame for MARI - ws_norm, ws_monitors = shift_frame_for_mari_lowE(origEi, wsname='ws_norm', wsmon='ws_monitors') - # this section shifts the time-of-flight such that the monitor2 peak # in the current monitor workspace (post monochromator) is at t=0 and L=0 # note that the offest is not the predicted value due to energy dependence of the source position - if m3spec is not None and not fixei: - (Ei,mon2_peak,_,_) = GetEi(ws_monitors,Monitor1Spec=m2spec,Monitor2Spec=m3spec,EnergyEstimate=Ei) + if m3spec is not None: + (Ei,mon2_peak,_,_) = GetEi(ws_monitors,Monitor1Spec=m2spec,Monitor2Spec=m3spec,EnergyEstimate=Ei,FixEi=fixei) print(f'... refined Ei={Ei:.2f} meV') else: (Ei,mon2_peak,_,_) = GetEi(ws_monitors,Monitor2Spec=m2spec,EnergyEstimate=Ei,FixEi=fixei) print(f'... m2 tof={mon2_peak:.2f} mus, m2 pos={m2pos:.2f} m') - ws_norm = ScaleX(ws_norm, Factor=-mon2_peak, Operation='Add', InstrumentParameter='DelayTime', Combine=True) - MoveInstrumentComponent(ws_norm, ComponentName=source, Z=m2pos, RelativePosition=False) + ws_rep = ScaleX(ws_rep, Factor=-mon2_peak, Operation='Add', InstrumentParameter='DelayTime', Combine=True) + MoveInstrumentComponent(ws_rep, ComponentName=source, Z=m2pos, RelativePosition=False) - ws_out = ConvertUnits(ws_norm, 'DeltaE', EMode='Direct', EFixed=Ei) - ws_out = Rebin(ws_out, [x*Ei for x in Erange], PreserveEvents=False) + ws_rep = ConvertUnits(ws_rep, 'DeltaE', EMode='Direct', EFixed=Ei) + ws_out = Rebin(ws_rep, [x*origEi for x in Erange], PreserveEvents=False) ws_out = DetectorEfficiencyCor(ws_out, IncidentEnergy=Ei) ws_out = CorrectKiKf(ws_out, Efixed=Ei, EMode='Direct') + ADS.remove('ws_rep') # monovan scaling if mv_file is not None: print(f'... applying mono van calibration factor {mvf:.1f}') ws_out = Scale('ws_out', mvf, 'Multiply') + # Sets output file name + if cs_block and cs_bin_size > 0: + ofile_prefix = f'{inst[:3]}{irun_orig}' + ofile_suffix = f"_{val:.1f}{unit}" + else: + ofile_prefix = f'{inst[:3]}{irun}' + ofile_suffix = '' + # rings grouping if desired - ofile_suffix='_1to1' if powder or inst == 'MARI' or QENS: ws_out=GroupDetectors(ws_out, MapFile=powdermap, Behaviour='Average') - ofile_suffix = '_powder' + ofile_suffix += '_powder' if inst == 'MARI' or QENS: ofile_suffix = '' print(f'... powder grouping using {powdermap}') - if sample_bg is not None and inst == 'MARI': - ofile_suffix += '_sub' + else: + ofile_suffix += '_1to1' # output nxspe file - ofile = f'{inst[:3]}{irun}_{origEi:g}meV{ofile_suffix}' + ofile = f'{ofile_prefix}_{origEi:g}meV{ofile_suffix}' # check elastic line (debug mode) if idebug: @@ -325,12 +440,12 @@ def load_sum(run_list): print(f'{inst}: Writing {ofile}{saveformat}') if saveformat.lower() == '.nxspe': SaveNXSPE('ws_out', ofile+saveformat, Efixed=Ei, KiOverKfScaling=True) - #if utils_loaded and not powder: - # copy_inst_info(ofile+saveformat, 'ws_out') # Copies instrument info for Horace + if utils_loaded: + copy_inst_info(ofile+saveformat, 'ws_out') # Copies instrument info for Horace elif saveformat.lower() == '.nxs': rmlogs = {'events_log', 'frame_log', 'good_frame_log', 'period_log', 'proton_charge', 'raw_events_log'} RemoveLogs('ws_out', KeepLogs=','.join(set(mtd['ws_out'].run().keys()).difference(rmlogs))) - SaveNexus('ws_out', ofile+saveformat) + SaveNexusProcessed('ws_out', ofile+saveformat, PreserveEvents=False, CompressNexus=True) if QENS: print('... outputting QENS "_red" format') theta = np.array([theta_range[0], theta_range[1]])*np.pi/180. diff --git a/direct_inelastic/ISIS/DG_whitevan.py b/direct_inelastic/ISIS/DG_whitevan.py index c4455ee..692209d 100644 --- a/direct_inelastic/ISIS/DG_whitevan.py +++ b/direct_inelastic/ISIS/DG_whitevan.py @@ -20,6 +20,7 @@ t = time.time() #start the clock +#!begin_params #=======================User Inputs====================== whitevan = 78715 # white vanadium run whitevan_bg = None # background for the white vanadium @@ -34,6 +35,7 @@ idebug = False # keep itermediate workspaces for debugging save_dir = f'/instrument/{inst}/RBNumber/USER_RB_FOLDER' # Set to None to avoid resetting #======================================================== +#!end_params config['default.instrument'] = inst if save_dir is not None: @@ -100,7 +102,7 @@ def load_sum(run_list): WV_normalised_integrals = Scale(WV_normalised_integrals, scale_factor, 'Multiply') # ===========================output integrals file================================ -wv_name = whitevan if (isinstance(whitevan, str) or len(whitevan)==1) else whitevan[0] +wv_name = whitevan if (isinstance(whitevan, (str, int, float)) or len(whitevan)==1) else whitevan[0] ofile = f'WV_{wv_name}.txt' print(f'WHITE_VAN {inst}: Saving integrals in {ofile}') print(f'... average intensity = {1/scale_factor:.2f}') diff --git a/direct_inelastic/ISIS/reduction_utils.py b/direct_inelastic/ISIS/reduction_utils.py new file mode 100644 index 0000000..12ae233 --- /dev/null +++ b/direct_inelastic/ISIS/reduction_utils.py @@ -0,0 +1,491 @@ +from mantid.simpleapi import * +from mantid.api import AnalysisDataService as ADS +import numpy as np +import re, os, sys, h5py +import json +import warnings +import importlib +import types +import scipy.optimize +from os.path import abspath, dirname +from mantid.kernel.funcinspect import lhs_info + +#======================================================== +# General utility functions + +def rename_existing_ws(ws_name): + # Gets an existing workspace and ensures we have its monitors for reduction + ws = CloneWorkspace(ws_name, OutputWorkspace='ws') + try: + mon_ws = ws.getMonitorWorkspace() + except RuntimeError: + specInfo = ws.spectrumInfo() + try: + mon_list = [i for i in range(ws.getNumberHistograms()) if specInfo.isMonitor(i)] + except RuntimeError: + mon_list = [] + if len(mon_list) > 0: + ExtractMonitors(ws_name, DetectorWorkspace='ws', MonitorWorkspace='ws_monitors') + else: + _get_mon_from_history(ws_name) + else: + CloneWorkspace(mon_ws, OutputWorkspace='ws_monitors') + return ws + +def _get_mon_from_history(ws_name): + # Tries to look in the history of a workspace for its original raw file + # loads the monitors from that raw file. + orig_file = None + for hist in mtd[ws_name].getHistory().getAlgorithmHistories(): + if hist.name().startswith('Load') and 'Filename' in [pp.name() for pp in hist.getProperties()]: + orig_file = hist.getPropertyValue('Filename') + break + if orig_file is None: + raise RuntimeError(f'Cannot find original file from workspace {ws_name} to load logs from') + try: + Load(orig_file, SpectrumMax=10, LoadMonitors=True, OutputWorkspace='tmp_mons') + except TypeError: + Load(orig_file, SpectrumMax=10, LoadMonitors='Separate', OutputWorkspace='tmp_mons') + ws_mon_name = f'{ws_name}_monitors' + RenameWorkspace('tmp_mons_monitors', ws_mon_name) + DeleteWorkspace('tmp_mons') + mtd[ws_name].setMonitorWorkspace(mtd[ws_mon_name]) + CloneWorkspace(ws_mon_name, OutputWorkspace='ws_monitors') + +def get_angle(irun, angle_workspace='angle_ws', psi_motor_name='rot', tryload=None): + # Checks if a workspace with previous angles exists and if we've seen this run before + if angle_workspace not in mtd: + CreateWorkspace(OutputWorkspace=angle_workspace, DataX=0, DataY=0) + prev_angles = {} + if mtd[angle_workspace].getRun().hasProperty('angles_seen'): + prev_angles = json.loads(mtd[angle_workspace].getRun().getProperty('angles_seen').value) + if irun in sum(prev_angles.values(), []): + angle = float([k for k, v in prev_angles.items() if irun in v][0]) + else: + if tryload is not None: + ws = tryload(irun) + else: + try: + ws = Load(orig_file, SpectrumMax=10, LoadMonitors=True, StoreInADS=False) + except TypeError: + ws = Load(orig_file, SpectrumMax=10, LoadMonitors='Separate', StoreInADS=False) + angle = ws.getRun().getLogData(psi_motor_name).value[-1] + # Reduce to 0.2 degree accuracy to check for equivalent angles + angle = np.round(angle * 5.) / 5. + print(f'Read logs from run {irun}, at rotation {angle} degrees') + angles = str(angle) + if angles not in prev_angles.keys(): + prev_angles[angles] = [] + if irun not in prev_angles[angles]: + prev_angles[angles].append(irun) + # Save previous angle information to workspace + loginfo = json.dumps(prev_angles) + mtd[angle_workspace].getRun().addProperty('angles_seen', loginfo, '', True) + return prev_angles[angles] + +def build_angles_ws(run_list, angle_workspace, psi_motor_name): + for irun in run_list: + # Just run through the list to build up list of angles in the angles_workspace + try: + get_angle(irun, angle_workspace=angle_workspace, psi_motor_name=psi_motor_name) + except: + pass + +#======================================================== +# Functions to copy instrument info needed by HORACE +# Resolution convolution if it exists in the raw file +def get_raw_file_from_ws(ws): + for alg in [h for h in ws.getHistory().getAlgorithmHistories() if 'Load' in h.name()]: + for prp in [a for a in alg.getProperties() if 'Filename' in a.name()]: + if re.search('[0-9]*.nxs', prp.value()) is not None: + return prp.value() + raise RuntimeError('Could not find raw NeXus file in workspace history') + + +def copy_inst_info(outfile, in_ws): + print(f'Copying Instrument Info to file {outfile}') + try: + raw_file_name = get_raw_file_from_ws(mtd[in_ws]) + except RuntimeError: + return + if not os.path.exists(outfile): + outfile = os.path.join(mantid.simpleapi.config['defaultsave.directory'], os.path.basename(outfile)) + with h5py.File(raw_file_name, 'r') as raw: + exclude = ['dae', 'detector_1', 'name'] + to_copy = [k for k in raw['/raw_data_1/instrument'] if not any([x in k for x in exclude])] + if 'aperture' not in to_copy and 'mono_chopper' not in to_copy: + return + n_spec = len(raw['/raw_data_1/instrument/detector_1/spectrum_index']) + with h5py.File(outfile, 'r+') as spe: + spe_root = list(spe.keys())[0] + en0 = spe[f'{spe_root}/instrument/fermi/energy'][()] + if 'fermi' in to_copy: + del spe[f'{spe_root}/instrument/fermi'] + for grp in to_copy: + src = raw[f'/raw_data_1/instrument/{grp}'] + h5py.Group.copy(src, src, spe[f'{spe_root}/instrument/']) + if 'fermi' in to_copy: + spe[f'{spe_root}/instrument/fermi/energy'][()] = en0 + detroot = f'{spe_root}/instrument/detector_elements_1' + spe.create_group(detroot) + spe[detroot].attrs['NX_class'] = np.array('NXdetector', dtype='S') + for df0, df1 in zip(['SPEC', 'UDET', 'DELT', 'LEN2', 'CODE', 'TTHE', 'UT01'], \ + ['det2spec', 'detector_number', 'delt', 'distance', 'detector_code', 'polar_angle', 'azimuthal_angle']): + src = raw[f'/raw_data_1/isis_vms_compat/{df0}'] + h5py.Group.copy(src, src, spe[detroot], df1) + for nn in range(raw['/raw_data_1/isis_vms_compat/NUSE'][0]): + src = raw[f'/raw_data_1/isis_vms_compat/UT{nn+1:02d}'] + h5py.Group.copy(src, src, spe[detroot], f'user_table{nn+1:02d}') + spec2work = f'{spe_root}/instrument/detector_elements_1/spec2work' + ws = mtd[in_ws] + if n_spec == ws.getNumberHistograms(): + s2w = np.arange(n_spec) + else: + nmon = np.array(raw['/raw_data_1/isis_vms_compat/NMON'])[0] + spec = np.array(raw['/raw_data_1/isis_vms_compat/SPEC'])[nmon:] + udet = np.array(raw['/raw_data_1/isis_vms_compat/UDET'])[nmon:] + _, iq = np.unique(spec, return_index=True) + s2 = np.hstack([np.array(ws.getSpectrum(ii).getDetectorIDs()) for ii in range(ws.getNumberHistograms())]) + _, c1, _ = np.intersect1d(udet[iq], s2, return_indices=True) + s2w = -np.ones(iq.shape, dtype=np.int32) + s2w[c1] = np.array(ws.getIndicesFromDetectorIDs(s2[c1].tolist())) + spe[detroot].create_dataset('spec2work', s2w.shape, dtype='i4', data=s2w) + +#======================================================== +# MARI specific functions + +def remove_extra_spectra_if_mari(wsname='ws'): + ws = mtd[wsname] + if ws.getNumberHistograms() == 919: + RemoveSpectra(ws, [0], OutputWorkspace=wsname) + try: + ws_mon = ws.getMonitorWorkspace() + except RuntimeError: + try: + ws_mon = mtd[f'{wsname}_monitors'] + except KeyError: + ws_mon = None + if ws_mon: + if ws_mon.getNumberHistograms() > 3: + RemoveSpectra(ws_mon, [3], OutputWorkspace=ws_mon.name()) + +def shift_frame_for_mari_lowE(origEi, wsname='ws_norm', wsmon='ws_monitors'): + ws_norm, ws_monitors = (mtd[wsname], mtd[wsmon]) + ws_out, wsmon_out = lhs_info('names') + if origEi < 4.01: + # If Ei < 4, mon 3 is in 2nd frame so need to shift it by 20ms + ws_monitors = ScaleX(wsmon, 20000, Operation='Add', IndexMin=2, IndexMax=2, OutputWorkspace=wsmon_out) + if origEi < 3.1: + # Additionally if Ei<3.1, data will also be in 2nd frame, shift all ToF by 20ms + ws_norm = ScaleX(wsname, 20000, Operation='Add', IndexMin=0, IndexMax=ws_norm.getNumberHistograms()-1, OutputWorkspace=ws_out) + return ws_norm, ws_monitors + +def gen_ana_bkg(quietws='MAR28952', target_ws=None): + # Generates an analytic background workspace from the quiet counts data + # by fitting each spectra with a decaying exponential a*exp(-b*ToF) + if quietws not in mtd: + Load(Filename=f'{quietws}.nxs', OutputWorkspace=quietws) + remove_extra_spectra_if_mari(quietws) + NormaliseByCurrent(InputWorkspace=quietws, OutputWorkspace=quietws) + bw, bx = (100, 1) + ws = mtd[quietws] + if target_ws: + wx = RebinToWorkspace(ws, target_ws, PreserveEvents=False, OutputWorkspace='bkg_his') + else: + wx = Rebin(ws, f'1700,{bx},19000', PreserveEvents=False, OutputWorkspace='bkg_his') + ws = Rebin(ws, f'1700,{bw},19000', OutputWorkspace='ws_ana_tmp') + current = ws.run().getProtonCharge() + xx = ws.extractX(); xx = (xx[:, 1:] + xx[:, :-1]) / 2. + yy = ws.extractY() + def fitfu(x, *pars): + return pars[0] * np.exp(-x * pars[1]) + fdi = {} + for spn in range(yy.shape[0]): + if np.max(yy[spn,:]) > (0.002): + popt, _ = scipy.optimize.curve_fit(fitfu, xx[spn,:], yy[spn,:], p0=[6e-4, 1./4000.]) + popt[0] /= bw + fdi[spn] = popt + x1 = wx.readX(spn) + y1 = fitfu((x1[1:] + x1[:-1]) / 2., *popt) + wx.setY(spn, y1) + wx.setE(spn, np.sqrt(y1*bw*current) / (bw*current)) + else: + wx.setY(spn, wx.readY(spn)*0) + wx.setE(spn, wx.readY(spn)*0) + if '28952' in quietws: + for spn in range(697 - 4, 759 - 4): + bb = spn - 259 + if bb in fdi.keys(): + x1 = wx.readX(bb) + y1 = fitfu((x1[1:] + x1[:-1]) / 2., *fdi[bb]) + wx.setY(spn, y1) + wx.setE(spn, np.sqrt(y1*bw*current) / (bw*current)) + bkg_ev = ConvertToEventWorkspace(wx) + return bkg_ev, wx + +#======================================================== +# Auto-Ei routine + + +def mode(inp): + return max(set(inp), key=list(inp).count) + + +def roundlog10(val): + expn = int(min(0, np.floor(np.log10(val)) - 2)) + return np.round(val, decimals=-expn) + + +def autoei(ws): + assert hasattr(ws, 'getInstrument') and hasattr(ws, 'getRun'), 'Input must be a Mantid workspace' + inst = ws.getInstrument().getName() + run = ws.getRun() + + def getLog(logname): + logv = run.getProperty(logname) + if 'Filtered' in str(type(logv)): + return logv.filtered_value + t0 = run.startTime().to_datetime64() + t1 = run.endTime().to_datetime64() + t = logv.times + mask = np.logical_and(t0 <= t, t <= t1) + return logv.value[mask] + + + def getfracLog(logname, frac=0.25, expn=1.): + val = getLog(logname) + if expn == 1.: + return mode(val[int(len(val) * frac):]) + else: + return mode(np.round(val[int(len(val) * frac):] * expn)) / expn + + + if inst == 'LET': + try: + c1_freq = getfracLog('Chopper1_Disk1_speed') + except ValueError: + return [] + c5_mode = run.getProperty('Chopper5_slits').value[-1] + if 'open' in c5_mode.lower(): # Whitebeam mode + return [] + f1, f2, c4_freq = (getfracLog(nm) for nm in ['Chopper5_Disk1_speed', 'Chopper3_speed', 'Chopper4_speed']) + ph1a, ph1b, ph2, ph3, ph4, ph5a, ph5b = (getfracLog(nm) for nm in ['Chopper1_Disk1_phase', 'Chopper1_Disk2_phase', + 'Chopper2_phase', 'Chopper3_phase', 'Chopper4_phase', 'Chopper5_Disk1_phase', 'Chopper5_Disk2_phase']) + # Moderator Chopper distances in m + lm1a, lm1b, lm2, lm3, lm4, lm5a, lm5b = (7.839, 7.828, 8.4, 11.749, 15.664, 23.504, 23.496) + # Chopper time offsets in us + tc1a, tc1b, tc3, tc4 = ((slp / frq) + yc for slp, frq, yc in zip([1.3333e5, 0.7694e5, 186944, 115000], \ + [c1_freq, c1_freq, f2, c4_freq], [7.3, 2.4, -87.2, -82.607])) + if 'resolution' in c5_mode.lower(): + tc5a = (1.5631e6 / f1) + 0.70 + tc5b = ( 1.3149e6 / f1) + 3.13 + elif 'flux' in c5_mode.lower(): + tc5a = (1.4520e6 / f1) + 2.56 + tc5b = (1.1204e6 / f1) + 3.33 + elif 'intermediate' in c5_mode.lower(): + tc5a = (1.0630e6 / f1) + 5.88 + tc5b = (1.8149e6 / f1) + 2.17 + else: + warnings.warning(f'Unknown mode "{c5_mode}" using High Flux settings') + tc5a = (1.4520e6 / f1) + 2.56 + tc5b = (1.1204e6 / f1) + 3.33 + # Determines the allowed reps through each set of choppers + tfoc = 70. # LET sets tfoc=80 for Ei<5 and 60 otherwise, we don't know which it is so use the average + periods = [1.e6 / nslit / frq for nslit, frq in zip([6, 2, 6, 2], [c1_freq, f2, c4_freq, f1])] + delays = [phase + tc - tfoc for phase, tc in zip([ph1a, ph3, ph4, ph5a], [tc1a, tc3, tc4, tc5a])] + # Chopper 2 is the bandwidth chopper and its phase determines the min and max ToF. + tf0 = 2500 if ph2 > 50000 else ph2 # Sets a minimum such that Ei~=50 is max Ei + tf1 = (ph2 + 27000) % 100000 # Opening width of chopper 2 is ~26ms [TODO: tune this value!] + def inrange(tf, l): + lim = [2286.26 * l / sqrte for sqrte in [2286.26 * lm2 / t for t in [tf0, tf1]]] + return tf >= lim[0] and tf < lim[1] + # Determine the reps which make it through the other choppers 1, 3, 4 and 5 + eis, eisv = ([], []) + for l, d, p, in zip([lm1a, lm3, lm4, lm5a], delays, periods): + eisv.append(np.array([(((2286.26 * l) / tf)**2) for tf in [(d + s*p) for s in range(-30, 30)] if inrange(tf, l)])) + # Calculates which reps which passes Chopper 5 also make it through all the others within 5% tolerance + for ei in eisv[3]: + istrans = [np.any((np.abs(eisn - ei) / eisn) < 0.05) for eisn in eisv[:3]] + if all(istrans): + eis.append(roundlog10(ei)) + return eis + + + elif inst == 'MARI': + try: + freq = mode(getLog('Fermi_Speed')) + except ValueError: + return [] + phase1, phase2 = (mode(getLog(nam)) for nam in ['Phase_Thick_1', 'Phase_Thick_2']) + delay = getfracLog('Fermi_delay') + lmc = 10.04 # Moderator-Fermi distance + period = 1.e6 / freq + try: + ei_nominal = mode(getLog('Ei_nominal')) + except RuntimeError: # Old file + ei_nominal = ((2286.26 * lmc) / delay)**2 + sqrt_ei = np.sqrt(ei_nominal) + delay_calc = ((2286.26 * lmc) / sqrt_ei) + t_offset_ref = {'S':2033.3/freq-5.4, 'G':1339.9/freq-7.3} + t_offset = delay - (delay_calc % period) + chopper_type = min(t_offset_ref.keys(), key=lambda x:np.abs(t_offset - t_offset_ref[x])) + nom_disk1, nom_disk2 = (((2286.26 * l) / sqrt_ei) - c for l, c in zip([7.861, 7.904], [5879., 6041.])) + delt_disk1, delt_disk2 = (ph - nom for ph, nom in zip([phase1, phase2], [nom_disk1, nom_disk2])) + disk_delta = delt_disk2 - delt_disk1 + slots_delta = np.round(disk_delta / 202.11) / 10 + assert slots_delta % 1.0 < 0.2, 'Bad slots calculation' + slots = {0:[0,1,2,4], 1:[0,1], 2:[0,2], 4:[0]}[abs(int(slots_delta))] + disk_ref = 6 - (np.round(delt_disk1 / 202.11) / 10) + assert disk_ref % 1.0 < 0.2, f'Bad disk calculation' + disk = {0:disk_ref, 1:disk_ref-1, 2:1 if disk_ref==2 else 0, 4:0}[abs(int(slots_delta))] + reps = [d-disk for d in slots] + eis_disk = {((2286.26*lmc) / (delay_calc + s*2500.))**2 for s in reps} + period = period / 2. if 'G' in chopper_type.upper() else period + eis = {((2286.26*lmc) / (delay_calc + s*period))**2 for s in range(-10, 10)} + inrange = lambda x: (x > 1 and x < 2.9) or (x > 4 and x < 1000) + return [roundlog10(ei) for ei in np.sort(list(eis.intersection(eis_disk)))[::-1] if inrange(ei)] + + + elif inst == 'MAPS': + try: + freq = mode(getLog('Fermi_Speed')) + except ValueError: + return [] + period = 1.e6 / freq + delay_angle = getfracLog('Fermi_Delay', expn=10.) + disk_delay = mode(getLog('Disc_Delay')) + lmc = 10.143 # Mod-Fermi distance + ldc = 8.831 # Mod-Disk distance + pickup_position = 180. + t_offsets = [20000 * (slit - pickup_position) / 360. for slit in [180., -39.1, 0., 39.1]] + disk_eis = [((2286.372 * ldc) / tof)**2 for tof in [disk_delay + t_offset for t_offset in t_offsets] if tof > 0] + fermi_eis = [] + for pickup_pos, t_off in zip([195.1, 14.429], [26.2, 24.53]): # 'S', and 'A' chopper respectively + delay = period * (delay_angle - pickup_pos) / 360. + fermi_eis += [((2286.372 * lmc) / tf)**2 for tf in [(delay + s*period) for s in range(-10, 10)] if tf > 100 and tf < 12000] + fermi_eis = np.array(fermi_eis) + eis = [] + for ei in disk_eis: + ed = np.abs(fermi_eis - ei) / fermi_eis + if np.any(ed < 0.1): + eis.append(roundlog10(ei)) + return eis + + + elif inst == 'MERLIN': + run = ws.getRun() + try: + freq = mode(getLog('Chopper_Speed')) + except ValueError: + return [] + delay = getfracLog('Chopper_Delay') + disk_delay = mode(getLog('Disc_Delay')) + rrm_mode = np.abs(disk_delay - 13700) < 10 or np.abs(disk_delay - 12400) < 10 + lmc = 9.995 # Mod-Fermi distance + ldc = 9.2176 # Mod-Disk distance + if not rrm_mode: + return [roundlog10(((2286.26*ldc) / disk_delay)**2)] + else: # Assume using Gd chopper (t_offset = 2000/freq-5, for 'S' it is 300/freq-8) + period = 20000 * (25. / freq) + tof = (2286.26 * lmc) / np.sqrt(roundlog10(((2286.26 * lmc) / (delay - 2000/freq - 2))**2)) + tfmx = 7500 if np.abs(disk_delay - 12400) < 10 else 8500 + return [roundlog10(((2286.26*lmc) / tf)**2) for tf in [(tof + s*period) for s in range(-10, 10)] if tf > 1500 and tf < tfmx] + + + else: + raise RuntimeError(f'Instrument {inst} not supported') + + +#======================================================== +# Iliad driver routines + +_DGRED = None + +class DG_reduction_wrapper: + + subsdict = {'\nconfig':'\n#config', + 'save_dir = ':'save_dir = None #', + 'INSTRUMENT_NAME':'MARI', + 'MASK_FILE_XML':'mari_mask2023_1.xml', + 'RINGS_MAP_XML':'mari_res2013.map', + 'whitevan\s*=\s*[0-9]*':'whitevan = 28580', + 'sample\s*=\s*\\[*[\\]0-9,]+':'sample = [28727, 28728]', + 'sample_bg\s*=\s*\\[*[\\]0-9,]+':'sample_bg = None', + 'wv_file\s*=\s*[\\\'A-z0-9\\.]*':'wv_file = \'WV_28580.txt\'', + 'wv_detrange\s*=\s*[\\[\\]0-9,]*':'wv_detrange = None', + 'Ei_list\s*=\s*[\\[\\]\\.0-9,]+.*':'Ei_list = [1.84, 1.1]'} + + def __init__(self): + self.curdir = abspath(dirname(__file__)) + self.reduction = self.load_code(os.path.join(self.curdir, 'DG_reduction.py')) + self.whitevan = self.load_code(os.path.join(self.curdir, 'DG_whitevan.py')) + self.monovan = self.load_code(os.path.join(self.curdir, 'DG_monovan.py')) + + def load_code(self, filename): + loader = importlib.machinery.SourceFileLoader('', filename) + # Loads the main reduction script and compiles it to bytecode + src = loader.get_data(filename).decode() + if 'INSTRUMENT_NAME' in src: + for ky, val in self.subsdict.items(): + src = re.sub(ky, val, src) + x0, x1 = (src.find('#!begin_params'), src.find('#!end_params')) + src_param = src[x0:x1] + src_body = (src[:x0] + re.sub('\n', '\n#', src_param) + src[x1:]).encode() + code = loader.source_to_code(src_body, filename) + # Parses the variables section + tmp_mod = types.ModuleType('dgredwrapper_param') + exec(src_param, tmp_mod.__dict__) + params = {k:getattr(tmp_mod, k) for k in dir(tmp_mod) if not k.startswith('_')} + return code, params + + def __call__(self, mod='reduction', **kwargs): + code, params = getattr(self, mod) + env = params + tmp_mod = types.ModuleType('DG_red_exec') + tmp_mod.__file__ = os.path.join(self.curdir, f'DG_{mod}.py') + env.update(tmp_mod.__dict__) + env.update(kwargs) + exec(code, env) + +def run_reduction(mod='reduction', **kwargs): + global _DGRED + if _DGRED is None: + _DGRED = DG_reduction_wrapper() + if 'inst' in kwargs: + config['default.instrument'] = kwargs['inst'] + _DGRED(mod=mod, **kwargs) + +def run_whitevan(**kwargs): + run_reduction(mod='whitevan', **kwargs) + +def run_monovan(**kwargs): + run_reduction(mod='monovan', **kwargs) + +def iliad(runno, ei, wbvan, monovan=None, sam_mass=None, sam_rmm=None, sum_runs=False, **kwargs): + wv_name = wbvan if (isinstance(wbvan, (str, int, float)) or len(wbvan)==1) else wbvan[0] + wv_file = f'WV_{wv_name}.txt' + wv_args = {} + if 'inst' in kwargs: + config['default.instrument'] = kwargs['inst'] + wv_args = {'inst':kwargs['inst']} + try: + LoadAscii(wv_file, OutputWorkspace=wv_file) + ReplaceSpecialValues(wv_file, SmallNumberThreshold=1e-20, SmallNumberValue='NaN', OutputWorkspace=wv_file) + except ValueError: + run_whitevan(whitevan=wbvan, **wv_args) + Ei_list = ei if hasattr(ei, '__iter__') else [ei] + if 'hard_mask_file' in kwargs: + kwargs['mask'] = kwargs.pop('hard_mask_file') + if 'sumruns' not in kwargs and sum_runs is True: + kwargs['sumruns'] = True + if monovan is not None and isinstance(monovan, (int, float)) and monovan > 0: + mv_file = f'MV_{monovan}.txt' + mvkw = {'mask':kwargs['mask']} if 'mask' in kwargs else {} + if 'inst' in kwargs: mvkw['inst'] = kwargs['inst'] + try: + LoadAscii(wv_file, OutputWorkspace=wv_file) + except ValueError: + run_monovan(monovan=monovan, Ei_list=Ei_list, wv_file=wv_file, **mvkw) + kwargs['sample_mass'] = sam_mass + kwargs['sample_fwt'] = sam_rmm + run_reduction(sample=runno, Ei_list=ei if hasattr(ei, '__iter__') else [ei], wv_file=wv_file, **kwargs) diff --git a/direct_inelastic/LET/USER_Files_description.xml b/direct_inelastic/LET/USER_Files_description.xml index c470ee5..d1028c6 100644 --- a/direct_inelastic/LET/USER_Files_description.xml +++ b/direct_inelastic/LET/USER_Files_description.xml @@ -17,6 +17,7 @@ + + diff --git a/direct_inelastic/MARI/DG_reduction.py b/direct_inelastic/MARI/DG_reduction.py new file mode 100644 index 0000000..32216b0 --- /dev/null +++ b/direct_inelastic/MARI/DG_reduction.py @@ -0,0 +1,516 @@ +# ======================================================= +# +# EXCITATIONS INSTRUMENTS REDUCTION SCRIPT +# JRS & MDL 3/7/23 +# +# Reads sample run(s), corrects for backgound, +# normalises to a white vanadium run, and the +# absolute normalisation factor for each Ei (if MV +# file is specified). +# - Masks bad detectors with hard mask only. +# - Converts time-of-flight to energy transfer +# - performs Q-rebinning for QENS data ('_red.nxs' format) +# - Outputs .nxspe, .nxs or _red.nxs files +# +#======================================================== + +# import mantid algorithms, numpy and matplotlib +from mantid.simpleapi import * +from mantid.api import AnalysisDataService as ADS +from mantid.kernel.funcinspect import lhs_info +import numpy as np +import os, sys +import time +from importlib import reload + +#!begin_params +#=======================User Inputs======================= +powder = True # powder or 1to1 map +sumruns = False # set to True to sum sample runs +sample = [93338] # sample runs (list) +sample_bg = 93329 # single background run +sample_cd = None # Cadmium run for instrument background subtraction +wv_file = 'WV_91329.txt' # white vanadium integral file (mandatory) +Ei_list = ['auto'] # incident energies (Ei) - from PyChop +Erange = [-0.8,0.0025,0.8] # energy transfer range to output in fractions of Ei +trans = [0.95,0.95,0.95] # elastic line transmission factors for each Ei +mask = 'MASK_FILE_XML' # hard mask +# what to do if see run with same angle as previous run (and sumruns==False) +# - 'replace': sum and replace first nxspe file created +# - 'ignore': ignore and create an nxspe for each run +same_angle_action = 'ignore' +# to enable creating multiple reduced data files from a single +# "continuous scan" set the following variables to a valid log "block" +# name and the bin size (width) to something larger than zero. +# An optional unit can be given. +# If sumruns=True, the sample list is summed and then divided +# If sumruns=False and multiple sample runs are listed, an error is raised. +# This is also only compatible with same_angle_action='ignore' +cs_block = None +cs_block_unit = '' +cs_bin_size = 0 +# MARI only: used to subtract an analytic approx of instrument background +sub_ana = False +#======================================================== + +#==================Absolute Units Inputs================= +mv_file = None # pre-processed MV calibration +sample_mass = 1 # mass of the sample +sample_fwt = 50.9415 # formula weight of sample +monovan_mass = None # Mass of vanadium sample (ask local contact) +#======================================================== + +#==================Local Contact Inputs================== +inst = 'INSTRUMENT_NAME' +cycle = 'CYCLE_ID' # cycle number +fixei = True # True for LET since no monitor 3 +powdermap = 'RINGS_MAP_XML' # rings mapping file - must be .xml format +file_wait = 30 # wait for data file to appear (seconds) +keepworkspaces = False # should be false for Horace scans +saveformat = '.nxspe' # format of output, '.nxspe', '.nxs' +QENS = False # output Q-binned data for QENS data analysis '_red.nxs' +Qbins = 20 # approximate number of Q-bins for QENS +theta_range = [5., 65.] # useful theta range for Q-binning (QENS) +idebug = False # keep workspaces and check absolute units on elastic line +save_dir = f'/instrument/{inst}/RBNumber/USER_RB_FOLDER' # Set to None to avoid reseting +psi_motor_name = 'rot' # name of the rotation motor in the logs +angles_workspace = 'angles_ws' # name of workspace to store previously seen angles +sumruns_savemem = False # Compresses event in summed ws to save memory + # (causes some loss of data so cannot use event filtering) +cs_smidge = 0.001 # Tolerance on continuous scan bins size +#======================================================== +#!end_params + +# ==============================setup directroties================================ +config['default.instrument'] = inst +if save_dir is not None: + config['defaultsave.directory'] = save_dir +if inst == 'MARI': + source = 'Moderator' + m2spec = 2 # specID of monitor2 (pre-sample) + m3spec = 3 # specID of monitor3 (post-sample) + #monovan_mass = 32.58 # mass of vanadium cylinder + same_angle_action = 'ignore' +elif inst == 'MERLIN': + source = 'undulator' + m2spec = 69636 # specID of monitor2 (pre-sample) + m3spec = 69640 # specID of monitor3 (post-sample) + #monovan_mass = 32.62 # mass of vanadium cylinder +elif inst == 'MAPS': + source = 'undulator' + m2spec = 36867 # specID of monitor2 (pre-sample) + m3spec = 36868 # specID of monitor3 (post-sample) + #monovan_mass = 30.1 # mass of vanadium cylinder +elif inst == 'LET': + source = 'undulator' + m2spec = 98310 # specID of monitor2 (pre-sample) + m3spec = None # specID of monitor3 (post-sample) + #monovan_mass = 3.56 # mass of vanadium cylinder +else: + raise RuntimeError(f'Unrecognised instrument: {inst}') + +cycle_shortform = cycle[2:] if cycle.startswith('20') else cycle +datadir = f'/archive/NDX{inst}/Instrument/data/cycle_{cycle_shortform}/' +instfiles = '/usr/local/mprogs/InstrumentFiles/' +mapdir = f'{instfiles}/{inst.swapcase()}/' +config.appendDataSearchDir(mapdir) +config.appendDataSearchDir(datadir) +ws_list = ADS.getObjectNames() + +# Try to load subsidiary utilities +sys.path.append(os.path.dirname(__file__)) +try: + from reduction_utils import * + utils_loaded = True +except ImportError: + utils_loaded = False + def rename_existing_ws(*a, **k): raise RuntimeError('Unable to load utils') + def remove_extra_spectra_if_mari(*a, **k): raise RuntimeError('Unable to load utils') + def autoei(*a, **k): + raise RuntimeError('Cannot use Auto-Ei: ' \ + 'You need to copy the reduction_utils.py file somewhere on the Python path') + +#======================================================== +# Helper functions +def tryload(irun): # loops till data file appears + if isinstance(irun, str) and irun in mtd: + # If users give a string which matches an existing workspace, use that instead + return rename_existing_ws(irun) + while True: + try: + ws = Load(str(irun),LoadMonitors=True) + except TypeError: + ws = Load(str(irun),LoadMonitors='Separate') + except ValueError: + print(f'...waiting for run #{irun}') + Pause(file_wait) + continue + break + if inst == 'MARI': + remove_extra_spectra_if_mari('ws') + if sub_ana is True and not sample_cd: + if 'bkg_ev' not in mtd: + gen_ana_bkg() + current = mtd['ws'].run().getProtonCharge() + if isinstance(mtd['ws'], mantid.dataobjects.EventWorkspace): + ws = mtd['ws'] - mtd['bkg_ev'] * current + else: + try: + ws = mtd['ws'] - mtd['bkg_his'] * current + except ValueError: + gen_ana_bkg(target_ws=mtd['ws']) + ws = mtd['ws'] - mtd['bkg_his'] * current + return mtd['ws'] + +def load_sum(run_list, block_name=None): + for ii, irun in enumerate(run_list): + tryload(irun) + if ii == 0: + w_buf = CloneWorkspace('ws') + w_buf_monitors = CloneWorkspace('ws_monitors') + if block_name: + bval = mtd['ws'].getRun().getLogData(block_name).filtered_value + print(f'run #{irun} loaded') + else: + w_buf = Plus('w_buf', 'ws') + w_buf_monitors = Plus('w_buf_monitors', 'ws_monitors') + if block_name: + bv2 = mtd['ws'].getRun().getLogData(block_name).filtered_value + bval = np.concatenate((bval, bv2)) + print(f'run #{irun} added') + ADS.remove('ws') + ADS.remove('ws_monitors') + wsout_name = lhs_info('names')[0] + if wsout_name != 'w_buf': + RenameWorkspace('w_buf_monitors', wsout_name+'_monitors') + RenameWorkspace('w_buf', wsout_name) + return (mtd[wsout_name], bval) if block_name else mtd[wsout_name] + +#======================================================== + + +print(f'\n======= {inst} data reduction =======') +print(f'Working directory... {ConfigService.Instance().getString("defaultsave.directory")}\n') + +# ============================create lists if necessary========================== +if isinstance(sample, str) or not hasattr(sample, '__iter__'): + sample = [sample] +if sample_bg is not None and (isinstance(sample_bg, str) or not hasattr(sample_bg, '__iter__')): + sample_bg = [sample_bg] +if sample_cd is not None and (isinstance(sample_cd, str) or not hasattr(sample_cd, '__iter__')): + sample_cd = [sample_cd] + +#==================================load hard mask================================ +if mask is None: + print(f'{inst}: WARNING - No hard mask! Bad detectors wont be removed') +if mask not in ws_list and mask is not None: + print(f'{inst}: Loading hard mask - {mask}') + LoadMask(inst,mask,OutputWorkspace=mask) +else: + print(f'{inst}: Using previously loaded hard mask - {mask}') + +# ===============================load whitevan file============================== +if wv_file is None: + raise RuntimeError(f'{inst}: ERROR - white vanadium calibration file missing') +if wv_file not in ws_list: + print(f'{inst}: Loading white vanadium - {wv_file}') + LoadAscii(Filename=wv_file, OutputWorkspace=wv_file) + ReplaceSpecialValues(wv_file, SmallNumberThreshold=1e-20, SmallNumberValue='NaN', OutputWorkspace=wv_file) +else: + print(f'{inst}: Using previously loaded white vanadium - {wv_file}') + +# =======================load background runs and sum========================= +if sample_cd is not None: + ws_cd = load_sum(sample_cd) +if sample_bg is not None: + ws_bg = load_sum(sample_bg) + if sample_cd is not None: + ws_bg = ws_bg - ws_cd + ws_bg = NormaliseByCurrent(ws_bg) + +# ==========================continuous scan stuff============================= +if cs_block and cs_bin_size > 0: + if len(sample) > 1 and not sumruns: + raise RuntimeError('Continous scans for multiple (non-summed) files not supported') + if same_angle_action.lower() != 'ignore': + raise RuntimeError(f'Continous scans not compatible with same_angle_action="{same_angle_action}"') + # Sets sumruns false so don't redo the sum below. Instead sum here to get filtered block values. + # This is due to a bug in the Mantid "Plus" algorithm, which gives incorrect filtered values + # https://github.com/mantidproject/mantid/issues/36194 + sumruns = False + ws_full, bval = load_sum(sample, cs_block) + ws_monitors = CloneWorkspace('ws_full_monitors') + bval_range = max(bval) - min(bval) + bval_nbins = int(bval_range / cs_bin_size) + bval_remainder = bval_range - cs_bin_size * bval_nbins + # Overrides the sample list in order to hijack the loop over run numbers below + irun_orig = sample[0] + sample = [x*cs_bin_size + min(bval) for x in range(bval_nbins)] + unit = cs_block_unit + if not unit: + unit = ws_full.getRun().getLogData(cs_block).units + print(f'{cs_block} = {min(bval):.1f} {unit} to {max(bval):.1f} {unit}') + print(f'... filtering in {cs_bin_size:.1f} {unit} steps') + print(f'... N={bval_nbins} bins with {bval_remainder:.2f} {unit} remainder') + +# =======================sum sample runs if required========================= +sumsuf = sumruns and len(sample) > 1 +if sumruns: + ws = load_sum(sample) + sample = [sample[0]] + +# =============================auto-Ei stuff================================= +is_auto = lambda x: isinstance(x, str) and 'auto' in x.lower() +if is_auto(Ei_list) or hasattr(Ei_list, '__iter__') and is_auto(Ei_list[0]): + use_auto_ei = True + try: + Ei_list = autoei(ws) + except NameError: + fn = str(sample[0]) + if not fn.startswith(inst[:3]): fn = f'{inst[:3]}{fn}' + if fn.endswith('.raw'): fn = fn[:-4] + if fn[-4:-2] == '.s': fn = fn[:-4] + '.n0' + fn[-2:] + if not fn.endswith('.nxs') and fn[-5:-2] != '.n0': fn += '.nxs' + Ei_list = autoei(LoadNexusMonitors(fn, OutputWorkspace='ws_tmp_mons')) + print(f"Automatically determined Ei's: {Ei_list}") + if len(trans) < len(Ei_list): + print(f'{inst}: WARNING - not enough transmision values for auto-Ei. ' \ + 'Extending list with last (end) value') + trans += [trans[-1]]*(len(Ei_list) - len(trans)) +else: + use_auto_ei = False + +# ============================load monovan file============================== +mv_fac = [] +if mv_file is not None and monovan_mass is not None: + use_auto_ei = False + if mv_file not in ws_list: + print(f'{inst}: Loading monovan calibration factors - {mv_file}') + LoadAscii(mv_file,OutputWorkspace=mv_file) + mv_ws = mtd[mv_file] + mv_eis = mv_ws.readX(0) + mv_cal = mv_ws.readY(0) +# check that monovan is compatible with Ei_list) + Ei_diff = sum([x-y for x,y in zip(mv_eis,Ei_list)]) + if (abs(Ei_diff) > 0): + raise RuntimeError('----ERROR: Monovan file Eis not compatible with Ei_list') + for Ei in Ei_list: + ii = np.where(mv_eis == Ei) + mvf = mv_cal[ii][0] + van_fwt = 50.9415 # vanadium molar mass (g) + van_xs = 5080. / 4. / np.pi # vanadium differential cross-section (mbarns/st) + mvf = (monovan_mass / van_fwt) / (sample_mass / sample_fwt) * van_xs / mvf + mv_fac.append(mvf) +else: + print(f'{inst}: Skipping absolute calibration') + mv_fac = [x/x for x in Ei_list] # monovan factors = 1 by default + +# =====================angles cache stuff==================================== +if utils_loaded: + build_angles_ws(sample, angles_workspace, psi_motor_name) + runs_with_angles_already_seen = [] +else: + same_angle_action = 'ignore' + +# =======================sample run loop===================================== +for irun in sample: + t = time.time() # start the clock + + # check loaded workspaces and remove old nxspe workspaces + if not keepworkspaces: + ws_list = ADS.getObjectNames() + nx_list = [ss for ss in ws_list if saveformat in ss] + for ss in nx_list: + ADS.remove(ss) + + # For continuous scans generate the workspace based on binned log values + # NB. we replaced the sample (run) list with list of block min bin boundaries + if cs_block and cs_bin_size > 0: + val = round(irun + cs_bin_size/2,1) + print(f"Filtering: {cs_block}={val:.1f} ± {round(cs_bin_size/2,2):.2f} {unit}") + ws = FilterByLogValue(ws_full, cs_block, irun, irun + cs_bin_size - cs_smidge) + + # if this run is at an angle already seen and we are doing 'replace', skip + if same_angle_action.lower() != 'ignore' and irun in runs_with_angles_already_seen: + continue + + print('============') + if not sumruns and not cs_block: + # Checks rotation angle + if same_angle_action.lower() != 'ignore': + runs_with_same_angles = get_angle(irun, angles_workspace, psi_motor_name, tryload) + if len(runs_with_same_angles) > 1: + ws = load_sum(runs_with_same_angles) + irun = runs_with_same_angles[0] + runs_with_angles_already_seen += runs_with_same_angles + else: + tryload(irun) + print(f'Loading run# {irun}') + + ws = NormaliseByCurrent('ws') + if sumruns and sumruns_savemem: + ws = CompressEvents(ws, Tolerance=1e-5) # Tolerance in microseconds + + # instrument geometry to work out ToF ranges + sampos = ws.getInstrument().getSample().getPos() + l1 = (sampos - ws.getInstrument().getSource().getPos()).norm() + l2 = (ws.getDetector(0).getPos() - sampos).norm() + + # Updates ei_loop if auto _and_ not using monovan + if use_auto_ei: + Ei_list = autoei(ws) + if len(trans) < len(Ei_list): trans += [trans[-1]]*(len(Ei_list) - len(trans)) + if len(mv_fac) < len(Ei_list): mv_fac += [1]*(len(Ei_list) - len(mv_fac)) + +# ============================= Ei loop ===================================== + for ienergy in range(len(Ei_list)): + Ei = Ei_list[ienergy] + origEi = Ei + tr = trans[ienergy] + mvf = mv_fac[ienergy] + print(f'\n{inst}: Reducing data for Ei={Ei:.2f} meV') + + # Low energy reps on MARI end up in the second frame + t_shift = 20000 if inst == 'MARI' and origEi < 3.1 else 0 + + tofs = ws.readX(0) + tof_min = np.sqrt(l1**2 * 5.227e6 / Ei) - t_shift + tof_max = tof_min + np.sqrt(l2**2 * 5.226e6 / (Ei*(1-Erange[-1]))) + ws_rep = CropWorkspace(ws, max(min(tofs), tof_min), min(max(tofs), tof_max)) + + if sample_cd is not None: + print(f'... subtracting Cd background') + ws_rep = ws_rep - ws_cd + + if sample_bg is not None: + print(f'... subtracting background - transmission factor = {tr:.2f}') + ws_rep = ws_rep/tr - ws_bg + else: + ws_rep = Scale('ws_rep', 1 / tr, 'Multiply') + + # normalise to WB vanadium and apply fixed mask + print('... normalising/masking data') + ws_rep = Divide('ws_rep', wv_file) # white beam normalisation + MaskDetectors(ws_rep, MaskedWorkspace=mask, ForceInstrumentMasking=True) + + # t2e section + print('... t2e section') + ws_monitors = mtd['ws_monitors'] + spectra = ws_monitors.getSpectrumNumbers() + index = spectra.index(m2spec) + m2pos = ws.detectorInfo().position(index)[2] + + if inst == 'MARI' and utils_loaded and origEi < 4.01: + # Shifts data / monitors into second frame for MARI + ws_rep, ws_monitors = shift_frame_for_mari_lowE(origEi, wsname='ws_rep', wsmon='ws_monitors') + + # this section shifts the time-of-flight such that the monitor2 peak + # in the current monitor workspace (post monochromator) is at t=0 and L=0 + # note that the offest is not the predicted value due to energy dependence of the source position + + if m3spec is not None: + (Ei,mon2_peak,_,_) = GetEi(ws_monitors,Monitor1Spec=m2spec,Monitor2Spec=m3spec,EnergyEstimate=Ei,FixEi=fixei) + print(f'... refined Ei={Ei:.2f} meV') + else: + (Ei,mon2_peak,_,_) = GetEi(ws_monitors,Monitor2Spec=m2spec,EnergyEstimate=Ei,FixEi=fixei) + + print(f'... m2 tof={mon2_peak:.2f} mus, m2 pos={m2pos:.2f} m') + + ws_rep = ScaleX(ws_rep, Factor=-mon2_peak, Operation='Add', InstrumentParameter='DelayTime', Combine=True) + MoveInstrumentComponent(ws_rep, ComponentName=source, Z=m2pos, RelativePosition=False) + + ws_rep = ConvertUnits(ws_rep, 'DeltaE', EMode='Direct', EFixed=Ei) + ws_out = Rebin(ws_rep, [x*origEi for x in Erange], PreserveEvents=False) + ws_out = DetectorEfficiencyCor(ws_out, IncidentEnergy=Ei) + ws_out = CorrectKiKf(ws_out, Efixed=Ei, EMode='Direct') + ADS.remove('ws_rep') + + # monovan scaling + if mv_file is not None: + print(f'... applying mono van calibration factor {mvf:.1f}') + ws_out = Scale('ws_out', mvf, 'Multiply') + + # Sets output file name + if cs_block and cs_bin_size > 0: + ofile_prefix = f'{inst[:3]}{irun_orig}' + ofile_suffix = f"_{val:.1f}{unit}" + else: + ofile_prefix = f'{inst[:3]}{irun}' + ofile_suffix = '' + + # rings grouping if desired + if powder or inst == 'MARI' or QENS: + ws_out=GroupDetectors(ws_out, MapFile=powdermap, Behaviour='Average') + ofile_suffix += '_powder' + if inst == 'MARI' or QENS: + ofile_suffix = '' + print(f'... powder grouping using {powdermap}') + else: + ofile_suffix += '_1to1' + if inst == 'MARI': + if sumsuf: + ofile_suffix += 'sum' + if sample_bg is not None: + ofile_suffix += '_sub' + if sub_ana: + ofile_suffix += '_sa' + + # output nxspe file + ofile = f'{ofile_prefix}_{origEi:g}meV{ofile_suffix}' + + # check elastic line (debug mode) + if idebug: + Rebin('ws_out', [-0.05*Ei, 100, 0.05*Ei], PreserveEvents=False, OutputWorkspace=ofile+'_elastic') + Transpose(f'{ofile}_elastic', OutputWorkspace=f'{ofile}_elastic') + + # output appropriate formats + ws_out = ConvertToDistribution(ws_out) + if QENS: + saveformat = '.nxs' + print(f'{inst}: Writing {ofile}{saveformat}') + if saveformat.lower() == '.nxspe': + SaveNXSPE('ws_out', ofile+saveformat, Efixed=Ei, KiOverKfScaling=True) + if utils_loaded: + copy_inst_info(ofile+saveformat, 'ws_out') # Copies instrument info for Horace + elif saveformat.lower() == '.nxs': + rmlogs = {'events_log', 'frame_log', 'good_frame_log', 'period_log', 'proton_charge', 'raw_events_log'} + RemoveLogs('ws_out', KeepLogs=','.join(set(mtd['ws_out'].run().keys()).difference(rmlogs))) + SaveNexusProcessed('ws_out', ofile+saveformat, PreserveEvents=False, CompressNexus=True) + if QENS: + print('... outputting QENS "_red" format') + theta = np.array([theta_range[0], theta_range[1]])*np.pi/180. + Q = 1.39 * np.sqrt(Ei) * np.sin(theta) + Q = np.around(Q*10) / 10. + Qbin = int((Q[1] - Q[0])) / Qbins + print(f'... Q-axis = [{Q[0]+Qbin:g},{Qbin:g},{Q[1]-Qbin:g}]') + ws_out = SofQW3('ws_out', [Q[0]+Qbin, Qbin, Q[1]-Qbin], 'Direct', Efixed=Ei) + # these lines are Anthony Lim's method of removing NaNs + # NaN are changed to zeros, and then placed at the end of the energy range + spectra = range(ws_out.getNumberHistograms()) + for ispec in spectra: + x = ws_out.dataX(ispec) + y = ws_out.dataY(ispec) + e = ws_out.dataE(ispec) + smidge = 1e-6 + for iq in range(len(x)-1): + if np.isnan(y[iq]): + x[iq] = np.max(x) + smidge + y[iq] = 0.0 + e[iq] = 0.0 + smidge += 1e-6 + SaveNexus('ws_out', ofile+"_red"+saveformat) + + CloneWorkspace('ws_out',OutputWorkspace=ofile+saveformat) + +# ============================= End of Ei loop ================================ + + print(f'\n{inst}: Reduction complete in {time.time() - t:.1f} seconds\n') + +# ============================= End of run loop ================================ + +# cleanup +if not idebug: + ws_list = ADS.getObjectNames() + nx_list = [ss for ss in ws_list if 'w_buf' in ss or 'ws' in ss] + for ss in nx_list: + ADS.remove(ss) diff --git a/direct_inelastic/MARI/USER_Files_description.xml b/direct_inelastic/MARI/USER_Files_description.xml index baffd0f..d28b0ca 100644 --- a/direct_inelastic/MARI/USER_Files_description.xml +++ b/direct_inelastic/MARI/USER_Files_description.xml @@ -13,13 +13,17 @@ not full rb folder path (/home/wkc26243/RB1501020), The values of these variables are taken from archive for current cycle and user --> - - + + + + - + diff --git a/direct_inelastic/MERLIN/USER_Files_description.xml b/direct_inelastic/MERLIN/USER_Files_description.xml index 1e2502b..39fed7e 100644 --- a/direct_inelastic/MERLIN/USER_Files_description.xml +++ b/direct_inelastic/MERLIN/USER_Files_description.xml @@ -17,6 +17,7 @@ +