diff --git a/direct_inelastic/ISIS/DG_monovan.py b/direct_inelastic/ISIS/DG_monovan.py index 49c1885..7545ee9 100644 --- a/direct_inelastic/ISIS/DG_monovan.py +++ b/direct_inelastic/ISIS/DG_monovan.py @@ -15,7 +15,6 @@ # import mantid algorithms, numpy and matplotlib from mantid.simpleapi import * from mantid.api import AnalysisDataService as ADS -import matplotlib.pyplot as plt import numpy as np import time diff --git a/direct_inelastic/ISIS/DG_reduction.py b/direct_inelastic/ISIS/DG_reduction.py index 5fa9ce7..f8aebab 100644 --- a/direct_inelastic/ISIS/DG_reduction.py +++ b/direct_inelastic/ISIS/DG_reduction.py @@ -1,25 +1,24 @@ # ======================================================= # # EXCITATIONS INSTRUMENTS REDUCTION SCRIPT -# JRS 20/6/23 +# 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). s +# file is specified). # - Masks bad detectors with hard mask only. # - Converts time-of-flight to energy transfer -# - performs Q-rebinning for QENS data only -# - Outputs .nxspe, .nxs or _red.nxs fiile +# - 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 -import matplotlib.pyplot as plt import numpy as np -import re, os, h5py +import os, sys import time from importlib import reload @@ -33,12 +32,18 @@ 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' #======================================================== #==================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================== @@ -49,65 +54,13 @@ file_wait = 30 # wait for data file to appear (seconds) keepworkspaces = True # should be false for Horace scans saveformat = '.nxspe' # format of output, '.nxspe', '.nxs' -QENS = True # output Q-binned data for QENS data analysis '_red.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 -#======================================================== - -#======================================================== -def tryload(irun): # loops till data file appears - 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}') - time.sleep(file_wait) - continue - break -#======================================================== - -#======================================================== -# 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(get_raw_file_from_ws(mtd[in_ws])) - if not os.path.exists(outfile): - outfile = mantid.simpleapi.config['defaultsave.directory'] + '/' + os.path.basename(outfile) - print(outfile) - with h5py.File(get_raw_file_from_ws(mtd[in_ws]), 'r') as raw: - with h5py.File(outfile, 'r+') as spe: - 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])] - print(spe.keys()) - 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: - print(grp) - 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) - for df0, df1 in zip(['SPEC', 'UDET', 'DELT', 'LEN2', 'CODE', 'TTHE', 'UT01'], \ - ['spectrum_number', '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}') +psi_motor_name = 'rot' # name of the rotation motor in the logs +angles_workspace = 'angles_ws' # name of workspace to store previously seen angles #======================================================== # ==============================setup directroties================================ @@ -118,39 +71,85 @@ def copy_inst_info(outfile, in_ws): source = 'Moderator' m2spec = 2 # specID of monitor2 (pre-sample) m3spec = 3 # specID of monitor3 (post-sample) - monovan_mass = 32.58 # mass of vanadium cylinder + #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 + #monovan_mass = 32.62 # mass of vanadium cylinder elif inst == 'MAPS': source = 'undulator' - m2spec = 41475 # specID of monitor2 (pre-sample) - m3spec = 41476 # specID of monitor3 (post-sample) - monovan_mass = 30.1 # mass of vanadium cylinder + 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 + #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}/' -mapdir = f'/usr/local/mprogs/InstrumentFiles/{inst.swapcase()}/' +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(instfiles) +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') + +#======================================================== +# 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') + return mtd['ws'] + +def load_sum(run_list): + for ii, irun in enumerate(run_list): + tryload(irun) + if ii == 0: + w_buf = CloneWorkspace('ws') + w_buf_monitors = CloneWorkspace('ws_monitors') + print(f'run #{irun} loaded') + else: + w_buf = Plus('w_buf', 'ws') + w_buf_monitors = Plus('w_buf_monitors', 'ws_monitors') + print(f'run #{irun} added') +#======================================================== + + print(f'\n======= {inst} data reduction =======') print(f'Working directory... {ConfigService.Instance().getString("defaultsave.directory")}\n') # ============================create lists if necessary========================== -if type(sample) is not list: +if isinstance(sample, str) or not hasattr(sample, '__iter__'): sample = [sample] -if type(sample_bg) is not list and sample_bg is not None: +if sample_bg is not None and (isinstance(sample_bg, str) or not hasattr(sample_bg, '__iter__')): sample_bg = [sample_bg] #==================================load hard mask================================ @@ -163,6 +162,8 @@ def copy_inst_info(outfile, in_ws): 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) @@ -172,7 +173,7 @@ def copy_inst_info(outfile, in_ws): # ===============================load monovan file=============================== mv_fac = [] -if mv_file is not None: +if mv_file is not None and monovan_mass is not None: if mv_file not in ws_list: print(f'{inst}: Loading monovan calibration factors - {mv_file}') LoadAscii(mv_file,OutputWorkspace=mv_file) @@ -196,33 +197,24 @@ def copy_inst_info(outfile, in_ws): # =======================load background runs and sum========================= if sample_bg is not None: - for irun in sample_bg: - ws_bg = Load(str(irun),LoadMonitors=False) - if (irun == sample_bg[0]): - print(f'background run #{irun} loaded') - w_buf = CloneWorkspace('ws_bg') - else: - print(f'background run #{irun} added') - w_buf = Plus('w_buf', 'ws_bg') + load_sum(sample_bg) ws_bg = CloneWorkspace('w_buf') ws_bg = NormaliseByCurrent('ws_bg') # =======================sum sample runs if required========================= if sumruns: - for irun in sample: - tryload(irun) - if (irun == sample[0]): - print(f'run #{irun} loaded') - w_buf = CloneWorkspace('ws') - w_buf_monitors = CloneWorkspace('ws_monitors') - else: - print(f'run #{irun} added') - w_buf = Plus('w_buf', 'ws') - w_buf_monitors = Plus('w_buf_monitors', 'ws_monitors') + 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) + runs_with_angles_already_seen = [] +else: + same_angle_action = 'ignore' + # =======================sample run loop===================================== for irun in sample: t = time.time() # start the clock @@ -234,12 +226,23 @@ def copy_inst_info(outfile, in_ws): for ss in nx_list: ADS.remove(ss) + # 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: - tryload(irun) - print(f'Loading run# {irun}') - if inst == 'MARI' and mtd['ws'].getNumberHistograms() > 918: - ws = RemoveSpectra('ws',[0]) + # 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 + else: + tryload(irun) + print(f'Loading run# {irun}') ws = NormaliseByCurrent('ws') # ============================= Ei loop ===================================== @@ -267,6 +270,10 @@ def copy_inst_info(outfile, in_ws): 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 @@ -297,11 +304,11 @@ def copy_inst_info(outfile, in_ws): if powder or inst == 'MARI' or QENS: ws_out=GroupDetectors(ws_out, MapFile=powdermap, Behaviour='Average') ofile_suffix = '_powder' - if inst == 'MARI': + if inst == 'MARI' or QENS: ofile_suffix = '' - if QENS: - ofile_suffix = '_red' print(f'... powder grouping using {powdermap}') + if sample_bg is not None and inst == 'MARI': + ofile_suffix += '_sub' # output nxspe file ofile = f'{inst[:3]}{irun}_{origEi:g}meV{ofile_suffix}' @@ -318,12 +325,15 @@ def copy_inst_info(outfile, in_ws): print(f'{inst}: Writing {ofile}{saveformat}') if saveformat.lower() == '.nxspe': SaveNXSPE('ws_out', ofile+saveformat, Efixed=Ei, KiOverKfScaling=True) - copy_inst_info(ofile, 'ws_out') - elif saveformat.lower() == '.nxs' and not QENS: + #if utils_loaded and not powder: + # 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) if QENS: print('... outputting QENS "_red" format') - theta = np.array([7.5,65.])*np.pi/180. + 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 @@ -343,7 +353,7 @@ def copy_inst_info(outfile, in_ws): y[iq] = 0.0 e[iq] = 0.0 smidge += 1e-6 - SaveNexus('ws_out', ofile+saveformat) + SaveNexus('ws_out', ofile+"_red"+saveformat) CloneWorkspace('ws_out',OutputWorkspace=ofile+saveformat) @@ -354,7 +364,7 @@ def copy_inst_info(outfile, in_ws): # ============================= End of run loop ================================ # cleanup -if (not idebug): +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: diff --git a/direct_inelastic/ISIS/DG_whitevan.py b/direct_inelastic/ISIS/DG_whitevan.py index 048bf52..c4455ee 100644 --- a/direct_inelastic/ISIS/DG_whitevan.py +++ b/direct_inelastic/ISIS/DG_whitevan.py @@ -14,7 +14,7 @@ # import mantid algorithms, numpy and matplotlib from mantid.simpleapi import * from mantid.api import AnalysisDataService as ADS -import matplotlib.pyplot as plt +from mantid.kernel.funcinspect import lhs_info import numpy as np import time @@ -24,7 +24,6 @@ whitevan = 78715 # white vanadium run whitevan_bg = None # background for the white vanadium whitevan_trans = 1 # transmission factor -mask_file = 'MASK_FILE_XML' # hardmask file (str or None) #======================================================== #==================Local Contact Inputs================== @@ -33,7 +32,6 @@ wv_lrange = [1,5] # wavelength integration limits for output of vanadium integrals wv_detrange = [30000,60000] # spectrum index range for average intensity calculation idebug = False # keep itermediate workspaces for debugging -mask_dir = None # folder for mask file (str or None) save_dir = f'/instrument/{inst}/RBNumber/USER_RB_FOLDER' # Set to None to avoid resetting #======================================================== @@ -43,18 +41,49 @@ cycle_shortform = cycle[2:] if cycle.startswith('20') else cycle data_dir = f'/archive/NDX{inst}/Instrument/data/cycle_{cycle_shortform}/' config.appendDataSearchDir(data_dir) -config.appendDataSearchDir(save_dir) +if save_dir is not None: + config.appendDataSearchDir(save_dir) + +#======================================================== +def try_load_no_mon(irun, ws_name=None): + # LoadISISNexus needs LoadMonitors='Exclude' + # LoadEventNexus needs LoadMonitors=False + if ws_name is None: + nret, ws_name = lhs_info('both') + if nret > 1: + raise RuntimeError('load: Only one output workspace name can be given') + ws_name = ws_name[0] + try: + return Load(str(irun), LoadMonitors='Exclude', OutputWorkspace=ws_name) + except ValueError: + return Load(str(irun), LoadMonitors=False, OutputWorkspace=ws_name) + +def load_sum(run_list): + if isinstance(run_list, str) or not hasattr(run_list, '__iter__'): + run_list = [run_list] + for ii, irun in enumerate(run_list): + try_load_no_mon(irun, 'ws') + if ii == 0: + w_buf = CloneWorkspace('ws') + print(f'run #{irun} loaded') + else: + w_buf = Plus('w_buf', 'ws') + print(f'run #{irun} added') + nret, ws_name = lhs_info('both') + if nret == 1: + CloneWorkspace('w_buf', OutputWorkspace=ws_name[0]) +#======================================================== print(f'\n======= {inst} white van reduction =======') # =================load white van and background and subtract===================== -print(f'WHITE_VAN {inst}: Loading white vanadium run# {whitevan}') -wv = Load(str(whitevan), LoadMonitors='Exclude') +print(f'WHITE_VAN {inst}: Loading white vanadium run(s)') +wv = load_sum(whitevan) wv = NormaliseByCurrent('wv') wv_corrected = Scale('wv', 1/whitevan_trans, 'Multiply') if (whitevan_bg is not None): - print(f'... subtracting white vanadium background run# {whitevan_bg}') - wv_bg = Load(str(whitevan_bg), LoadMonitors='Exclude') + print(f'... subtracting white vanadium background run(s)') + wv_bg = load_sum(whitevan_bg) wv_bg = NormaliseByCurrent('wv_bg') wv_corrected = wv/whitevan_trans - wv_bg @@ -62,12 +91,7 @@ print('... Calculating white vanadium integrals') print(f'... Summing between {wv_lrange[0]:.1f} < \u03BB < {wv_lrange[1]:.1f} \u212B') WV_normalised_integrals = ConvertUnits(wv_corrected, 'Wavelength') -WV_normalised_integrals = Rebin(WV_normalised_integrals, f'{wv_lrange[0]},100,{wv_lrange[1]}') -if mask_file is not None: - if mask_dir is not None: - config.appendDataSearchDir(mask_dir) - LoadMask(inst, mask, OutputWorkspace=mask) - MaskDetectors(wv_normalised_integrals, MaskedWorkspace=mask) +WV_normalised_integrals = Rebin(WV_normalised_integrals, f'{wv_lrange[0]},100,{wv_lrange[1]}', PreserveEvents=False) wv_normt = Transpose(WV_normalised_integrals) if wv_detrange is not None: wv_normt = CropWorkspace(wv_normt, XMin=wv_detrange[0], Xmax=wv_detrange[1]) @@ -76,7 +100,8 @@ WV_normalised_integrals = Scale(WV_normalised_integrals, scale_factor, 'Multiply') # ===========================output integrals file================================ -ofile = f'WV_{whitevan}.txt' +wv_name = whitevan if (isinstance(whitevan, str) 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}') SaveAscii(WV_normalised_integrals, ofile)